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