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