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