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