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