View Javadoc

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