    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.dfp;
    019
    020    import java.util.Arrays;
    021
    022    import org.apache.commons.math3.FieldElement;
    023
    024    /**
    025     *  Decimal floating point library for Java
    026     *
    027     *  <p>Another floating point class.  This one is built using radix 10000
    028     *  which is 10<sup>4</sup>, so its almost decimal.</p>
    029     *
    030     *  <p>The design goals here are:
    031     *  <ol>
    032     *    <li>Decimal math, or close to it</li>
    033     *    <li>Settable precision (but no mix between numbers using different settings)</li>
    034     *    <li>Portability.  Code should be kept as portable as possible.</li>
    035     *    <li>Performance</li>
    036     *    <li>Accuracy  - Results should always be +/- 1 ULP for basic
    037     *         algebraic operation</li>
    038     *    <li>Comply with IEEE 854-1987 as much as possible.
    039     *         (See IEEE 854-1987 notes below)</li>
    040     *  </ol></p>
    041     *
    042     *  <p>Trade offs:
    043     *  <ol>
    044     *    <li>Memory foot print.  I'm using more memory than necessary to
    045     *         represent numbers to get better performance.</li>
    046     *    <li>Digits are bigger, so rounding is a greater loss.  So, if you
    047     *         really need 12 decimal digits, better use 4 base 10000 digits
    048     *         there can be one partially filled.</li>
    049     *  </ol></p>
    050     *
    051     *  <p>Numbers are represented  in the following form:
    052     *  <pre>
    053     *  n  =  sign &times; mant &times; (radix)<sup>exp</sup>;</p>
    054     *  </pre>
    055     *  where sign is &plusmn;1, mantissa represents a fractional number between
    056     *  zero and one.  mant[0] is the least significant digit.
    057     *  exp is in the range of -32767 to 32768</p>
    058     *
    059     *  <p>IEEE 854-1987  Notes and differences</p>
    060     *
    061     *  <p>IEEE 854 requires the radix to be either 2 or 10.  The radix here is
    062     *  10000, so that requirement is not met, but  it is possible that a
    063     *  subclassed can be made to make it behave as a radix 10
    064     *  number.  It is my opinion that if it looks and behaves as a radix
    065     *  10 number then it is one and that requirement would be met.</p>
    066     *
    067     *  <p>The radix of 10000 was chosen because it should be faster to operate
    068     *  on 4 decimal digits at once instead of one at a time.  Radix 10 behavior
    069     *  can be realized by adding an additional rounding step to ensure that
    070     *  the number of decimal digits represented is constant.</p>
    071     *
    072     *  <p>The IEEE standard specifically leaves out internal data encoding,
    073     *  so it is reasonable to conclude that such a subclass of this radix
    074     *  10000 system is merely an encoding of a radix 10 system.</p>
    075     *
    076     *  <p>IEEE 854 also specifies the existence of "sub-normal" numbers.  This
    077     *  class does not contain any such entities.  The most significant radix
    078     *  10000 digit is always non-zero.  Instead, we support "gradual underflow"
    079     *  by raising the underflow flag for numbers less with exponent less than
    080     *  expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
    081     *  Thus the smallest number we can represent would be:
    082     *  1E(-(MIN_EXP-digits-1)*4),  eg, for digits=5, MIN_EXP=-32767, that would
    083     *  be 1e-131092.</p>
    084     *
    085     *  <p>IEEE 854 defines that the implied radix point lies just to the right
    086     *  of the most significant digit and to the left of the remaining digits.
    087     *  This implementation puts the implied radix point to the left of all
    088     *  digits including the most significant one.  The most significant digit
    089     *  here is the one just to the right of the radix point.  This is a fine
    090     *  detail and is really only a matter of definition.  Any side effects of
    091     *  this can be rendered invisible by a subclass.</p>
    092     * @see DfpField
    093     * @version $Id: Dfp.java 1416643 2012-12-03 19:37:14Z tn $
    094     * @since 2.2
    095     */
    096    public class Dfp implements FieldElement<Dfp> {
    097
    098        /** The radix, or base of this system.  Set to 10000 */
    099        public static final int RADIX = 10000;
    100
    101        /** The minimum exponent before underflow is signaled.  Flush to zero
    102         *  occurs at minExp-DIGITS */
    103        public static final int MIN_EXP = -32767;
    104
    105        /** The maximum exponent before overflow is signaled and results flushed
    106         *  to infinity */
    107        public static final int MAX_EXP =  32768;
    108
    109        /** The amount under/overflows are scaled by before going to trap handler */
    110        public static final int ERR_SCALE = 32760;
    111
    112        /** Indicator value for normal finite numbers. */
    113        public static final byte FINITE = 0;
    114
    115        /** Indicator value for Infinity. */
    116        public static final byte INFINITE = 1;
    117
    118        /** Indicator value for signaling NaN. */
    119        public static final byte SNAN = 2;
    120
    121        /** Indicator value for quiet NaN. */
    122        public static final byte QNAN = 3;
    123
    124        /** String for NaN representation. */
    125        private static final String NAN_STRING = "NaN";
    126
    127        /** String for positive infinity representation. */
    128        private static final String POS_INFINITY_STRING = "Infinity";
    129
    130        /** String for negative infinity representation. */
    131        private static final String NEG_INFINITY_STRING = "-Infinity";
    132
    133        /** Name for traps triggered by addition. */
    134        private static final String ADD_TRAP = "add";
    135
    136        /** Name for traps triggered by multiplication. */
    137        private static final String MULTIPLY_TRAP = "multiply";
    138
    139        /** Name for traps triggered by division. */
    140        private static final String DIVIDE_TRAP = "divide";
    141
    142        /** Name for traps triggered by square root. */
    143        private static final String SQRT_TRAP = "sqrt";
    144
    145        /** Name for traps triggered by alignment. */
    146        private static final String ALIGN_TRAP = "align";
    147
    148        /** Name for traps triggered by truncation. */
    149        private static final String TRUNC_TRAP = "trunc";
    150
    151        /** Name for traps triggered by nextAfter. */
    152        private static final String NEXT_AFTER_TRAP = "nextAfter";
    153
    154        /** Name for traps triggered by lessThan. */
    155        private static final String LESS_THAN_TRAP = "lessThan";
    156
    157        /** Name for traps triggered by greaterThan. */
    158        private static final String GREATER_THAN_TRAP = "greaterThan";
    159
    160        /** Name for traps triggered by newInstance. */
    161        private static final String NEW_INSTANCE_TRAP = "newInstance";
    162
    163        /** Mantissa. */
    164        protected int[] mant;
    165
    166        /** Sign bit: 1 for positive, -1 for negative. */
    167        protected byte sign;
    168
    169        /** Exponent. */
    170        protected int exp;
    171
    172        /** Indicator for non-finite / non-number values. */
    173        protected byte nans;
    174
    175        /** Factory building similar Dfp's. */
    176        private final DfpField field;
    177
    178        /** Makes an instance with a value of zero.
    179         * @param field field to which this instance belongs
    180         */
    181        protected Dfp(final DfpField field) {
    182            mant = new int[field.getRadixDigits()];
    183            sign = 1;
    184            exp = 0;
    185            nans = FINITE;
    186            this.field = field;
    187        }
    188
    189        /** Create an instance from a byte value.
    190         * @param field field to which this instance belongs
    191         * @param x value to convert to an instance
    192         */
    193        protected Dfp(final DfpField field, byte x) {
    194            this(field, (long) x);
    195        }
    196
    197        /** Create an instance from an int value.
    198         * @param field field to which this instance belongs
    199         * @param x value to convert to an instance
    200         */
    201        protected Dfp(final DfpField field, int x) {
    202            this(field, (long) x);
    203        }
    204
    205        /** Create an instance from a long value.
    206         * @param field field to which this instance belongs
    207         * @param x value to convert to an instance
    208         */
    209        protected Dfp(final DfpField field, long x) {
    210
    211            // initialize as if 0
    212            mant = new int[field.getRadixDigits()];
    213            nans = FINITE;
    214            this.field = field;
    215
    216            boolean isLongMin = false;
    217            if (x == Long.MIN_VALUE) {
    218                // special case for Long.MIN_VALUE (-9223372036854775808)
    219                // we must shift it before taking its absolute value
    220                isLongMin = true;
    221                ++x;
    222            }
    223
    224            // set the sign
    225            if (x < 0) {
    226                sign = -1;
    227                x = -x;
    228            } else {
    229                sign = 1;
    230            }
    231
    232            exp = 0;
    233            while (x != 0) {
    234                System.arraycopy(mant, mant.length - exp, mant, mant.length - 1 - exp, exp);
    235                mant[mant.length - 1] = (int) (x % RADIX);
    236                x /= RADIX;
    237                exp++;
    238            }
    239
    240            if (isLongMin) {
    241                // remove the shift added for Long.MIN_VALUE
    242                // we know in this case that fixing the last digit is sufficient
    243                for (int i = 0; i < mant.length - 1; i++) {
    244                    if (mant[i] != 0) {
    245                        mant[i]++;
    246                        break;
    247                    }
    248                }
    249            }
    250        }
    251
    252        /** Create an instance from a double value.
    253         * @param field field to which this instance belongs
    254         * @param x value to convert to an instance
    255         */
    256        protected Dfp(final DfpField field, double x) {
    257
    258            // initialize as if 0
    259            mant = new int[field.getRadixDigits()];
    260            sign = 1;
    261            exp = 0;
    262            nans = FINITE;
    263            this.field = field;
    264
    265            long bits = Double.doubleToLongBits(x);
    266            long mantissa = bits & 0x000fffffffffffffL;
    267            int exponent = (int) ((bits & 0x7ff0000000000000L) >> 52) - 1023;
    268
    269            if (exponent == -1023) {
    270                // Zero or sub-normal
    271                if (x == 0) {
    272                    // make sure 0 has the right sign
    273                    if ((bits & 0x8000000000000000L) != 0) {
    274                        sign = -1;
    275                    }
    276                    return;
    277                }
    278
    279                exponent++;
    280
    281                // Normalize the subnormal number
    282                while ( (mantissa & 0x0010000000000000L) == 0) {
    283                    exponent--;
    284                    mantissa <<= 1;
    285                }
    286                mantissa &= 0x000fffffffffffffL;
    287            }
    288
    289            if (exponent == 1024) {
    290                // infinity or NAN
    291                if (x != x) {
    292                    sign = (byte) 1;
    293                    nans = QNAN;
    294                } else if (x < 0) {
    295                    sign = (byte) -1;
    296                    nans = INFINITE;
    297                } else {
    298                    sign = (byte) 1;
    299                    nans = INFINITE;
    300                }
    301                return;
    302            }
    303
    304            Dfp xdfp = new Dfp(field, mantissa);
    305            xdfp = xdfp.divide(new Dfp(field, 4503599627370496l)).add(field.getOne());  // Divide by 2^52, then add one
    306            xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));
    307
    308            if ((bits & 0x8000000000000000L) != 0) {
    309                xdfp = xdfp.negate();
    310            }
    311
    312            System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
    313            sign = xdfp.sign;
    314            exp  = xdfp.exp;
    315            nans = xdfp.nans;
    316
    317        }
    318
    319        /** Copy constructor.
    320         * @param d instance to copy
    321         */
    322        public Dfp(final Dfp d) {
    323            mant  = d.mant.clone();
    324            sign  = d.sign;
    325            exp   = d.exp;
    326            nans  = d.nans;
    327            field = d.field;
    328        }
    329
    330        /** Create an instance from a String representation.
    331         * @param field field to which this instance belongs
    332         * @param s string representation of the instance
    333         */
    334        protected Dfp(final DfpField field, final String s) {
    335
    336            // initialize as if 0
    337            mant = new int[field.getRadixDigits()];
    338            sign = 1;
    339            exp = 0;
    340            nans = FINITE;
    341            this.field = field;
    342
    343            boolean decimalFound = false;
    344            final int rsize = 4;   // size of radix in decimal digits
    345            final int offset = 4;  // Starting offset into Striped
    346            final char[] striped = new char[getRadixDigits() * rsize + offset * 2];
    347
    348            // Check some special cases
    349            if (s.equals(POS_INFINITY_STRING)) {
    350                sign = (byte) 1;
    351                nans = INFINITE;
    352                return;
    353            }
    354
    355            if (s.equals(NEG_INFINITY_STRING)) {
    356                sign = (byte) -1;
    357                nans = INFINITE;
    358                return;
    359            }
    360
    361            if (s.equals(NAN_STRING)) {
    362                sign = (byte) 1;
    363                nans = QNAN;
    364                return;
    365            }
    366
    367            // Check for scientific notation
    368            int p = s.indexOf("e");
    369            if (p == -1) { // try upper case?
    370                p = s.indexOf("E");
    371            }
    372
    373            final String fpdecimal;
    374            int sciexp = 0;
    375            if (p != -1) {
    376                // scientific notation
    377                fpdecimal = s.substring(0, p);
    378                String fpexp = s.substring(p+1);
    379                boolean negative = false;
    380
    381                for (int i=0; i<fpexp.length(); i++)
    382                {
    383                    if (fpexp.charAt(i) == '-')
    384                    {
    385                        negative = true;
    386                        continue;
    387                    }
    388                    if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') {
    389                        sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
    390                    }
    391                }
    392
    393                if (negative) {
    394                    sciexp = -sciexp;
    395                }
    396            } else {
    397                // normal case
    398                fpdecimal = s;
    399            }
    400
    401            // If there is a minus sign in the number then it is negative
    402            if (fpdecimal.indexOf("-") !=  -1) {
    403                sign = -1;
    404            }
    405
    406            // First off, find all of the leading zeros, trailing zeros, and significant digits
    407            p = 0;
    408
    409            // Move p to first significant digit
    410            int decimalPos = 0;
    411            for (;;) {
    412                if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
    413                    break;
    414                }
    415
    416                if (decimalFound && fpdecimal.charAt(p) == '0') {
    417                    decimalPos--;
    418                }
    419
    420                if (fpdecimal.charAt(p) == '.') {
    421                    decimalFound = true;
    422                }
    423
    424                p++;
    425
    426                if (p == fpdecimal.length()) {
    427                    break;
    428                }
    429            }
    430
    431            // Copy the string onto Stripped
    432            int q = offset;
    433            striped[0] = '0';
    434            striped[1] = '0';
    435            striped[2] = '0';
    436            striped[3] = '0';
    437            int significantDigits=0;
    438            for(;;) {
    439                if (p == (fpdecimal.length())) {
    440                    break;
    441                }
    442
    443                // Don't want to run pass the end of the array
    444                if (q == mant.length*rsize+offset+1) {
    445                    break;
    446                }
    447
    448                if (fpdecimal.charAt(p) == '.') {
    449                    decimalFound = true;
    450                    decimalPos = significantDigits;
    451                    p++;
    452                    continue;
    453                }
    454
    455                if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
    456                    p++;
    457                    continue;
    458                }
    459
    460                striped[q] = fpdecimal.charAt(p);
    461                q++;
    462                p++;
    463                significantDigits++;
    464            }
    465
    466
    467            // If the decimal point has been found then get rid of trailing zeros.
    468            if (decimalFound && q != offset) {
    469                for (;;) {
    470                    q--;
    471                    if (q == offset) {
    472                        break;
    473                    }
    474                    if (striped[q] == '0') {
    475                        significantDigits--;
    476                    } else {
    477                        break;
    478                    }
    479                }
    480            }
    481
    482            // special case of numbers like "0.00000"
    483            if (decimalFound && significantDigits == 0) {
    484                decimalPos = 0;
    485            }
    486
    487            // Implicit decimal point at end of number if not present
    488            if (!decimalFound) {
    489                decimalPos = q-offset;
    490            }
    491
    492            // Find the number of significant trailing zeros
    493            q = offset;  // set q to point to first sig digit
    494            p = significantDigits-1+offset;
    495
    496            while (p > q) {
    497                if (striped[p] != '0') {
    498                    break;
    499                }
    500                p--;
    501            }
    502
    503            // Make sure the decimal is on a mod 10000 boundary
    504            int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
    505            q -= i;
    506            decimalPos += i;
    507
    508            // Make the mantissa length right by adding zeros at the end if necessary
    509            while ((p - q) < (mant.length * rsize)) {
    510                for (i = 0; i < rsize; i++) {
    511                    striped[++p] = '0';
    512                }
    513            }
    514
    515            // Ok, now we know how many trailing zeros there are,
    516            // and where the least significant digit is
    517            for (i = mant.length - 1; i >= 0; i--) {
    518                mant[i] = (striped[q]   - '0') * 1000 +
    519                          (striped[q+1] - '0') * 100  +
    520                          (striped[q+2] - '0') * 10   +
    521                          (striped[q+3] - '0');
    522                q += 4;
    523            }
    524
    525
    526            exp = (decimalPos+sciexp) / rsize;
    527
    528            if (q < striped.length) {
    529                // Is there possible another digit?
    530                round((striped[q] - '0')*1000);
    531            }
    532
    533        }
    534
    535        /** Creates an instance with a non-finite value.
    536         * @param field field to which this instance belongs
    537         * @param sign sign of the Dfp to create
    538         * @param nans code of the value, must be one of {@link #INFINITE},
    539         * {@link #SNAN},  {@link #QNAN}
    540         */
    541        protected Dfp(final DfpField field, final byte sign, final byte nans) {
    542            this.field = field;
    543            this.mant    = new int[field.getRadixDigits()];
    544            this.sign    = sign;
    545            this.exp     = 0;
    546            this.nans    = nans;
    547        }
    548
    549        /** Create an instance with a value of 0.
    550         * Use this internally in preference to constructors to facilitate subclasses
    551         * @return a new instance with a value of 0
    552         */
    553        public Dfp newInstance() {
    554            return new Dfp(getField());
    555        }
    556
    557        /** Create an instance from a byte value.
    558         * @param x value to convert to an instance
    559         * @return a new instance with value x
    560         */
    561        public Dfp newInstance(final byte x) {
    562            return new Dfp(getField(), x);
    563        }
    564
    565        /** Create an instance from an int value.
    566         * @param x value to convert to an instance
    567         * @return a new instance with value x
    568         */
    569        public Dfp newInstance(final int x) {
    570            return new Dfp(getField(), x);
    571        }
    572
    573        /** Create an instance from a long value.
    574         * @param x value to convert to an instance
    575         * @return a new instance with value x
    576         */
    577        public Dfp newInstance(final long x) {
    578            return new Dfp(getField(), x);
    579        }
    580
    581        /** Create an instance from a double value.
    582         * @param x value to convert to an instance
    583         * @return a new instance with value x
    584         */
    585        public Dfp newInstance(final double x) {
    586            return new Dfp(getField(), x);
    587        }
    588
    589        /** Create an instance by copying an existing one.
    590         * Use this internally in preference to constructors to facilitate subclasses.
    591         * @param d instance to copy
    592         * @return a new instance with the same value as d
    593         */
    594        public Dfp newInstance(final Dfp d) {
    595
    596            // make sure we don't mix number with different precision
    597            if (field.getRadixDigits() != d.field.getRadixDigits()) {
    598                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    599                final Dfp result = newInstance(getZero());
    600                result.nans = QNAN;
    601                return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
    602            }
    603
    604            return new Dfp(d);
    605
    606        }
    607
    608        /** Create an instance from a String representation.
    609         * Use this internally in preference to constructors to facilitate subclasses.
    610         * @param s string representation of the instance
    611         * @return a new instance parsed from specified string
    612         */
    613        public Dfp newInstance(final String s) {
    614            return new Dfp(field, s);
    615        }
    616
    617        /** Creates an instance with a non-finite value.
    618         * @param sig sign of the Dfp to create
    619         * @param code code of the value, must be one of {@link #INFINITE},
    620         * {@link #SNAN},  {@link #QNAN}
    621         * @return a new instance with a non-finite value
    622         */
    623        public Dfp newInstance(final byte sig, final byte code) {
    624            return field.newDfp(sig, code);
    625        }
    626
    627        /** Get the {@link org.apache.commons.math3.Field Field} (really a {@link DfpField}) to which the instance belongs.
    628         * <p>
    629         * The field is linked to the number of digits and acts as a factory
    630         * for {@link Dfp} instances.
    631         * </p>
    632         * @return {@link org.apache.commons.math3.Field Field} (really a {@link DfpField}) to which the instance belongs
    633         */
    634        public DfpField getField() {
    635            return field;
    636        }
    637
    638        /** Get the number of radix digits of the instance.
    639         * @return number of radix digits
    640         */
    641        public int getRadixDigits() {
    642            return field.getRadixDigits();
    643        }
    644
    645        /** Get the constant 0.
    646         * @return a Dfp with value zero
    647         */
    648        public Dfp getZero() {
    649            return field.getZero();
    650        }
    651
    652        /** Get the constant 1.
    653         * @return a Dfp with value one
    654         */
    655        public Dfp getOne() {
    656            return field.getOne();
    657        }
    658
    659        /** Get the constant 2.
    660         * @return a Dfp with value two
    661         */
    662        public Dfp getTwo() {
    663            return field.getTwo();
    664        }
    665
    666        /** Shift the mantissa left, and adjust the exponent to compensate.
    667         */
    668        protected void shiftLeft() {
    669            for (int i = mant.length - 1; i > 0; i--) {
    670                mant[i] = mant[i-1];
    671            }
    672            mant[0] = 0;
    673            exp--;
    674        }
    675
    676        /* Note that shiftRight() does not call round() as that round() itself
    677         uses shiftRight() */
    678        /** Shift the mantissa right, and adjust the exponent to compensate.
    679         */
    680        protected void shiftRight() {
    681            for (int i = 0; i < mant.length - 1; i++) {
    682                mant[i] = mant[i+1];
    683            }
    684            mant[mant.length - 1] = 0;
    685            exp++;
    686        }
    687
    688        /** Make our exp equal to the supplied one, this may cause rounding.
    689         *  Also causes de-normalized numbers.  These numbers are generally
    690         *  dangerous because most routines assume normalized numbers.
    691         *  Align doesn't round, so it will return the last digit destroyed
    692         *  by shifting right.
    693         *  @param e desired exponent
    694         *  @return last digit destroyed by shifting right
    695         */
    696        protected int align(int e) {
    697            int lostdigit = 0;
    698            boolean inexact = false;
    699
    700            int diff = exp - e;
    701
    702            int adiff = diff;
    703            if (adiff < 0) {
    704                adiff = -adiff;
    705            }
    706
    707            if (diff == 0) {
    708                return 0;
    709            }
    710
    711            if (adiff > (mant.length + 1)) {
    712                // Special case
    713                Arrays.fill(mant, 0);
    714                exp = e;
    715
    716                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
    717                dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
    718
    719                return 0;
    720            }
    721
    722            for (int i = 0; i < adiff; i++) {
    723                if (diff < 0) {
    724                    /* Keep track of loss -- only signal inexact after losing 2 digits.
    725                     * the first lost digit is returned to add() and may be incorporated
    726                     * into the result.
    727                     */
    728                    if (lostdigit != 0) {
    729                        inexact = true;
    730                    }
    731
    732                    lostdigit = mant[0];
    733
    734                    shiftRight();
    735                } else {
    736                    shiftLeft();
    737                }
    738            }
    739
    740            if (inexact) {
    741                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
    742                dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
    743            }
    744
    745            return lostdigit;
    746
    747        }
    748
    749        /** Check if instance is less than x.
    750         * @param x number to check instance against
    751         * @return true if instance is less than x and neither are NaN, false otherwise
    752         */
    753        public boolean lessThan(final Dfp x) {
    754
    755            // make sure we don't mix number with different precision
    756            if (field.getRadixDigits() != x.field.getRadixDigits()) {
    757                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    758                final Dfp result = newInstance(getZero());
    759                result.nans = QNAN;
    760                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
    761                return false;
    762            }
    763
    764            /* if a nan is involved, signal invalid and return false */
    765            if (isNaN() || x.isNaN()) {
    766                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    767                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
    768                return false;
    769            }
    770
    771            return compare(this, x) < 0;
    772        }
    773
    774        /** Check if instance is greater than x.
    775         * @param x number to check instance against
    776         * @return true if instance is greater than x and neither are NaN, false otherwise
    777         */
    778        public boolean greaterThan(final Dfp x) {
    779
    780            // make sure we don't mix number with different precision
    781            if (field.getRadixDigits() != x.field.getRadixDigits()) {
    782                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    783                final Dfp result = newInstance(getZero());
    784                result.nans = QNAN;
    785                dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
    786                return false;
    787            }
    788
    789            /* if a nan is involved, signal invalid and return false */
    790            if (isNaN() || x.isNaN()) {
    791                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    792                dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
    793                return false;
    794            }
    795
    796            return compare(this, x) > 0;
    797        }
    798
    799        /** Check if instance is less than or equal to 0.
    800         * @return true if instance is not NaN and less than or equal to 0, false otherwise
    801         */
    802        public boolean negativeOrNull() {
    803
    804            if (isNaN()) {
    805                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    806                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
    807                return false;
    808            }
    809
    810            return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
    811
    812        }
    813
    814        /** Check if instance is strictly less than 0.
    815         * @return true if instance is not NaN and less than or equal to 0, false otherwise
    816         */
    817        public boolean strictlyNegative() {
    818
    819            if (isNaN()) {
    820                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    821                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
    822                return false;
    823            }
    824
    825            return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite());
    826
    827        }
    828
    829        /** Check if instance is greater than or equal to 0.
    830         * @return true if instance is not NaN and greater than or equal to 0, false otherwise
    831         */
    832        public boolean positiveOrNull() {
    833
    834            if (isNaN()) {
    835                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    836                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
    837                return false;
    838            }
    839
    840            return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
    841
    842        }
    843
    844        /** Check if instance is strictly greater than 0.
    845         * @return true if instance is not NaN and greater than or equal to 0, false otherwise
    846         */
    847        public boolean strictlyPositive() {
    848
    849            if (isNaN()) {
    850                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    851                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
    852                return false;
    853            }
    854
    855            return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite());
    856
    857        }
    858
    859        /** Get the absolute value of instance.
    860         * @return absolute value of instance
    861         */
    862        public Dfp abs() {
    863            Dfp result = newInstance(this);
    864            result.sign = 1;
    865            return result;
    866        }
    867
    868        /** Check if instance is infinite.
    869         * @return true if instance is infinite
    870         */
    871        public boolean isInfinite() {
    872            return nans == INFINITE;
    873        }
    874
    875        /** Check if instance is not a number.
    876         * @return true if instance is not a number
    877         */
    878        public boolean isNaN() {
    879            return (nans == QNAN) || (nans == SNAN);
    880        }
    881
    882        /** Check if instance is equal to zero.
    883         * @return true if instance is equal to zero
    884         */
    885        public boolean isZero() {
    886
    887            if (isNaN()) {
    888                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    889                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
    890                return false;
    891            }
    892
    893            return (mant[mant.length - 1] == 0) && !isInfinite();
    894
    895        }
    896
    897        /** Check if instance is equal to x.
    898         * @param other object to check instance against
    899         * @return true if instance is equal to x and neither are NaN, false otherwise
    900         */
    901        @Override
    902        public boolean equals(final Object other) {
    903
    904            if (other instanceof Dfp) {
    905                final Dfp x = (Dfp) other;
    906                if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
    907                    return false;
    908                }
    909
    910                return compare(this, x) == 0;
    911            }
    912
    913            return false;
    914
    915        }
    916
    917        /**
    918         * Gets a hashCode for the instance.
    919         * @return a hash code value for this object
    920         */
    921        @Override
    922        public int hashCode() {
    923            return 17 + (sign << 8) + (nans << 16) + exp + Arrays.hashCode(mant);
    924        }
    925
    926        /** Check if instance is not equal to x.
    927         * @param x number to check instance against
    928         * @return true if instance is not equal to x and neither are NaN, false otherwise
    929         */
    930        public boolean unequal(final Dfp x) {
    931            if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
    932                return false;
    933            }
    934
    935            return greaterThan(x) || lessThan(x);
    936        }
    937
    938        /** Compare two instances.
    939         * @param a first instance in comparison
    940         * @param b second instance in comparison
    941         * @return -1 if a<b, 1 if a>b and 0 if a==b
    942         *  Note this method does not properly handle NaNs or numbers with different precision.
    943         */
    944        private static int compare(final Dfp a, final Dfp b) {
    945            // Ignore the sign of zero
    946            if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
    947                a.nans == FINITE && b.nans == FINITE) {
    948                return 0;
    949            }
    950
    951            if (a.sign != b.sign) {
    952                if (a.sign == -1) {
    953                    return -1;
    954                } else {
    955                    return 1;
    956                }
    957            }
    958
    959            // deal with the infinities
    960            if (a.nans == INFINITE && b.nans == FINITE) {
    961                return a.sign;
    962            }
    963
    964            if (a.nans == FINITE && b.nans == INFINITE) {
    965                return -b.sign;
    966            }
    967
    968            if (a.nans == INFINITE && b.nans == INFINITE) {
    969                return 0;
    970            }
    971
    972            // Handle special case when a or b is zero, by ignoring the exponents
    973            if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) {
    974                if (a.exp < b.exp) {
    975                    return -a.sign;
    976                }
    977
    978                if (a.exp > b.exp) {
    979                    return a.sign;
    980                }
    981            }
    982
    983            // compare the mantissas
    984            for (int i = a.mant.length - 1; i >= 0; i--) {
    985                if (a.mant[i] > b.mant[i]) {
    986                    return a.sign;
    987                }
    988
    989                if (a.mant[i] < b.mant[i]) {
    990                    return -a.sign;
    991                }
    992            }
    993
    994            return 0;
    995
    996        }
    997
    998        /** Round to nearest integer using the round-half-even method.
    999         *  That is round to nearest integer unless both are equidistant.
    1000         *  In which case round to the even one.
    1001         *  @return rounded value
    1002         */
    1003        public Dfp rint() {
    1004            return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
    1005        }
    1006
    1007        /** Round to an integer using the round floor mode.
    1008         * That is, round toward -Infinity
    1009         *  @return rounded value
    1010         */
    1011        public Dfp floor() {
    1012            return trunc(DfpField.RoundingMode.ROUND_FLOOR);
    1013        }
    1014
    1015        /** Round to an integer using the round ceil mode.
    1016         * That is, round toward +Infinity
    1017         *  @return rounded value
    1018         */
    1019        public Dfp ceil() {
    1020            return trunc(DfpField.RoundingMode.ROUND_CEIL);
    1021        }
    1022
    1023        /** Returns the IEEE remainder.
    1024         * @param d divisor
    1025         * @return this less n &times; d, where n is the integer closest to this/d
    1026         */
    1027        public Dfp remainder(final Dfp d) {
    1028
    1029            final Dfp result = this.subtract(this.divide(d).rint().multiply(d));
    1030
    1031            // IEEE 854-1987 says that if the result is zero, then it carries the sign of this
    1032            if (result.mant[mant.length-1] == 0) {
    1033                result.sign = sign;
    1034            }
    1035
    1036            return result;
    1037
    1038        }
    1039
    1040        /** Does the integer conversions with the specified rounding.
    1041         * @param rmode rounding mode to use
    1042         * @return truncated value
    1043         */
    1044        protected Dfp trunc(final DfpField.RoundingMode rmode) {
    1045            boolean changed = false;
    1046
    1047            if (isNaN()) {
    1048                return newInstance(this);
    1049            }
    1050
    1051            if (nans == INFINITE) {
    1052                return newInstance(this);
    1053            }
    1054
    1055            if (mant[mant.length-1] == 0) {
    1056                // a is zero
    1057                return newInstance(this);
    1058            }
    1059
    1060            /* If the exponent is less than zero then we can certainly
    1061             * return zero */
    1062            if (exp < 0) {
    1063                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
    1064                Dfp result = newInstance(getZero());
    1065                result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
    1066                return result;
    1067            }
    1068
    1069            /* If the exponent is greater than or equal to digits, then it
    1070             * must already be an integer since there is no precision left
    1071             * for any fractional part */
    1072
    1073            if (exp >= mant.length) {
    1074                return newInstance(this);
    1075            }
    1076
    1077            /* General case:  create another dfp, result, that contains the
    1078             * a with the fractional part lopped off.  */
    1079
    1080            Dfp result = newInstance(this);
    1081            for (int i = 0; i < mant.length-result.exp; i++) {
    1082                changed |= result.mant[i] != 0;
    1083                result.mant[i] = 0;
    1084            }
    1085
    1086            if (changed) {
    1087                switch (rmode) {
    1088                    case ROUND_FLOOR:
    1089                        if (result.sign == -1) {
    1090                            // then we must increment the mantissa by one
    1091                            result = result.add(newInstance(-1));
    1092                        }
    1093                        break;
    1094
    1095                    case ROUND_CEIL:
    1096                        if (result.sign == 1) {
    1097                            // then we must increment the mantissa by one
    1098                            result = result.add(getOne());
    1099                        }
    1100                        break;
    1101
    1102                    case ROUND_HALF_EVEN:
    1103                    default:
    1104                        final Dfp half = newInstance("0.5");
    1105                        Dfp a = subtract(result);  // difference between this and result
    1106                        a.sign = 1;            // force positive (take abs)
    1107                        if (a.greaterThan(half)) {
    1108                            a = newInstance(getOne());
    1109                            a.sign = sign;
    1110                            result = result.add(a);
    1111                        }
    1112
    1113                        /** If exactly equal to 1/2 and odd then increment */
    1114                        if (a.equals(half) && result.exp > 0 && (result.mant[mant.length-result.exp]&1) != 0) {
    1115                            a = newInstance(getOne());
    1116                            a.sign = sign;
    1117                            result = result.add(a);
    1118                        }
    1119                        break;
    1120                }
    1121
    1122                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);  // signal inexact
    1123                result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
    1124                return result;
    1125            }
    1126
    1127            return result;
    1128        }
    1129
    1130        /** Convert this to an integer.
    1131         * If greater than 2147483647, it returns 2147483647. If less than -2147483648 it returns -2147483648.
    1132         * @return converted number
    1133         */
    1134        public int intValue() {
    1135            Dfp rounded;
    1136            int result = 0;
    1137
    1138            rounded = rint();
    1139
    1140            if (rounded.greaterThan(newInstance(2147483647))) {
    1141                return 2147483647;
    1142            }
    1143
    1144            if (rounded.lessThan(newInstance(-2147483648))) {
    1145                return -2147483648;
    1146            }
    1147
    1148            for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
    1149                result = result * RADIX + rounded.mant[i];
    1150            }
    1151
    1152            if (rounded.sign == -1) {
    1153                result = -result;
    1154            }
    1155
    1156            return result;
    1157        }
    1158
    1159        /** Get the exponent of the greatest power of 10000 that is
    1160         *  less than or equal to the absolute value of this.  I.E.  if
    1161         *  this is 10<sup>6</sup> then log10K would return 1.
    1162         *  @return integer base 10000 logarithm
    1163         */
    1164        public int log10K() {
    1165            return exp - 1;
    1166        }
    1167
    1168        /** Get the specified  power of 10000.
    1169         * @param e desired power
    1170         * @return 10000<sup>e</sup>
    1171         */
    1172        public Dfp power10K(final int e) {
    1173            Dfp d = newInstance(getOne());
    1174            d.exp = e + 1;
    1175            return d;
    1176        }
    1177
    1178        /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
    1179         *  @return integer base 10 logarithm
    1180         */
    1181        public int log10()  {
    1182            if (mant[mant.length-1] > 1000) {
    1183                return exp * 4 - 1;
    1184            }
    1185            if (mant[mant.length-1] > 100) {
    1186                return exp * 4 - 2;
    1187            }
    1188            if (mant[mant.length-1] > 10) {
    1189                return exp * 4 - 3;
    1190            }
    1191            return exp * 4 - 4;
    1192        }
    1193
    1194        /** Return the specified  power of 10.
    1195         * @param e desired power
    1196         * @return 10<sup>e</sup>
    1197         */
    1198        public Dfp power10(final int e) {
    1199            Dfp d = newInstance(getOne());
    1200
    1201            if (e >= 0) {
    1202                d.exp = e / 4 + 1;
    1203            } else {
    1204                d.exp = (e + 1) / 4;
    1205            }
    1206
    1207            switch ((e % 4 + 4) % 4) {
    1208                case 0:
    1209                    break;
    1210                case 1:
    1211                    d = d.multiply(10);
    1212                    break;
    1213                case 2:
    1214                    d = d.multiply(100);
    1215                    break;
    1216                default:
    1217                    d = d.multiply(1000);
    1218            }
    1219
    1220            return d;
    1221        }
    1222
    1223        /** Negate the mantissa of this by computing the complement.
    1224         *  Leaves the sign bit unchanged, used internally by add.
    1225         *  Denormalized numbers are handled properly here.
    1226         *  @param extra ???
    1227         *  @return ???
    1228         */
    1229        protected int complement(int extra) {
    1230
    1231            extra = RADIX-extra;
    1232            for (int i = 0; i < mant.length; i++) {
    1233                mant[i] = RADIX-mant[i]-1;
    1234            }
    1235
    1236            int rh = extra / RADIX;
    1237            extra = extra - rh * RADIX;
    1238            for (int i = 0; i < mant.length; i++) {
    1239                final int r = mant[i] + rh;
    1240                rh = r / RADIX;
    1241                mant[i] = r - rh * RADIX;
    1242            }
    1243
    1244            return extra;
    1245        }
    1246
    1247        /** Add x to this.
    1248         * @param x number to add
    1249         * @return sum of this and x
    1250         */
    1251        public Dfp add(final Dfp x) {
    1252
    1253            // make sure we don't mix number with different precision
    1254            if (field.getRadixDigits() != x.field.getRadixDigits()) {
    1255                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1256                final Dfp result = newInstance(getZero());
    1257                result.nans = QNAN;
    1258                return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
    1259            }
    1260
    1261            /* handle special cases */
    1262            if (nans != FINITE || x.nans != FINITE) {
    1263                if (isNaN()) {
    1264                    return this;
    1265                }
    1266
    1267                if (x.isNaN()) {
    1268                    return x;
    1269                }
    1270
    1271                if (nans == INFINITE && x.nans == FINITE) {
    1272                    return this;
    1273                }
    1274
    1275                if (x.nans == INFINITE && nans == FINITE) {
    1276                    return x;
    1277                }
    1278
    1279                if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
    1280                    return x;
    1281                }
    1282
    1283                if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
    1284                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1285                    Dfp result = newInstance(getZero());
    1286                    result.nans = QNAN;
    1287                    result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
    1288                    return result;
    1289                }
    1290            }
    1291
    1292            /* copy this and the arg */
    1293            Dfp a = newInstance(this);
    1294            Dfp b = newInstance(x);
    1295
    1296            /* initialize the result object */
    1297            Dfp result = newInstance(getZero());
    1298
    1299            /* Make all numbers positive, but remember their sign */
    1300            final byte asign = a.sign;
    1301            final byte bsign = b.sign;
    1302
    1303            a.sign = 1;
    1304            b.sign = 1;
    1305
    1306            /* The result will be signed like the arg with greatest magnitude */
    1307            byte rsign = bsign;
    1308            if (compare(a, b) > 0) {
    1309                rsign = asign;
    1310            }
    1311
    1312            /* Handle special case when a or b is zero, by setting the exponent
    1313           of the zero number equal to the other one.  This avoids an alignment
    1314           which would cause catastropic loss of precision */
    1315            if (b.mant[mant.length-1] == 0) {
    1316                b.exp = a.exp;
    1317            }
    1318
    1319            if (a.mant[mant.length-1] == 0) {
    1320                a.exp = b.exp;
    1321            }
    1322
    1323            /* align number with the smaller exponent */
    1324            int aextradigit = 0;
    1325            int bextradigit = 0;
    1326            if (a.exp < b.exp) {
    1327                aextradigit = a.align(b.exp);
    1328            } else {
    1329                bextradigit = b.align(a.exp);
    1330            }
    1331
    1332            /* complement the smaller of the two if the signs are different */
    1333            if (asign != bsign) {
    1334                if (asign == rsign) {
    1335                    bextradigit = b.complement(bextradigit);
    1336                } else {
    1337                    aextradigit = a.complement(aextradigit);
    1338                }
    1339            }
    1340
    1341            /* add the mantissas */
    1342            int rh = 0; /* acts as a carry */
    1343            for (int i = 0; i < mant.length; i++) {
    1344                final int r = a.mant[i]+b.mant[i]+rh;
    1345                rh = r / RADIX;
    1346                result.mant[i] = r - rh * RADIX;
    1347            }
    1348            result.exp = a.exp;
    1349            result.sign = rsign;
    1350
    1351            /* handle overflow -- note, when asign!=bsign an overflow is
    1352             * normal and should be ignored.  */
    1353
    1354            if (rh != 0 && (asign == bsign)) {
    1355                final int lostdigit = result.mant[0];
    1356                result.shiftRight();
    1357                result.mant[mant.length-1] = rh;
    1358                final int excp = result.round(lostdigit);
    1359                if (excp != 0) {
    1360                    result = dotrap(excp, ADD_TRAP, x, result);
    1361                }
    1362            }
    1363
    1364            /* normalize the result */
    1365            for (int i = 0; i < mant.length; i++) {
    1366                if (result.mant[mant.length-1] != 0) {
    1367                    break;
    1368                }
    1369                result.shiftLeft();
    1370                if (i == 0) {
    1371                    result.mant[0] = aextradigit+bextradigit;
    1372                    aextradigit = 0;
    1373                    bextradigit = 0;
    1374                }
    1375            }
    1376
    1377            /* result is zero if after normalization the most sig. digit is zero */
    1378            if (result.mant[mant.length-1] == 0) {
    1379                result.exp = 0;
    1380
    1381                if (asign != bsign) {
    1382                    // Unless adding 2 negative zeros, sign is positive
    1383                    result.sign = 1;  // Per IEEE 854-1987 Section 6.3
    1384                }
    1385            }
    1386
    1387            /* Call round to test for over/under flows */
    1388            final int excp = result.round(aextradigit + bextradigit);
    1389            if (excp != 0) {
    1390                result = dotrap(excp, ADD_TRAP, x, result);
    1391            }
    1392
    1393            return result;
    1394        }
    1395
    1396        /** Returns a number that is this number with the sign bit reversed.
    1397         * @return the opposite of this
    1398         */
    1399        public Dfp negate() {
    1400            Dfp result = newInstance(this);
    1401            result.sign = (byte) - result.sign;
    1402            return result;
    1403        }
    1404
    1405        /** Subtract x from this.
    1406         * @param x number to subtract
    1407         * @return difference of this and a
    1408         */
    1409        public Dfp subtract(final Dfp x) {
    1410            return add(x.negate());
    1411        }
    1412
    1413        /** Round this given the next digit n using the current rounding mode.
    1414         * @param n ???
    1415         * @return the IEEE flag if an exception occurred
    1416         */
    1417        protected int round(int n) {
    1418            boolean inc = false;
    1419            switch (field.getRoundingMode()) {
    1420                case ROUND_DOWN:
    1421                    inc = false;
    1422                    break;
    1423
    1424                case ROUND_UP:
    1425                    inc = n != 0;       // round up if n!=0
    1426                    break;
    1427
    1428                case ROUND_HALF_UP:
    1429                    inc = n >= 5000;  // round half up
    1430                    break;
    1431
    1432                case ROUND_HALF_DOWN:
    1433                    inc = n > 5000;  // round half down
    1434                    break;
    1435
    1436                case ROUND_HALF_EVEN:
    1437                    inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1);  // round half-even
    1438                    break;
    1439
    1440                case ROUND_HALF_ODD:
    1441                    inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0);  // round half-odd
    1442                    break;
    1443
    1444                case ROUND_CEIL:
    1445                    inc = sign == 1 && n != 0;  // round ceil
    1446                    break;
    1447
    1448                case ROUND_FLOOR:
    1449                default:
    1450                    inc = sign == -1 && n != 0;  // round floor
    1451                    break;
    1452            }
    1453
    1454            if (inc) {
    1455                // increment if necessary
    1456                int rh = 1;
    1457                for (int i = 0; i < mant.length; i++) {
    1458                    final int r = mant[i] + rh;
    1459                    rh = r / RADIX;
    1460                    mant[i] = r - rh * RADIX;
    1461                }
    1462
    1463                if (rh != 0) {
    1464                    shiftRight();
    1465                    mant[mant.length-1] = rh;
    1466                }
    1467            }
    1468
    1469            // check for exceptional cases and raise signals if necessary
    1470            if (exp < MIN_EXP) {
    1471                // Gradual Underflow
    1472                field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
    1473                return DfpField.FLAG_UNDERFLOW;
    1474            }
    1475
    1476            if (exp > MAX_EXP) {
    1477                // Overflow
    1478                field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
    1479                return DfpField.FLAG_OVERFLOW;
    1480            }
    1481
    1482            if (n != 0) {
    1483                // Inexact
    1484                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
    1485                return DfpField.FLAG_INEXACT;
    1486            }
    1487
    1488            return 0;
    1489
    1490        }
    1491
    1492        /** Multiply this by x.
    1493         * @param x multiplicand
    1494         * @return product of this and x
    1495         */
    1496        public Dfp multiply(final Dfp x) {
    1497
    1498            // make sure we don't mix number with different precision
    1499            if (field.getRadixDigits() != x.field.getRadixDigits()) {
    1500                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1501                final Dfp result = newInstance(getZero());
    1502                result.nans = QNAN;
    1503                return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
    1504            }
    1505
    1506            Dfp result = newInstance(getZero());
    1507
    1508            /* handle special cases */
    1509            if (nans != FINITE || x.nans != FINITE) {
    1510                if (isNaN()) {
    1511                    return this;
    1512                }
    1513
    1514                if (x.isNaN()) {
    1515                    return x;
    1516                }
    1517
    1518                if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] != 0) {
    1519                    result = newInstance(this);
    1520                    result.sign = (byte) (sign * x.sign);
    1521                    return result;
    1522                }
    1523
    1524                if (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] != 0) {
    1525                    result = newInstance(x);
    1526                    result.sign = (byte) (sign * x.sign);
    1527                    return result;
    1528                }
    1529
    1530                if (x.nans == INFINITE && nans == INFINITE) {
    1531                    result = newInstance(this);
    1532                    result.sign = (byte) (sign * x.sign);
    1533                    return result;
    1534                }
    1535
    1536                if ( (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] == 0) ||
    1537                        (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] == 0) ) {
    1538                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1539                    result = newInstance(getZero());
    1540                    result.nans = QNAN;
    1541                    result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
    1542                    return result;
    1543                }
    1544            }
    1545
    1546            int[] product = new int[mant.length*2];  // Big enough to hold even the largest result
    1547
    1548            for (int i = 0; i < mant.length; i++) {
    1549                int rh = 0;  // acts as a carry
    1550                for (int j=0; j<mant.length; j++) {
    1551                    int r = mant[i] * x.mant[j];    // multiply the 2 digits
    1552                    r = r + product[i+j] + rh;  // add to the product digit with carry in
    1553
    1554                    rh = r / RADIX;
    1555                    product[i+j] = r - rh * RADIX;
    1556                }
    1557                product[i+mant.length] = rh;
    1558            }
    1559
    1560            // Find the most sig digit
    1561            int md = mant.length * 2 - 1;  // default, in case result is zero
    1562            for (int i = mant.length * 2 - 1; i >= 0; i--) {
    1563                if (product[i] != 0) {
    1564                    md = i;
    1565                    break;
    1566                }
    1567            }
    1568
    1569            // Copy the digits into the result
    1570            for (int i = 0; i < mant.length; i++) {
    1571                result.mant[mant.length - i - 1] = product[md - i];
    1572            }
    1573
    1574            // Fixup the exponent.
    1575            result.exp = exp + x.exp + md - 2 * mant.length + 1;
    1576            result.sign = (byte)((sign == x.sign)?1:-1);
    1577
    1578            if (result.mant[mant.length-1] == 0) {
    1579                // if result is zero, set exp to zero
    1580                result.exp = 0;
    1581            }
    1582
    1583            final int excp;
    1584            if (md > (mant.length-1)) {
    1585                excp = result.round(product[md-mant.length]);
    1586            } else {
    1587                excp = result.round(0); // has no effect except to check status
    1588            }
    1589
    1590            if (excp != 0) {
    1591                result = dotrap(excp, MULTIPLY_TRAP, x, result);
    1592            }
    1593
    1594            return result;
    1595
    1596        }
    1597
    1598        /** Multiply this by a single digit x.
    1599         * @param x multiplicand
    1600         * @return product of this and x
    1601         */
    1602        public Dfp multiply(final int x) {
    1603            if (x >= 0 && x < RADIX) {
    1604                return multiplyFast(x);
    1605            } else {
    1606                return multiply(newInstance(x));
    1607            }
    1608        }
    1609
    1610        /** Multiply this by a single digit 0&lt;=x&lt;radix.
    1611         * There are speed advantages in this special case.
    1612         * @param x multiplicand
    1613         * @return product of this and x
    1614         */
    1615        private Dfp multiplyFast(final int x) {
    1616            Dfp result = newInstance(this);
    1617
    1618            /* handle special cases */
    1619            if (nans != FINITE) {
    1620                if (isNaN()) {
    1621                    return this;
    1622                }
    1623
    1624                if (nans == INFINITE && x != 0) {
    1625                    result = newInstance(this);
    1626                    return result;
    1627                }
    1628
    1629                if (nans == INFINITE && x == 0) {
    1630                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1631                    result = newInstance(getZero());
    1632                    result.nans = QNAN;
    1633                    result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
    1634                    return result;
    1635                }
    1636            }
    1637
    1638            /* range check x */
    1639            if (x < 0 || x >= RADIX) {
    1640                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1641                result = newInstance(getZero());
    1642                result.nans = QNAN;
    1643                result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
    1644                return result;
    1645            }
    1646
    1647            int rh = 0;
    1648            for (int i = 0; i < mant.length; i++) {
    1649                final int r = mant[i] * x + rh;
    1650                rh = r / RADIX;
    1651                result.mant[i] = r - rh * RADIX;
    1652            }
    1653
    1654            int lostdigit = 0;
    1655            if (rh != 0) {
    1656                lostdigit = result.mant[0];
    1657                result.shiftRight();
    1658                result.mant[mant.length-1] = rh;
    1659            }
    1660
    1661            if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
    1662                result.exp = 0;
    1663            }
    1664
    1665            final int excp = result.round(lostdigit);
    1666            if (excp != 0) {
    1667                result = dotrap(excp, MULTIPLY_TRAP, result, result);
    1668            }
    1669
    1670            return result;
    1671        }
    1672
    1673        /** Divide this by divisor.
    1674         * @param divisor divisor
    1675         * @return quotient of this by divisor
    1676         */
    1677        public Dfp divide(Dfp divisor) {
    1678            int dividend[]; // current status of the dividend
    1679            int quotient[]; // quotient
    1680            int remainder[];// remainder
    1681            int qd;         // current quotient digit we're working with
    1682            int nsqd;       // number of significant quotient digits we have
    1683            int trial=0;    // trial quotient digit
    1684            int minadj;     // minimum adjustment
    1685            boolean trialgood; // Flag to indicate a good trail digit
    1686            int md=0;       // most sig digit in result
    1687            int excp;       // exceptions
    1688
    1689            // make sure we don't mix number with different precision
    1690            if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
    1691                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1692                final Dfp result = newInstance(getZero());
    1693                result.nans = QNAN;
    1694                return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
    1695            }
    1696
    1697            Dfp result = newInstance(getZero());
    1698
    1699            /* handle special cases */
    1700            if (nans != FINITE || divisor.nans != FINITE) {
    1701                if (isNaN()) {
    1702                    return this;
    1703                }
    1704
    1705                if (divisor.isNaN()) {
    1706                    return divisor;
    1707                }
    1708
    1709                if (nans == INFINITE && divisor.nans == FINITE) {
    1710                    result = newInstance(this);
    1711                    result.sign = (byte) (sign * divisor.sign);
    1712                    return result;
    1713                }
    1714
    1715                if (divisor.nans == INFINITE && nans == FINITE) {
    1716                    result = newInstance(getZero());
    1717                    result.sign = (byte) (sign * divisor.sign);
    1718                    return result;
    1719                }
    1720
    1721                if (divisor.nans == INFINITE && nans == INFINITE) {
    1722                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1723                    result = newInstance(getZero());
    1724                    result.nans = QNAN;
    1725                    result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
    1726                    return result;
    1727                }
    1728            }
    1729
    1730            /* Test for divide by zero */
    1731            if (divisor.mant[mant.length-1] == 0) {
    1732                field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
    1733                result = newInstance(getZero());
    1734                result.sign = (byte) (sign * divisor.sign);
    1735                result.nans = INFINITE;
    1736                result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
    1737                return result;
    1738            }
    1739
    1740            dividend = new int[mant.length+1];  // one extra digit needed
    1741            quotient = new int[mant.length+2];  // two extra digits needed 1 for overflow, 1 for rounding
    1742            remainder = new int[mant.length+1]; // one extra digit needed
    1743
    1744            /* Initialize our most significant digits to zero */
    1745
    1746            dividend[mant.length] = 0;
    1747            quotient[mant.length] = 0;
    1748            quotient[mant.length+1] = 0;
    1749            remainder[mant.length] = 0;
    1750
    1751            /* copy our mantissa into the dividend, initialize the
    1752           quotient while we are at it */
    1753
    1754            for (int i = 0; i < mant.length; i++) {
    1755                dividend[i] = mant[i];
    1756                quotient[i] = 0;
    1757                remainder[i] = 0;
    1758            }
    1759
    1760            /* outer loop.  Once per quotient digit */
    1761            nsqd = 0;
    1762            for (qd = mant.length+1; qd >= 0; qd--) {
    1763                /* Determine outer limits of our quotient digit */
    1764
    1765                // r =  most sig 2 digits of dividend
    1766                final int divMsb = dividend[mant.length]*RADIX+dividend[mant.length-1];
    1767                int min = divMsb       / (divisor.mant[mant.length-1]+1);
    1768                int max = (divMsb + 1) / divisor.mant[mant.length-1];
    1769
    1770                trialgood = false;
    1771                while (!trialgood) {
    1772                    // try the mean
    1773                    trial = (min+max)/2;
    1774
    1775                    /* Multiply by divisor and store as remainder */
    1776                    int rh = 0;
    1777                    for (int i = 0; i < mant.length + 1; i++) {
    1778                        int dm = (i<mant.length)?divisor.mant[i]:0;
    1779                        final int r = (dm * trial) + rh;
    1780                        rh = r / RADIX;
    1781                        remainder[i] = r - rh * RADIX;
    1782                    }
    1783
    1784                    /* subtract the remainder from the dividend */
    1785                    rh = 1;  // carry in to aid the subtraction
    1786                    for (int i = 0; i < mant.length + 1; i++) {
    1787                        final int r = ((RADIX-1) - remainder[i]) + dividend[i] + rh;
    1788                        rh = r / RADIX;
    1789                        remainder[i] = r - rh * RADIX;
    1790                    }
    1791
    1792                    /* Lets analyze what we have here */
    1793                    if (rh == 0) {
    1794                        // trial is too big -- negative remainder
    1795                        max = trial-1;
    1796                        continue;
    1797                    }
    1798
    1799                    /* find out how far off the remainder is telling us we are */
    1800                    minadj = (remainder[mant.length] * RADIX)+remainder[mant.length-1];
    1801                    minadj = minadj / (divisor.mant[mant.length-1]+1);
    1802
    1803                    if (minadj >= 2) {
    1804                        min = trial+minadj;  // update the minimum
    1805                        continue;
    1806                    }
    1807
    1808                    /* May have a good one here, check more thoroughly.  Basically
    1809               its a good one if it is less than the divisor */
    1810                    trialgood = false;  // assume false
    1811                    for (int i = mant.length - 1; i >= 0; i--) {
    1812                        if (divisor.mant[i] > remainder[i]) {
    1813                            trialgood = true;
    1814                        }
    1815                        if (divisor.mant[i] < remainder[i]) {
    1816                            break;
    1817                        }
    1818                    }
    1819
    1820                    if (remainder[mant.length] != 0) {
    1821                        trialgood = false;
    1822                    }
    1823
    1824                    if (trialgood == false) {
    1825                        min = trial+1;
    1826                    }
    1827                }
    1828
    1829                /* Great we have a digit! */
    1830                quotient[qd] = trial;
    1831                if (trial != 0 || nsqd != 0) {
    1832                    nsqd++;
    1833                }
    1834
    1835                if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
    1836                    // We have enough for this mode
    1837                    break;
    1838                }
    1839
    1840                if (nsqd > mant.length) {
    1841                    // We have enough digits
    1842                    break;
    1843                }
    1844
    1845                /* move the remainder into the dividend while left shifting */
    1846                dividend[0] = 0;
    1847                for (int i = 0; i < mant.length; i++) {
    1848                    dividend[i + 1] = remainder[i];
    1849                }
    1850            }
    1851
    1852            /* Find the most sig digit */
    1853            md = mant.length;  // default
    1854            for (int i = mant.length + 1; i >= 0; i--) {
    1855                if (quotient[i] != 0) {
    1856                    md = i;
    1857                    break;
    1858                }
    1859            }
    1860
    1861            /* Copy the digits into the result */
    1862            for (int i=0; i<mant.length; i++) {
    1863                result.mant[mant.length-i-1] = quotient[md-i];
    1864            }
    1865
    1866            /* Fixup the exponent. */
    1867            result.exp = exp - divisor.exp + md - mant.length;
    1868            result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
    1869
    1870            if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
    1871                result.exp = 0;
    1872            }
    1873
    1874            if (md > (mant.length-1)) {
    1875                excp = result.round(quotient[md-mant.length]);
    1876            } else {
    1877                excp = result.round(0);
    1878            }
    1879
    1880            if (excp != 0) {
    1881                result = dotrap(excp, DIVIDE_TRAP, divisor, result);
    1882            }
    1883
    1884            return result;
    1885        }
    1886
    1887        /** Divide by a single digit less than radix.
    1888         *  Special case, so there are speed advantages. 0 &lt;= divisor &lt; radix
    1889         * @param divisor divisor
    1890         * @return quotient of this by divisor
    1891         */
    1892        public Dfp divide(int divisor) {
    1893
    1894            // Handle special cases
    1895            if (nans != FINITE) {
    1896                if (isNaN()) {
    1897                    return this;
    1898                }
    1899
    1900                if (nans == INFINITE) {
    1901                    return newInstance(this);
    1902                }
    1903            }
    1904
    1905            // Test for divide by zero
    1906            if (divisor == 0) {
    1907                field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
    1908                Dfp result = newInstance(getZero());
    1909                result.sign = sign;
    1910                result.nans = INFINITE;
    1911                result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
    1912                return result;
    1913            }
    1914
    1915            // range check divisor
    1916            if (divisor < 0 || divisor >= RADIX) {
    1917                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1918                Dfp result = newInstance(getZero());
    1919                result.nans = QNAN;
    1920                result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
    1921                return result;
    1922            }
    1923
    1924            Dfp result = newInstance(this);
    1925
    1926            int rl = 0;
    1927            for (int i = mant.length-1; i >= 0; i--) {
    1928                final int r = rl*RADIX + result.mant[i];
    1929                final int rh = r / divisor;
    1930                rl = r - rh * divisor;
    1931                result.mant[i] = rh;
    1932            }
    1933
    1934            if (result.mant[mant.length-1] == 0) {
    1935                // normalize
    1936                result.shiftLeft();
    1937                final int r = rl * RADIX;        // compute the next digit and put it in
    1938                final int rh = r / divisor;
    1939                rl = r - rh * divisor;
    1940                result.mant[0] = rh;
    1941            }
    1942
    1943            final int excp = result.round(rl * RADIX / divisor);  // do the rounding
    1944            if (excp != 0) {
    1945                result = dotrap(excp, DIVIDE_TRAP, result, result);
    1946            }
    1947
    1948            return result;
    1949
    1950        }
    1951
    1952        /** {@inheritDoc} */
    1953        public Dfp reciprocal() {
    1954            return field.getOne().divide(this);
    1955        }
    1956
    1957        /** Compute the square root.
    1958         * @return square root of the instance
    1959         */
    1960        public Dfp sqrt() {
    1961
    1962            // check for unusual cases
    1963            if (nans == FINITE && mant[mant.length-1] == 0) {
    1964                // if zero
    1965                return newInstance(this);
    1966            }
    1967
    1968            if (nans != FINITE) {
    1969                if (nans == INFINITE && sign == 1) {
    1970                    // if positive infinity
    1971                    return newInstance(this);
    1972                }
    1973
    1974                if (nans == QNAN) {
    1975                    return newInstance(this);
    1976                }
    1977
    1978                if (nans == SNAN) {
    1979                    Dfp result;
    1980
    1981                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1982                    result = newInstance(this);
    1983                    result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
    1984                    return result;
    1985                }
    1986            }
    1987
    1988            if (sign == -1) {
    1989                // if negative
    1990                Dfp result;
    1991
    1992                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    1993                result = newInstance(this);
    1994                result.nans = QNAN;
    1995                result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
    1996                return result;
    1997            }
    1998
    1999            Dfp x = newInstance(this);
    2000
    2001            /* Lets make a reasonable guess as to the size of the square root */
    2002            if (x.exp < -1 || x.exp > 1) {
    2003                x.exp = this.exp / 2;
    2004            }
    2005
    2006            /* Coarsely estimate the mantissa */
    2007            switch (x.mant[mant.length-1] / 2000) {
    2008                case 0:
    2009                    x.mant[mant.length-1] = x.mant[mant.length-1]/2+1;
    2010                    break;
    2011                case 2:
    2012                    x.mant[mant.length-1] = 1500;
    2013                    break;
    2014                case 3:
    2015                    x.mant[mant.length-1] = 2200;
    2016                    break;
    2017                default:
    2018                    x.mant[mant.length-1] = 3000;
    2019            }
    2020
    2021            Dfp dx = newInstance(x);
    2022
    2023            /* Now that we have the first pass estimate, compute the rest
    2024           by the formula dx = (y - x*x) / (2x); */
    2025
    2026            Dfp px  = getZero();
    2027            Dfp ppx = getZero();
    2028            while (x.unequal(px)) {
    2029                dx = newInstance(x);
    2030                dx.sign = -1;
    2031                dx = dx.add(this.divide(x));
    2032                dx = dx.divide(2);
    2033                ppx = px;
    2034                px = x;
    2035                x = x.add(dx);
    2036
    2037                if (x.equals(ppx)) {
    2038                    // alternating between two values
    2039                    break;
    2040                }
    2041
    2042                // if dx is zero, break.  Note testing the most sig digit
    2043                // is a sufficient test since dx is normalized
    2044                if (dx.mant[mant.length-1] == 0) {
    2045                    break;
    2046                }
    2047            }
    2048
    2049            return x;
    2050
    2051        }
    2052
    2053        /** Get a string representation of the instance.
    2054         * @return string representation of the instance
    2055         */
    2056        @Override
    2057        public String toString() {
    2058            if (nans != FINITE) {
    2059                // if non-finite exceptional cases
    2060                if (nans == INFINITE) {
    2061                    return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
    2062                } else {
    2063                    return NAN_STRING;
    2064                }
    2065            }
    2066
    2067            if (exp > mant.length || exp < -1) {
    2068                return dfp2sci();
    2069            }
    2070
    2071            return dfp2string();
    2072
    2073        }
    2074
    2075        /** Convert an instance to a string using scientific notation.
    2076         * @return string representation of the instance in scientific notation
    2077         */
    2078        protected String dfp2sci() {
    2079            char rawdigits[]    = new char[mant.length * 4];
    2080            char outputbuffer[] = new char[mant.length * 4 + 20];
    2081            int p;
    2082            int q;
    2083            int e;
    2084            int ae;
    2085            int shf;
    2086
    2087            // Get all the digits
    2088            p = 0;
    2089            for (int i = mant.length - 1; i >= 0; i--) {
    2090                rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
    2091                rawdigits[p++] = (char) (((mant[i] / 100) %10) + '0');
    2092                rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
    2093                rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
    2094            }
    2095
    2096            // Find the first non-zero one
    2097            for (p = 0; p < rawdigits.length; p++) {
    2098                if (rawdigits[p] != '0') {
    2099                    break;
    2100                }
    2101            }
    2102            shf = p;
    2103
    2104            // Now do the conversion
    2105            q = 0;
    2106            if (sign == -1) {
    2107                outputbuffer[q++] = '-';
    2108            }
    2109
    2110            if (p != rawdigits.length) {
    2111                // there are non zero digits...
    2112                outputbuffer[q++] = rawdigits[p++];
    2113                outputbuffer[q++] = '.';
    2114
    2115                while (p<rawdigits.length) {
    2116                    outputbuffer[q++] = rawdigits[p++];
    2117                }
    2118            } else {
    2119                outputbuffer[q++] = '0';
    2120                outputbuffer[q++] = '.';
    2121                outputbuffer[q++] = '0';
    2122                outputbuffer[q++] = 'e';
    2123                outputbuffer[q++] = '0';
    2124                return new String(outputbuffer, 0, 5);
    2125            }
    2126
    2127            outputbuffer[q++] = 'e';
    2128
    2129            // Find the msd of the exponent
    2130
    2131            e = exp * 4 - shf - 1;
    2132            ae = e;
    2133            if (e < 0) {
    2134                ae = -e;
    2135            }
    2136
    2137            // Find the largest p such that p < e
    2138            for (p = 1000000000; p > ae; p /= 10) {
    2139                // nothing to do
    2140            }
    2141
    2142            if (e < 0) {
    2143                outputbuffer[q++] = '-';
    2144            }
    2145
    2146            while (p > 0) {
    2147                outputbuffer[q++] = (char)(ae / p + '0');
    2148                ae = ae % p;
    2149                p = p / 10;
    2150            }
    2151
    2152            return new String(outputbuffer, 0, q);
    2153
    2154        }
    2155
    2156        /** Convert an instance to a string using normal notation.
    2157         * @return string representation of the instance in normal notation
    2158         */
    2159        protected String dfp2string() {
    2160            char buffer[] = new char[mant.length*4 + 20];
    2161            int p = 1;
    2162            int q;
    2163            int e = exp;
    2164            boolean pointInserted = false;
    2165
    2166            buffer[0] = ' ';
    2167
    2168            if (e <= 0) {
    2169                buffer[p++] = '0';
    2170                buffer[p++] = '.';
    2171                pointInserted = true;
    2172            }
    2173
    2174            while (e < 0) {
    2175                buffer[p++] = '0';
    2176                buffer[p++] = '0';
    2177                buffer[p++] = '0';
    2178                buffer[p++] = '0';
    2179                e++;
    2180            }
    2181
    2182            for (int i = mant.length - 1; i >= 0; i--) {
    2183                buffer[p++] = (char) ((mant[i] / 1000) + '0');
    2184                buffer[p++] = (char) (((mant[i] / 100) % 10) + '0');
    2185                buffer[p++] = (char) (((mant[i] / 10) % 10) + '0');
    2186                buffer[p++] = (char) (((mant[i]) % 10) + '0');
    2187                if (--e == 0) {
    2188                    buffer[p++] = '.';
    2189                    pointInserted = true;
    2190                }
    2191            }
    2192
    2193            while (e > 0) {
    2194                buffer[p++] = '0';
    2195                buffer[p++] = '0';
    2196                buffer[p++] = '0';
    2197                buffer[p++] = '0';
    2198                e--;
    2199            }
    2200
    2201            if (!pointInserted) {
    2202                // Ensure we have a radix point!
    2203                buffer[p++] = '.';
    2204            }
    2205
    2206            // Suppress leading zeros
    2207            q = 1;
    2208            while (buffer[q] == '0') {
    2209                q++;
    2210            }
    2211            if (buffer[q] == '.') {
    2212                q--;
    2213            }
    2214
    2215            // Suppress trailing zeros
    2216            while (buffer[p-1] == '0') {
    2217                p--;
    2218            }
    2219
    2220            // Insert sign
    2221            if (sign < 0) {
    2222                buffer[--q] = '-';
    2223            }
    2224
    2225            return new String(buffer, q, p - q);
    2226
    2227        }
    2228
    2229        /** Raises a trap.  This does not set the corresponding flag however.
    2230         *  @param type the trap type
    2231         *  @param what - name of routine trap occurred in
    2232         *  @param oper - input operator to function
    2233         *  @param result - the result computed prior to the trap
    2234         *  @return The suggested return value from the trap handler
    2235         */
    2236        public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
    2237            Dfp def = result;
    2238
    2239            switch (type) {
    2240                case DfpField.FLAG_INVALID:
    2241                    def = newInstance(getZero());
    2242                    def.sign = result.sign;
    2243                    def.nans = QNAN;
    2244                    break;
    2245
    2246                case DfpField.FLAG_DIV_ZERO:
    2247                    if (nans == FINITE && mant[mant.length-1] != 0) {
    2248                        // normal case, we are finite, non-zero
    2249                        def = newInstance(getZero());
    2250                        def.sign = (byte)(sign*oper.sign);
    2251                        def.nans = INFINITE;
    2252                    }
    2253
    2254                    if (nans == FINITE && mant[mant.length-1] == 0) {
    2255                        //  0/0
    2256                        def = newInstance(getZero());
    2257                        def.nans = QNAN;
    2258                    }
    2259
    2260                    if (nans == INFINITE || nans == QNAN) {
    2261                        def = newInstance(getZero());
    2262                        def.nans = QNAN;
    2263                    }
    2264
    2265                    if (nans == INFINITE || nans == SNAN) {
    2266                        def = newInstance(getZero());
    2267                        def.nans = QNAN;
    2268                    }
    2269                    break;
    2270
    2271                case DfpField.FLAG_UNDERFLOW:
    2272                    if ( (result.exp+mant.length) < MIN_EXP) {
    2273                        def = newInstance(getZero());
    2274                        def.sign = result.sign;
    2275                    } else {
    2276                        def = newInstance(result);  // gradual underflow
    2277                    }
    2278                    result.exp = result.exp + ERR_SCALE;
    2279                    break;
    2280
    2281                case DfpField.FLAG_OVERFLOW:
    2282                    result.exp = result.exp - ERR_SCALE;
    2283                    def = newInstance(getZero());
    2284                    def.sign = result.sign;
    2285                    def.nans = INFINITE;
    2286                    break;
    2287
    2288                default: def = result; break;
    2289            }
    2290
    2291            return trap(type, what, oper, def, result);
    2292
    2293        }
    2294
    2295        /** Trap handler.  Subclasses may override this to provide trap
    2296         *  functionality per IEEE 854-1987.
    2297         *
    2298         *  @param type  The exception type - e.g. FLAG_OVERFLOW
    2299         *  @param what  The name of the routine we were in e.g. divide()
    2300         *  @param oper  An operand to this function if any
    2301         *  @param def   The default return value if trap not enabled
    2302         *  @param result    The result that is specified to be delivered per
    2303         *                   IEEE 854, if any
    2304         *  @return the value that should be return by the operation triggering the trap
    2305         */
    2306        protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
    2307            return def;
    2308        }
    2309
    2310        /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN.
    2311         * @return type of the number
    2312         */
    2313        public int classify() {
    2314            return nans;
    2315        }
    2316
    2317        /** Creates an instance that is the same as x except that it has the sign of y.
    2318         * abs(x) = dfp.copysign(x, dfp.one)
    2319         * @param x number to get the value from
    2320         * @param y number to get the sign from
    2321         * @return a number with the value of x and the sign of y
    2322         */
    2323        public static Dfp copysign(final Dfp x, final Dfp y) {
    2324            Dfp result = x.newInstance(x);
    2325            result.sign = y.sign;
    2326            return result;
    2327        }
    2328
    2329        /** Returns the next number greater than this one in the direction of x.
    2330         * If this==x then simply returns this.
    2331         * @param x direction where to look at
    2332         * @return closest number next to instance in the direction of x
    2333         */
    2334        public Dfp nextAfter(final Dfp x) {
    2335
    2336            // make sure we don't mix number with different precision
    2337            if (field.getRadixDigits() != x.field.getRadixDigits()) {
    2338                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
    2339                final Dfp result = newInstance(getZero());
    2340                result.nans = QNAN;
    2341                return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
    2342            }
    2343
    2344            // if this is greater than x
    2345            boolean up = false;
    2346            if (this.lessThan(x)) {
    2347                up = true;
    2348            }
    2349
    2350            if (compare(this, x) == 0) {
    2351                return newInstance(x);
    2352            }
    2353
    2354            if (lessThan(getZero())) {
    2355                up = !up;
    2356            }
    2357
    2358            final Dfp inc;
    2359            Dfp result;
    2360            if (up) {
    2361                inc = newInstance(getOne());
    2362                inc.exp = this.exp-mant.length+1;
    2363                inc.sign = this.sign;
    2364
    2365                if (this.equals(getZero())) {
    2366                    inc.exp = MIN_EXP-mant.length;
    2367                }
    2368
    2369                result = add(inc);
    2370            } else {
    2371                inc = newInstance(getOne());
    2372                inc.exp = this.exp;
    2373                inc.sign = this.sign;
    2374
    2375                if (this.equals(inc)) {
    2376                    inc.exp = this.exp-mant.length;
    2377                } else {
    2378                    inc.exp = this.exp-mant.length+1;
    2379                }
    2380
    2381                if (this.equals(getZero())) {
    2382                    inc.exp = MIN_EXP-mant.length;
    2383                }
    2384
    2385                result = this.subtract(inc);
    2386            }
    2387
    2388            if (result.classify() == INFINITE && this.classify() != INFINITE) {
    2389                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
    2390                result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
    2391            }
    2392
    2393            if (result.equals(getZero()) && this.equals(getZero()) == false) {
    2394                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
    2395                result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
    2396            }
    2397
    2398            return result;
    2399
    2400        }
    2401
    2402        /** Convert the instance into a double.
    2403         * @return a double approximating the instance
    2404         * @see #toSplitDouble()
    2405         */
    2406        public double toDouble() {
    2407
    2408            if (isInfinite()) {
    2409                if (lessThan(getZero())) {
    2410                    return Double.NEGATIVE_INFINITY;
    2411                } else {
    2412                    return Double.POSITIVE_INFINITY;
    2413                }
    2414            }
    2415
    2416            if (isNaN()) {
    2417                return Double.NaN;
    2418            }
    2419
    2420            Dfp y = this;
    2421            boolean negate = false;
    2422            int cmp0 = compare(this, getZero());
    2423            if (cmp0 == 0) {
    2424                return sign < 0 ? -0.0 : +0.0;
    2425            } else if (cmp0 < 0) {
    2426                y = negate();
    2427                negate = true;
    2428            }
    2429
    2430            /* Find the exponent, first estimate by integer log10, then adjust.
    2431             Should be faster than doing a natural logarithm.  */
    2432            int exponent = (int)(y.log10() * 3.32);
    2433            if (exponent < 0) {
    2434                exponent--;
    2435            }
    2436
    2437            Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
    2438            while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
    2439                tempDfp = tempDfp.multiply(2);
    2440                exponent++;
    2441            }
    2442            exponent--;
    2443
    2444            /* We have the exponent, now work on the mantissa */
    2445
    2446            y = y.divide(DfpMath.pow(getTwo(), exponent));
    2447            if (exponent > -1023) {
    2448                y = y.subtract(getOne());
    2449            }
    2450
    2451            if (exponent < -1074) {
    2452                return 0;
    2453            }
    2454
    2455            if (exponent > 1023) {
    2456                return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
    2457            }
    2458
    2459
    2460            y = y.multiply(newInstance(4503599627370496l)).rint();
    2461            String str = y.toString();
    2462            str = str.substring(0, str.length()-1);
    2463            long mantissa = Long.parseLong(str);
    2464
    2465            if (mantissa == 4503599627370496L) {
    2466                // Handle special case where we round up to next power of two
    2467                mantissa = 0;
    2468                exponent++;
    2469            }
    2470
    2471            /* Its going to be subnormal, so make adjustments */
    2472            if (exponent <= -1023) {
    2473                exponent--;
    2474            }
    2475
    2476            while (exponent < -1023) {
    2477                exponent++;
    2478                mantissa >>>= 1;
    2479            }
    2480
    2481            long bits = mantissa | ((exponent + 1023L) << 52);
    2482            double x = Double.longBitsToDouble(bits);
    2483
    2484            if (negate) {
    2485                x = -x;
    2486            }
    2487
    2488            return x;
    2489
    2490        }
    2491
    2492        /** Convert the instance into a split double.
    2493         * @return an array of two doubles which sum represent the instance
    2494         * @see #toDouble()
    2495         */
    2496        public double[] toSplitDouble() {
    2497            double split[] = new double[2];
    2498            long mask = 0xffffffffc0000000L;
    2499
    2500            split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
    2501            split[1] = subtract(newInstance(split[0])).toDouble();
    2502
    2503            return split;
    2504        }
    2505
    2506    }