Dfp.java

  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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.dfp;

  25. import java.util.Arrays;

  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.FieldSinhCosh;
  32. import org.hipparchus.util.MathUtils;

  34. /**
  35.  *  Decimal floating point library for Java
  36.  *
  37.  *  <p>Another floating point class.  This one is built using radix 10000
  38.  *  which is 10<sup>4</sup>, so its almost decimal.</p>
  39.  *
  40.  *  <p>The design goals here are:</p>
  41.  *  <ol>
  42.  *    <li>Decimal math, or close to it</li>
  43.  *    <li>Settable precision (but no mix between numbers using different settings)</li>
  44.  *    <li>Portability.  Code should be kept as portable as possible.</li>
  45.  *    <li>Performance</li>
  46.  *    <li>Accuracy  - Results should always be +/- 1 ULP for basic
  47.  *         algebraic operation</li>
  48.  *    <li>Comply with IEEE 854-1987 as much as possible.
  49.  *         (See IEEE 854-1987 notes below)</li>
  50.  *  </ol>
  51.  *
  52.  *  <p>Trade offs:</p>
  53.  *  <ol>
  54.  *    <li>Memory foot print.  I'm using more memory than necessary to
  55.  *         represent numbers to get better performance.</li>
  56.  *    <li>Digits are bigger, so rounding is a greater loss.  So, if you
  57.  *         really need 12 decimal digits, better use 4 base 10000 digits
  58.  *         there can be one partially filled.</li>
  59.  *  </ol>
  60.  *
  61.  *  <p>Numbers are represented  in the following form:
  62.  *  \[
  63.  *  n  =  \mathrm{sign} \times \mathrm{mant} \times \mathrm{radix}^\mathrm{exp}
  64.  *  \]
  65.  *  where sign is &plusmn;1, mantissa represents a fractional number between
  66.  *  zero and one.  mant[0] is the least significant digit.
  67.  *  exp is in the range of -32767 to 32768</p>
  68.  *
  69.  *  <p>IEEE 854-1987  Notes and differences</p>
  70.  *
  71.  *  <p>IEEE 854 requires the radix to be either 2 or 10.  The radix here is
  72.  *  10000, so that requirement is not met, but  it is possible that a
  73.  *  subclassed can be made to make it behave as a radix 10
  74.  *  number.  It is my opinion that if it looks and behaves as a radix
  75.  *  10 number then it is one and that requirement would be met.</p>
  76.  *
  77.  *  <p>The radix of 10000 was chosen because it should be faster to operate
  78.  *  on 4 decimal digits at once instead of one at a time.  Radix 10 behavior
  79.  *  can be realized by adding an additional rounding step to ensure that
  80.  *  the number of decimal digits represented is constant.</p>
  81.  *
  82.  *  <p>The IEEE standard specifically leaves out internal data encoding,
  83.  *  so it is reasonable to conclude that such a subclass of this radix
  84.  *  10000 system is merely an encoding of a radix 10 system.</p>
  85.  *
  86.  *  <p>IEEE 854 also specifies the existence of "sub-normal" numbers.  This
  87.  *  class does not contain any such entities.  The most significant radix
  88.  *  10000 digit is always non-zero.  Instead, we support "gradual underflow"
  89.  *  by raising the underflow flag for numbers less with exponent less than
  90.  *  expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
  91.  *  Thus the smallest number we can represent would be:
  92.  *  1E(-(MIN_EXP-digits-1)*4),  eg, for digits=5, MIN_EXP=-32767, that would
  93.  *  be 1e-131092.</p>
  94.  *
  95.  *  <p>IEEE 854 defines that the implied radix point lies just to the right
  96.  *  of the most significant digit and to the left of the remaining digits.
  97.  *  This implementation puts the implied radix point to the left of all
  98.  *  digits including the most significant one.  The most significant digit
  99.  *  here is the one just to the right of the radix point.  This is a fine
  100.  *  detail and is really only a matter of definition.  Any side effects of
  101.  *  this can be rendered invisible by a subclass.</p>
  102.  * @see DfpField
  103.  */
  104. public class Dfp implements CalculusFieldElement<Dfp> {

  106.     /** The radix, or base of this system.  Set to 10000 */
  107.     public static final int RADIX = 10000;

  109.     /** The minimum exponent before underflow is signaled.  Flush to zero
  110.      *  occurs at minExp-DIGITS */
  111.     public static final int MIN_EXP = -32767;

  113.     /** The maximum exponent before overflow is signaled and results flushed
  114.      *  to infinity */
  115.     public static final int MAX_EXP =  32768;

  117.     /** The amount under/overflows are scaled by before going to trap handler */
  118.     public static final int ERR_SCALE = 32760;

  120.     /** Indicator value for normal finite numbers. */
  121.     public static final byte FINITE = 0;

  123.     /** Indicator value for Infinity. */
  124.     public static final byte INFINITE = 1;

  126.     /** Indicator value for signaling NaN. */
  127.     public static final byte SNAN = 2;

  129.     /** Indicator value for quiet NaN. */
  130.     public static final byte QNAN = 3;

  132.     /** String for NaN representation. */
  133.     private static final String NAN_STRING = "NaN";

  135.     /** String for positive infinity representation. */
  136.     private static final String POS_INFINITY_STRING = "Infinity";

  138.     /** String for negative infinity representation. */
  139.     private static final String NEG_INFINITY_STRING = "-Infinity";

  141.     /** Name for traps triggered by addition. */
  142.     private static final String ADD_TRAP = "add";

  144.     /** Name for traps triggered by multiplication. */
  145.     private static final String MULTIPLY_TRAP = "multiply";

  147.     /** Name for traps triggered by division. */
  148.     private static final String DIVIDE_TRAP = "divide";

  150.     /** Name for traps triggered by square root. */
  151.     private static final String SQRT_TRAP = "sqrt";

  153.     /** Name for traps triggered by alignment. */
  154.     private static final String ALIGN_TRAP = "align";

  156.     /** Name for traps triggered by truncation. */
  157.     private static final String TRUNC_TRAP = "trunc";

  159.     /** Name for traps triggered by nextAfter. */
  160.     private static final String NEXT_AFTER_TRAP = "nextAfter";

  162.     /** Name for traps triggered by lessThan. */
  163.     private static final String LESS_THAN_TRAP = "lessThan";

  165.     /** Name for traps triggered by greaterThan. */
  166.     private static final String GREATER_THAN_TRAP = "greaterThan";

  168.     /** Name for traps triggered by newInstance. */
  169.     private static final String NEW_INSTANCE_TRAP = "newInstance";

  171.     /** Multiplication factor for number of digits used to compute linear combinations. */
  172.     private static final int LINEAR_COMBINATION_DIGITS_FACTOR = 2;

  174.     /** Mantissa. */
  175.     protected int[] mant;

  177.     /** Sign bit: 1 for positive, -1 for negative. */
  178.     protected byte sign;

  180.     /** Exponent. */
  181.     protected int exp;

  183.     /** Indicator for non-finite / non-number values. */
  184.     protected byte nans;

  186.     /** Factory building similar Dfp's. */
  187.     private final DfpField field;

  189.     /** Makes an instance with a value of zero.
  190.      * @param field field to which this instance belongs
  191.      */
  192.     protected Dfp(final DfpField field) {
  193.         mant = new int[field.getRadixDigits()];
  194.         sign = 1;
  195.         exp = 0;
  196.         nans = FINITE;
  197.         this.field = field;
  198.     }

  200.     /** Create an instance from a byte value.
  201.      * @param field field to which this instance belongs
  202.      * @param x value to convert to an instance
  203.      */
  204.     protected Dfp(final DfpField field, byte x) {
  205.         this(field, (long) x);
  206.     }

  208.     /** Create an instance from an int value.
  209.      * @param field field to which this instance belongs
  210.      * @param x value to convert to an instance
  211.      */
  212.     protected Dfp(final DfpField field, int x) {
  213.         this(field, (long) x);
  214.     }

  216.     /** Create an instance from a long value.
  217.      * @param field field to which this instance belongs
  218.      * @param x value to convert to an instance
  219.      */
  220.     protected Dfp(final DfpField field, long x) {

  222.         // initialize as if 0
  223.         mant = new int[field.getRadixDigits()];
  224.         nans = FINITE;
  225.         this.field = field;

  227.         boolean isLongMin = false;
  228.         if (x == Long.MIN_VALUE) {
  229.             // special case for Long.MIN_VALUE (-9223372036854775808)
  230.             // we must shift it before taking its absolute value
  231.             isLongMin = true;
  232.             ++x;
  233.         }

  235.         // set the sign
  236.         if (x < 0) {
  237.             sign = -1;
  238.             x = -x;
  239.         } else {
  240.             sign = 1;
  241.         }

  243.         exp = 0;
  244.         while (x != 0) {
  245.             System.arraycopy(mant, mant.length - exp, mant, mant.length - 1 - exp, exp);
  246.             mant[mant.length - 1] = (int) (x % RADIX);
  247.             x /= RADIX;
  248.             exp++;
  249.         }

  251.         if (isLongMin) {
  252.             // remove the shift added for Long.MIN_VALUE
  253.             // we know in this case that fixing the last digit is sufficient
  254.             for (int i = 0; i < mant.length - 1; i++) {
  255.                 if (mant[i] != 0) {
  256.                     mant[i]++;
  257.                     break;
  258.                 }
  259.             }
  260.         }
  261.     }

  263.     /** Create an instance from a double value.
  264.      * @param field field to which this instance belongs
  265.      * @param x value to convert to an instance
  266.      */
  267.     protected Dfp(final DfpField field, double x) {

  269.         // initialize as if 0
  270.         mant = new int[field.getRadixDigits()];
  271.         this.field = field;

  273.         long bits = Double.doubleToLongBits(x);
  274.         long mantissa = bits & 0x000fffffffffffffL;
  275.         int exponent = (int) ((bits & 0x7ff0000000000000L) >> 52) - 1023;

  277.         if (exponent == -1023) {
  278.             // Zero or sub-normal
  279.             if (x == 0) {
  280.                 // make sure 0 has the right sign
  281.                 if ((bits & 0x8000000000000000L) != 0) {
  282.                     sign = -1;
  283.                 } else {
  284.                     sign = 1;
  285.                 }
  286.                 return;
  287.             }

  289.             exponent++;

  291.             // Normalize the subnormal number
  292.             while ( (mantissa & 0x0010000000000000L) == 0) {
  293.                 exponent--;
  294.                 mantissa <<= 1;
  295.             }
  296.             mantissa &= 0x000fffffffffffffL;
  297.         }

  299.         if (exponent == 1024) {
  300.             // infinity or NAN
  301.             if (x != x) {
  302.                 sign = (byte) 1;
  303.                 nans = QNAN;
  304.             } else if (x < 0) {
  305.                 sign = (byte) -1;
  306.                 nans = INFINITE;
  307.             } else {
  308.                 sign = (byte) 1;
  309.                 nans = INFINITE;
  310.             }
  311.             return;
  312.         }

  314.         Dfp xdfp = new Dfp(field, mantissa);
  315.         xdfp = xdfp.divide(new Dfp(field, 4503599627370496L)).add(field.getOne());  // Divide by 2^52, then add one
  316.         xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));

  318.         if ((bits & 0x8000000000000000L) != 0) {
  319.             xdfp = xdfp.negate();
  320.         }

  322.         System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
  323.         sign = xdfp.sign;
  324.         exp  = xdfp.exp;
  325.         nans = xdfp.nans;

  327.     }

  329.     /** Copy constructor.
  330.      * @param d instance to copy
  331.      */
  332.     public Dfp(final Dfp d) {
  333.         mant  = d.mant.clone();
  334.         sign  = d.sign;
  335.         exp   = d.exp;
  336.         nans  = d.nans;
  337.         field = d.field;
  338.     }

  340.     /** Create an instance from a String representation.
  341.      * @param field field to which this instance belongs
  342.      * @param s string representation of the instance
  343.      */
  344.     protected Dfp(final DfpField field, final String s) {

  346.         // initialize as if 0
  347.         mant = new int[field.getRadixDigits()];
  348.         sign = 1;
  349.         nans = FINITE;
  350.         this.field = field;

  352.         boolean decimalFound = false;
  353.         final int rsize = 4;   // size of radix in decimal digits
  354.         final int offset = 4;  // Starting offset into Striped
  355.         final char[] striped = new char[getRadixDigits() * rsize + offset * 2];

  357.         // Check some special cases
  358.         if (POS_INFINITY_STRING.equals(s)) {
  359.             sign = (byte) 1;
  360.             nans = INFINITE;
  361.             return;
  362.         }

  364.         if (NEG_INFINITY_STRING.equals(s)) {
  365.             sign = (byte) -1;
  366.             nans = INFINITE;
  367.             return;
  368.         }

  370.         if (NAN_STRING.equals(s)) {
  371.             sign = (byte) 1;
  372.             nans = QNAN;
  373.             return;
  374.         }

  376.         // Check for scientific notation
  377.         int p = s.indexOf('e');
  378.         if (p == -1) { // try upper case?
  379.             p = s.indexOf('E');
  380.         }

  382.         final String fpdecimal;
  383.         int sciexp = 0;
  384.         if (p != -1) {
  385.             // scientific notation
  386.             fpdecimal = s.substring(0, p);
  387.             String fpexp = s.substring(p+1);
  388.             boolean negative = false;

  390.             for (int i=0; i<fpexp.length(); i++)
  391.             {
  392.                 if (fpexp.charAt(i) == '-')
  393.                 {
  394.                     negative = true;
  395.                     continue;
  396.                 }
  397.                 if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') {
  398.                     sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
  399.                 }
  400.             }

  402.             if (negative) {
  403.                 sciexp = -sciexp;
  404.             }
  405.         } else {
  406.             // normal case
  407.             fpdecimal = s;
  408.         }

  410.         // If there is a minus sign in the number then it is negative
  411.         if (fpdecimal.indexOf('-') !=  -1) {
  412.             sign = -1;
  413.         }

  415.         // First off, find all of the leading zeros, trailing zeros, and significant digits
  416.         p = 0;

  418.         // Move p to first significant digit
  419.         int decimalPos = 0;
  420.         for (;;) {
  421.             if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
  422.                 break;
  423.             }

  425.             if (decimalFound && fpdecimal.charAt(p) == '0') {
  426.                 decimalPos--;
  427.             }

  429.             if (fpdecimal.charAt(p) == '.') {
  430.                 decimalFound = true;
  431.             }

  433.             p++;

  435.             if (p == fpdecimal.length()) {
  436.                 break;
  437.             }
  438.         }

  440.         // Copy the string onto Stripped
  441.         int q = offset;
  442.         striped[0] = '0';
  443.         striped[1] = '0';
  444.         striped[2] = '0';
  445.         striped[3] = '0';
  446.         int significantDigits=0;
  447.         for(;;) {
  448.             if (p == (fpdecimal.length())) {
  449.                 break;
  450.             }

  452.             // Don't want to run pass the end of the array
  453.             if (q == mant.length*rsize+offset+1) {
  454.                 break;
  455.             }

  457.             if (fpdecimal.charAt(p) == '.') {
  458.                 decimalFound = true;
  459.                 decimalPos = significantDigits;
  460.                 p++;
  461.                 continue;
  462.             }

  464.             if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
  465.                 p++;
  466.                 continue;
  467.             }

  469.             striped[q] = fpdecimal.charAt(p);
  470.             q++;
  471.             p++;
  472.             significantDigits++;
  473.         }


  476.         // If the decimal point has been found then get rid of trailing zeros.
  477.         if (decimalFound && q != offset) {
  478.             for (;;) {
  479.                 q--;
  480.                 if (q == offset) {
  481.                     break;
  482.                 }
  483.                 if (striped[q] == '0') {
  484.                     significantDigits--;
  485.                 } else {
  486.                     break;
  487.                 }
  488.             }
  489.         }

  491.         // special case of numbers like "0.00000"
  492.         if (decimalFound && significantDigits == 0) {
  493.             decimalPos = 0;
  494.         }

  496.         // Implicit decimal point at end of number if not present
  497.         if (!decimalFound) {
  498.             decimalPos = q-offset;
  499.         }

  501.         // Find the number of significant trailing zeros
  502.         q = offset;  // set q to point to first sig digit
  503.         p = significantDigits-1+offset;

  505.         while (p > q) {
  506.             if (striped[p] != '0') {
  507.                 break;
  508.             }
  509.             p--;
  510.         }

  512.         // Make sure the decimal is on a mod 10000 boundary
  513.         int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
  514.         q -= i;
  515.         decimalPos += i;

  517.         // Make the mantissa length right by adding zeros at the end if necessary
  518.         while ((p - q) < (mant.length * rsize)) {
  519.             for (i = 0; i < rsize; i++) {
  520.                 striped[++p] = '0';
  521.             }
  522.         }

  524.         // Ok, now we know how many trailing zeros there are,
  525.         // and where the least significant digit is
  526.         for (i = mant.length - 1; i >= 0; i--) {
  527.             mant[i] = (striped[q]   - '0') * 1000 +
  528.                       (striped[q+1] - '0') * 100  +
  529.                       (striped[q+2] - '0') * 10   +
  530.                       (striped[q+3] - '0');
  531.             q += 4;
  532.         }

  534.         exp = (decimalPos+sciexp) / rsize;

  536.         if (q < striped.length) {
  537.             // Is there possible another digit?
  538.             round((striped[q] - '0')*1000);
  539.         }

  541.     }

  543.     /** Creates an instance with a non-finite value.
  544.      * @param field field to which this instance belongs
  545.      * @param sign sign of the Dfp to create
  546.      * @param nans code of the value, must be one of {@link #INFINITE},
  547.      * {@link #SNAN},  {@link #QNAN}
  548.      */
  549.     protected Dfp(final DfpField field, final byte sign, final byte nans) {
  550.         this.field = field;
  551.         this.mant    = new int[field.getRadixDigits()];
  552.         this.sign    = sign;
  553.         this.exp     = 0;
  554.         this.nans    = nans;
  555.     }

  557.     /** Create an instance with a value of 0.
  558.      * Use this internally in preference to constructors to facilitate subclasses
  559.      * @return a new instance with a value of 0
  560.      */
  561.     public Dfp newInstance() {
  562.         return new Dfp(getField());
  563.     }

  565.     /** Create an instance from a byte value.
  566.      * @param x value to convert to an instance
  567.      * @return a new instance with value x
  568.      */
  569.     public Dfp newInstance(final byte x) {
  570.         return new Dfp(getField(), x);
  571.     }

  573.     /** Create an instance from an int value.
  574.      * @param x value to convert to an instance
  575.      * @return a new instance with value x
  576.      */
  577.     public Dfp newInstance(final int x) {
  578.         return new Dfp(getField(), x);
  579.     }

  581.     /** Create an instance from a long value.
  582.      * @param x value to convert to an instance
  583.      * @return a new instance with value x
  584.      */
  585.     public Dfp newInstance(final long x) {
  586.         return new Dfp(getField(), x);
  587.     }

  589.     /** {@inheritDoc} */
  590.     @Override
  591.     public Dfp newInstance(final double x) {
  592.         return new Dfp(getField(), x);
  593.     }

  595.     /** Create an instance by copying an existing one.
  596.      * Use this internally in preference to constructors to facilitate subclasses.
  597.      * @param d instance to copy
  598.      * @return a new instance with the same value as d
  599.      */
  600.     public Dfp newInstance(final Dfp d) {

  602.         // make sure we don't mix number with different precision
  603.         if (field.getRadixDigits() != d.field.getRadixDigits()) {
  604.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  605.             final Dfp result = newInstance(getZero());
  606.             result.nans = QNAN;
  607.             return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
  608.         }

  610.         return new Dfp(d);

  612.     }

  614.     /** Create an instance from a String representation.
  615.      * Use this internally in preference to constructors to facilitate subclasses.
  616.      * @param s string representation of the instance
  617.      * @return a new instance parsed from specified string
  618.      */
  619.     public Dfp newInstance(final String s) {
  620.         return new Dfp(field, s);
  621.     }

  623.     /** Creates an instance with a non-finite value.
  624.      * @param sig sign of the Dfp to create
  625.      * @param code code of the value, must be one of {@link #INFINITE},
  626.      * {@link #SNAN},  {@link #QNAN}
  627.      * @return a new instance with a non-finite value
  628.      */
  629.     public Dfp newInstance(final byte sig, final byte code) {
  630.         return field.newDfp(sig, code);
  631.     }

  633.     /** Creates an instance by converting the instance to a different field (i.e. different accuracy).
  634.      * <p>
  635.      * If the target field as a greater number of digits, the extra least significant digits
  636.      * will be set to zero.
  637.      * </p>
  638.      * @param targetField field to convert the instance to
  639.      * @param rmode rounding mode to use if target field as less digits than the instance, can be null otherwise
  640.      * @return converted instance (or the instance itself if it already has the required number of digits)
  641.      * @see DfpField#getExtendedField(int, boolean)
  642.      * @since 1.7
  643.      */
  644.     public Dfp newInstance(final DfpField targetField, final DfpField.RoundingMode rmode) {
  645.         final int deltaLength = targetField.getRadixDigits() - field.getRadixDigits();
  646.         if (deltaLength == 0) {

  648.             // no conversion, we return the instance itself
  649.             return this;

  651.         } else {

  653.             // create an instance (initially set to 0) with the expected number of digits
  654.             Dfp result = new Dfp(targetField);
  655.             result.sign = sign;
  656.             result.exp  = exp;
  657.             result.nans = nans;
  658.             if (nans == 0) {

  660.                 if (deltaLength < 0) {

  662.                     // copy only the most significant digits, dropping the least significant ones
  663.                     // the result corresponds to pure truncation, proper rounding will follow
  664.                     System.arraycopy(mant, -deltaLength, result.mant, 0, result.mant.length);

  666.                     // check if we have dropped any non-zero digits in the low part
  667.                     // (not counting the last dropped digit which will be handled specially)
  668.                     final int last = -(deltaLength + 1);
  669.                     boolean zeroLSB = true;
  670.                     for (int i = 0; i < last; ++i) {
  671.                         zeroLSB &= mant[i] == 0;
  672.                     }

  674.                     if (!(zeroLSB && mant[last] == 0)) {
  675.                         // there are some non-zero digits that have been discarded, perform rounding

  677.                         if (shouldIncrement(rmode, zeroLSB, mant[last], result.mant[0], sign)) {
  678.                             // rounding requires incrementing the mantissa
  679.                             result.incrementMantissa();
  680.                         }

  682.                         targetField.setIEEEFlagsBits(DfpField.FLAG_INEXACT);  // signal inexact
  683.                         result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);

  685.                     }

  687.                 } else {
  688.                     // copy all digits as the new most significant ones, leaving the least significant digits to zero
  689.                     System.arraycopy(mant, 0, result.mant, deltaLength, mant.length);
  690.                 }

  692.             }

  694.             return result;

  696.         }
  697.     }

  699.     /** Check if mantissa of a truncated number must be incremented.
  700.      * <p>
  701.      * This method must be called <em>only</em> when some non-zero digits have been
  702.      * discarded (i.e. when either {@code zeroLSB} is false or {@code lastDiscarded} is non-zero),
  703.      * otherwise it would generate false positive
  704.      * </p>
  705.      * @param rmode rounding mode to use if target field as less digits than the instance, can be null otherwise
  706.      * @param zeroLSB true is least significant discarded digits (except last) are all zero
  707.      * @param lastDiscarded last discarded digit
  708.      * @param firstNonDiscarded first non-discarded digit
  709.      * @param sign of the number
  710.      * @return true if the already truncated mantissa should be incremented to achieve correct rounding
  711.      * @since 1.7
  712.      */
  713.     private static boolean shouldIncrement(final DfpField.RoundingMode rmode,
  714.                                            final boolean zeroLSB, final int lastDiscarded,
  715.                                            final int firstNonDiscarded, final int sign) {
  716.         switch (rmode) {
  717.             case ROUND_DOWN :
  718.                 return false;

  720.             case ROUND_UP :
  721.                 return true;

  723.             case ROUND_HALF_UP :
  724.                 return lastDiscarded >= 5000;

  726.             case ROUND_HALF_DOWN :
  727.                 return isAboveHalfWay(zeroLSB, lastDiscarded);

  729.             case ROUND_HALF_EVEN :
  730.                 return (isHalfWay(zeroLSB, lastDiscarded) && (firstNonDiscarded & 0x1) == 0x1) ||
  731.                        isAboveHalfWay(zeroLSB, lastDiscarded);

  733.             case ROUND_HALF_ODD :
  734.                 return (isHalfWay(zeroLSB, lastDiscarded) && (firstNonDiscarded & 0x1) == 0x0) ||
  735.                        isAboveHalfWay(zeroLSB, lastDiscarded);

  737.             case ROUND_CEIL :
  738.                 return sign > 0;

  740.             case ROUND_FLOOR :
  741.                 return sign < 0;

  743.             default :
  744.                 // this should never happen
  745.                 throw MathRuntimeException.createInternalError();
  746.         }
  747.     }

  749.     /** Increment the mantissa of the instance
  750.      * @since 1.7
  751.      */
  752.     private void incrementMantissa() {
  753.         boolean carry = true;
  754.         for (int i = 0; carry && i < mant.length; ++i) {
  755.             ++mant[i];
  756.             if (mant[i] >= RADIX) {
  757.                 mant[i] -= RADIX;
  758.             } else {
  759.                 carry = false;
  760.             }
  761.         }
  762.         if (carry) {
  763.             // we have exceeded capacity, we need to drop one digit
  764.             for (int i = 0; i < mant.length - 1; i++) {
  765.                 mant[i] = mant[i+1];
  766.             }
  767.             mant[mant.length - 1] = 1;
  768.             exp++;
  769.         }
  770.     }

  772.     /** Check if discarded digits are exactly halfway between two rounder numbers.
  773.      * @param zeroLSB true is least significant discarded digits (except last) are all zero
  774.      * @param lastDiscarded last discarded digit
  775.      * @return true if discarded digits correspond to a number exactly halfway between two rounded numbers
  776.      * @since 1.7
  777.      */
  778.     private static boolean isHalfWay(final boolean zeroLSB, final int lastDiscarded) {
  779.         return lastDiscarded == 5000 && zeroLSB;
  780.     }

  782.     /** Check if discarded digits are strictly above halfway between two rounder numbers.
  783.      * @param zeroLSB true is least significant discarded digits (except last) are all zero
  784.      * @param lastDiscarded last discarded digit
  785.      * @return true if discarded digits correspond to a number strictly above halfway between two rounded numbers
  786.      * @since 1.7
  787.      */
  788.     private static boolean isAboveHalfWay(final boolean zeroLSB, final int lastDiscarded) {
  789.         return (lastDiscarded > 5000) || (lastDiscarded == 5000 && !zeroLSB);
  790.     }

  792.     /** Get the {@link org.hipparchus.Field Field} (really a {@link DfpField}) to which the instance belongs.
  793.      * <p>
  794.      * The field is linked to the number of digits and acts as a factory
  795.      * for {@link Dfp} instances.
  796.      * </p>
  797.      * @return {@link org.hipparchus.Field Field} (really a {@link DfpField}) to which the instance belongs
  798.      */
  799.     @Override
  800.     public DfpField getField() {
  801.         return field;
  802.     }

  804.     /** Get the number of radix digits of the instance.
  805.      * @return number of radix digits
  806.      */
  807.     public int getRadixDigits() {
  808.         return field.getRadixDigits();
  809.     }

  811.     /** Get the constant 0.
  812.      * @return a Dfp with value zero
  813.      */
  814.     public Dfp getZero() {
  815.         return field.getZero();
  816.     }

  818.     /** Get the constant 1.
  819.      * @return a Dfp with value one
  820.      */
  821.     public Dfp getOne() {
  822.         return field.getOne();
  823.     }

  825.     /** Get the constant 2.
  826.      * @return a Dfp with value two
  827.      */
  828.     public Dfp getTwo() {
  829.         return field.getTwo();
  830.     }

  832.     /** Shift the mantissa left, and adjust the exponent to compensate.
  833.      */
  834.     protected void shiftLeft() {
  835.         for (int i = mant.length - 1; i > 0; i--) {
  836.             mant[i] = mant[i-1];
  837.         }
  838.         mant[0] = 0;
  839.         exp--;
  840.     }

  842.     /* Note that shiftRight() does not call round() as that round() itself
  843.      uses shiftRight() */
  844.     /** Shift the mantissa right, and adjust the exponent to compensate.
  845.      */
  846.     protected void shiftRight() {
  847.         for (int i = 0; i < mant.length - 1; i++) {
  848.             mant[i] = mant[i+1];
  849.         }
  850.         mant[mant.length - 1] = 0;
  851.         exp++;
  852.     }

  854.     /** Make our exp equal to the supplied one, this may cause rounding.
  855.      *  Also causes de-normalized numbers.  These numbers are generally
  856.      *  dangerous because most routines assume normalized numbers.
  857.      *  Align doesn't round, so it will return the last digit destroyed
  858.      *  by shifting right.
  859.      *  @param e desired exponent
  860.      *  @return last digit destroyed by shifting right
  861.      */
  862.     protected int align(int e) {
  863.         int lostdigit = 0;
  864.         boolean inexact = false;

  866.         int diff = exp - e;

  868.         int adiff = diff;
  869.         if (adiff < 0) {
  870.             adiff = -adiff;
  871.         }

  873.         if (diff == 0) {
  874.             return 0;
  875.         }

  877.         if (adiff > (mant.length + 1)) {
  878.             // Special case
  879.             Arrays.fill(mant, 0);
  880.             exp = e;

  882.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
  883.             dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);

  885.             return 0;
  886.         }

  888.         for (int i = 0; i < adiff; i++) {
  889.             if (diff < 0) {
  890.                 /* Keep track of loss -- only signal inexact after losing 2 digits.
  891.                  * the first lost digit is returned to add() and may be incorporated
  892.                  * into the result.
  893.                  */
  894.                 if (lostdigit != 0) {
  895.                     inexact = true;
  896.                 }

  898.                 lostdigit = mant[0];

  900.                 shiftRight();
  901.             } else {
  902.                 shiftLeft();
  903.             }
  904.         }

  906.         if (inexact) {
  907.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
  908.             dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
  909.         }

  911.         return lostdigit;

  913.     }

  915.     /** Check if instance is less than x.
  916.      * @param x number to check instance against
  917.      * @return true if instance is less than x and neither are NaN, false otherwise
  918.      */
  919.     public boolean lessThan(final Dfp x) {

  921.         // make sure we don't mix number with different precision
  922.         if (field.getRadixDigits() != x.field.getRadixDigits()) {
  923.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  924.             final Dfp result = newInstance(getZero());
  925.             result.nans = QNAN;
  926.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
  927.             return false;
  928.         }

  930.         /* if a nan is involved, signal invalid and return false */
  931.         if (isNaN() || x.isNaN()) {
  932.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  933.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
  934.             return false;
  935.         }

  937.         return compare(this, x) < 0;
  938.     }

  940.     /** Check if instance is greater than x.
  941.      * @param x number to check instance against
  942.      * @return true if instance is greater than x and neither are NaN, false otherwise
  943.      */
  944.     public boolean greaterThan(final Dfp x) {

  946.         // make sure we don't mix number with different precision
  947.         if (field.getRadixDigits() != x.field.getRadixDigits()) {
  948.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  949.             final Dfp result = newInstance(getZero());
  950.             result.nans = QNAN;
  951.             dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
  952.             return false;
  953.         }

  955.         /* if a nan is involved, signal invalid and return false */
  956.         if (isNaN() || x.isNaN()) {
  957.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  958.             dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
  959.             return false;
  960.         }

  962.         return compare(this, x) > 0;
  963.     }

  965.     /** Check if instance is less than or equal to 0.
  966.      * @return true if instance is not NaN and less than or equal to 0, false otherwise
  967.      */
  968.     public boolean negativeOrNull() {

  970.         if (isNaN()) {
  971.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  972.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
  973.             return false;
  974.         }

  976.         return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite());

  978.     }

  980.     /** Check if instance is strictly less than 0.
  981.      * @return true if instance is not NaN and less than or equal to 0, false otherwise
  982.      */
  983.     public boolean strictlyNegative() {

  985.         if (isNaN()) {
  986.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  987.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
  988.             return false;
  989.         }

  991.         return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite());

  993.     }

  995.     /** Check if instance is greater than or equal to 0.
  996.      * @return true if instance is not NaN and greater than or equal to 0, false otherwise
  997.      */
  998.     public boolean positiveOrNull() {

  1000.         if (isNaN()) {
  1001.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1002.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
  1003.             return false;
  1004.         }

  1006.         return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite());

  1008.     }

  1010.     /** Check if instance is strictly greater than 0.
  1011.      * @return true if instance is not NaN and greater than or equal to 0, false otherwise
  1012.      */
  1013.     public boolean strictlyPositive() {

  1015.         if (isNaN()) {
  1016.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1017.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
  1018.             return false;
  1019.         }

  1021.         return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite());

  1023.     }

  1025.     /** {@inheritDoc} */
  1026.     @Override
  1027.     public Dfp abs() {
  1028.         Dfp result = newInstance(this);
  1029.         result.sign = 1;
  1030.         return result;
  1031.     }

  1033.     /** {@inheritDoc} */
  1034.     @Override
  1035.     public boolean isInfinite() {
  1036.         return nans == INFINITE;
  1037.     }

  1039.     /** {@inheritDoc} */
  1040.     @Override
  1041.     public boolean isNaN() {
  1042.         return (nans == QNAN) || (nans == SNAN);
  1043.     }

  1045.     /** Check if instance is equal to zero.
  1046.      * @return true if instance is equal to zero
  1047.      */
  1048.     @Override
  1049.     public boolean isZero() {

  1051.         if (isNaN()) {
  1052.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1053.             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
  1054.             return false;
  1055.         }

  1057.         return (mant[mant.length - 1] == 0) && !isInfinite();

  1059.     }

  1061.     /** Check if instance is equal to x.
  1062.      * @param other object to check instance against
  1063.      * @return true if instance is equal to x and neither are NaN, false otherwise
  1064.      */
  1065.     @Override
  1066.     public boolean equals(final Object other) {

  1068.         if (other instanceof Dfp) {
  1069.             final Dfp x = (Dfp) other;
  1070.             if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
  1071.                 return false;
  1072.             }

  1074.             return compare(this, x) == 0;
  1075.         }

  1077.         return false;

  1079.     }

  1081.     /**
  1082.      * Gets a hashCode for the instance.
  1083.      * @return a hash code value for this object
  1084.      */
  1085.     @Override
  1086.     public int hashCode() {
  1087.         return 17 + (isZero() ? 0 : (sign << 8)) + (nans << 16) + exp + Arrays.hashCode(mant);
  1088.     }

  1090.     /** Check if instance is not equal to x.
  1091.      * @param x number to check instance against
  1092.      * @return true if instance is not equal to x and neither are NaN, false otherwise
  1093.      */
  1094.     public boolean unequal(final Dfp x) {
  1095.         if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
  1096.             return false;
  1097.         }

  1099.         return greaterThan(x) || lessThan(x);
  1100.     }

  1102.     /** Compare two instances.
  1103.      * @param a first instance in comparison
  1104.      * @param b second instance in comparison
  1105.      * @return -1 if a<b, 1 if a>b and 0 if a==b
  1106.      *  Note this method does not properly handle NaNs or numbers with different precision.
  1107.      */
  1108.     private static int compare(final Dfp a, final Dfp b) {
  1109.         // Ignore the sign of zero
  1110.         if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
  1111.             a.nans == FINITE && b.nans == FINITE) {
  1112.             return 0;
  1113.         }

  1115.         if (a.sign != b.sign) {
  1116.             if (a.sign == -1) {
  1117.                 return -1;
  1118.             } else {
  1119.                 return 1;
  1120.             }
  1121.         }

  1123.         // deal with the infinities
  1124.         if (a.nans == INFINITE && b.nans == FINITE) {
  1125.             return a.sign;
  1126.         }

  1128.         if (a.nans == FINITE && b.nans == INFINITE) {
  1129.             return -b.sign;
  1130.         }

  1132.         if (a.nans == INFINITE && b.nans == INFINITE) {
  1133.             return 0;
  1134.         }

  1136.         // Handle special case when a or b is zero, by ignoring the exponents
  1137.         if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) {
  1138.             if (a.exp < b.exp) {
  1139.                 return -a.sign;
  1140.             }

  1142.             if (a.exp > b.exp) {
  1143.                 return a.sign;
  1144.             }
  1145.         }

  1147.         // compare the mantissas
  1148.         for (int i = a.mant.length - 1; i >= 0; i--) {
  1149.             if (a.mant[i] > b.mant[i]) {
  1150.                 return a.sign;
  1151.             }

  1153.             if (a.mant[i] < b.mant[i]) {
  1154.                 return -a.sign;
  1155.             }
  1156.         }

  1158.         return 0;

  1160.     }

  1162.     /** Round to nearest integer using the round-half-even method.
  1163.      *  That is round to nearest integer unless both are equidistant.
  1164.      *  In which case round to the even one.
  1165.      *  @return rounded value
  1166.      */
  1167.     @Override
  1168.     public Dfp rint() {
  1169.         return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
  1170.     }

  1172.     /** Round to an integer using the round floor mode.
  1173.      * That is, round toward -Infinity
  1174.      *  @return rounded value
  1175.      */
  1176.     @Override
  1177.     public Dfp floor() {
  1178.         return trunc(DfpField.RoundingMode.ROUND_FLOOR);
  1179.     }

  1181.     /** Round to an integer using the round ceil mode.
  1182.      * That is, round toward +Infinity
  1183.      *  @return rounded value
  1184.      */
  1185.     @Override
  1186.     public Dfp ceil() {
  1187.         return trunc(DfpField.RoundingMode.ROUND_CEIL);
  1188.     }

  1190.     /** Returns the IEEE remainder.
  1191.      * @param d divisor
  1192.      * @return this less n &times; d, where n is the integer closest to this/d
  1193.      */
  1194.     @Override
  1195.     public Dfp remainder(final Dfp d) {

  1197.         final Dfp result = this.subtract(this.divide(d).rint().multiply(d));

  1199.         // IEEE 854-1987 says that if the result is zero, then it carries the sign of this
  1200.         if (result.mant[mant.length-1] == 0) {
  1201.             result.sign = sign;
  1202.         }

  1204.         return result;

  1206.     }

  1208.     /** Does the integer conversions with the specified rounding.
  1209.      * @param rmode rounding mode to use
  1210.      * @return truncated value
  1211.      */
  1212.     protected Dfp trunc(final DfpField.RoundingMode rmode) {
  1213.         boolean changed = false;

  1215.         if (isNaN()) {
  1216.             return newInstance(this);
  1217.         }

  1219.         if (nans == INFINITE) {
  1220.             return newInstance(this);
  1221.         }

  1223.         if (mant[mant.length-1] == 0) {
  1224.             // a is zero
  1225.             return newInstance(this);
  1226.         }

  1228.         /* If the exponent is less than zero then we can certainly
  1229.          * return -1, 0 or +1 depending on sign and rounding mode */
  1230.         if (exp < 0) {
  1231.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
  1232.             final Dfp result;
  1233.             if (sign == -1 && rmode == DfpField.RoundingMode.ROUND_FLOOR) {
  1234.                 result = newInstance(-1);
  1235.             } else if (sign == +1 && rmode == DfpField.RoundingMode.ROUND_CEIL) {
  1236.                 result = newInstance(+1);
  1237.             } else {
  1238.                 // for all other combinations of sign and mode, zero is the correct rounding
  1239.                 result = newInstance(0);
  1240.             }
  1241.             return dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
  1242.         }

  1244.         /* If the exponent is greater than or equal to digits, then it
  1245.          * must already be an integer since there is no precision left
  1246.          * for any fractional part */

  1248.         if (exp >= mant.length) {
  1249.             return newInstance(this);
  1250.         }

  1252.         /* General case:  create another dfp, result, that contains the
  1253.          * a with the fractional part lopped off.  */

  1255.         Dfp result = newInstance(this);
  1256.         for (int i = 0; i < mant.length-result.exp; i++) {
  1257.             changed |= result.mant[i] != 0;
  1258.             result.mant[i] = 0;
  1259.         }

  1261.         if (changed) {
  1262.             switch (rmode) {
  1263.                 case ROUND_FLOOR:
  1264.                     if (result.sign == -1) {
  1265.                         // then we must increment the mantissa by one
  1266.                         result = result.add(newInstance(-1));
  1267.                     }
  1268.                     break;

  1270.                 case ROUND_CEIL:
  1271.                     if (result.sign == 1) {
  1272.                         // then we must increment the mantissa by one
  1273.                         result = result.add(getOne());
  1274.                     }
  1275.                     break;

  1277.                 case ROUND_HALF_EVEN:
  1278.                 default:
  1279.                     final Dfp half = newInstance("0.5");
  1280.                     Dfp a = subtract(result);  // difference between this and result
  1281.                     a.sign = 1;            // force positive (take abs)
  1282.                     if (a.greaterThan(half)) {
  1283.                         a = newInstance(getOne());
  1284.                         a.sign = sign;
  1285.                         result = result.add(a);
  1286.                     }

  1288.                     // If exactly equal to 1/2 and odd then increment
  1289.                     if (a.equals(half) && result.exp > 0 && (result.mant[mant.length-result.exp]&1) != 0) {
  1290.                         a = newInstance(getOne());
  1291.                         a.sign = sign;
  1292.                         result = result.add(a);
  1293.                     }
  1294.                     break;
  1295.             }

  1297.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);  // signal inexact
  1298.             result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
  1299.             return result;
  1300.         }

  1302.         return result;
  1303.     }

  1305.     /** Convert this to an integer.
  1306.      * If greater than 2147483647, it returns 2147483647. If less than -2147483648 it returns -2147483648.
  1307.      * @return converted number
  1308.      */
  1309.     public int intValue() {
  1310.         Dfp rounded;
  1311.         int result = 0;

  1313.         rounded = rint();

  1315.         if (rounded.greaterThan(newInstance(2147483647))) {
  1316.             return 2147483647;
  1317.         }

  1319.         if (rounded.lessThan(newInstance(-2147483648))) {
  1320.             return -2147483648;
  1321.         }

  1323.         for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
  1324.             result = result * RADIX + rounded.mant[i];
  1325.         }

  1327.         if (rounded.sign == -1) {
  1328.             result = -result;
  1329.         }

  1331.         return result;
  1332.     }

  1334.     /** Get the exponent of the greatest power of 10000 that is
  1335.      *  less than or equal to the absolute value of this.  I.E.  if
  1336.      *  this is 10<sup>6</sup> then log10K would return 1.
  1337.      *  @return integer base 10000 logarithm
  1338.      */
  1339.     public int log10K() {
  1340.         return exp - 1;
  1341.     }

  1343.     /** Get the specified  power of 10000.
  1344.      * @param e desired power
  1345.      * @return 10000<sup>e</sup>
  1346.      */
  1347.     public Dfp power10K(final int e) {
  1348.         Dfp d = newInstance(getOne());
  1349.         d.exp = e + 1;
  1350.         return d;
  1351.     }

  1353.     /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
  1354.      *  @return integer base 10 logarithm
  1355.      */
  1356.     public int intLog10()  {
  1357.         if (mant[mant.length-1] > 1000) {
  1358.             return exp * 4 - 1;
  1359.         }
  1360.         if (mant[mant.length-1] > 100) {
  1361.             return exp * 4 - 2;
  1362.         }
  1363.         if (mant[mant.length-1] > 10) {
  1364.             return exp * 4 - 3;
  1365.         }
  1366.         return exp * 4 - 4;
  1367.     }

  1369.     /** Return the specified  power of 10.
  1370.      * @param e desired power
  1371.      * @return 10<sup>e</sup>
  1372.      */
  1373.     public Dfp power10(final int e) {
  1374.         Dfp d = newInstance(getOne());

  1376.         if (e >= 0) {
  1377.             d.exp = e / 4 + 1;
  1378.         } else {
  1379.             d.exp = (e + 1) / 4;
  1380.         }

  1382.         switch ((e % 4 + 4) % 4) {
  1383.             case 0:
  1384.                 break;
  1385.             case 1:
  1386.                 d = d.multiply(10);
  1387.                 break;
  1388.             case 2:
  1389.                 d = d.multiply(100);
  1390.                 break;
  1391.             default:
  1392.                 d = d.multiply(1000);
  1393.                 break;
  1394.         }

  1396.         return d;
  1397.     }

  1399.     /** Negate the mantissa of this by computing the complement.
  1400.      *  Leaves the sign bit unchanged, used internally by add.
  1401.      *  Denormalized numbers are handled properly here.
  1402.      *  @param extra ???
  1403.      *  @return ???
  1404.      */
  1405.     protected int complement(int extra) {

  1407.         extra = RADIX-extra;
  1408.         for (int i = 0; i < mant.length; i++) {
  1409.             mant[i] = RADIX-mant[i]-1;
  1410.         }

  1412.         int rh = extra / RADIX;
  1413.         extra -= rh * RADIX;
  1414.         for (int i = 0; i < mant.length; i++) {
  1415.             final int r = mant[i] + rh;
  1416.             rh = r / RADIX;
  1417.             mant[i] = r - rh * RADIX;
  1418.         }

  1420.         return extra;
  1421.     }

  1423.     /** Add x to this.
  1424.      * @param x number to add
  1425.      * @return sum of this and x
  1426.      */
  1427.     @Override
  1428.     public Dfp add(final Dfp x) {

  1430.         // make sure we don't mix number with different precision
  1431.         if (field.getRadixDigits() != x.field.getRadixDigits()) {
  1432.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1433.             final Dfp result = newInstance(getZero());
  1434.             result.nans = QNAN;
  1435.             return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
  1436.         }

  1438.         /* handle special cases */
  1439.         if (nans != FINITE || x.nans != FINITE) {
  1440.             if (isNaN()) {
  1441.                 return this;
  1442.             }

  1444.             if (x.isNaN()) {
  1445.                 return x;
  1446.             }

  1448.             if (nans == INFINITE && x.nans == FINITE) {
  1449.                 return this;
  1450.             }

  1452.             if (x.nans == INFINITE && nans == FINITE) {
  1453.                 return x;
  1454.             }

  1456.             if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
  1457.                 return x;
  1458.             }

  1460.             if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
  1461.                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1462.                 Dfp result = newInstance(getZero());
  1463.                 result.nans = QNAN;
  1464.                 result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
  1465.                 return result;
  1466.             }
  1467.         }

  1469.         /* copy this and the arg */
  1470.         Dfp a = newInstance(this);
  1471.         Dfp b = newInstance(x);

  1473.         /* initialize the result object */
  1474.         Dfp result = newInstance(getZero());

  1476.         /* Make all numbers positive, but remember their sign */
  1477.         final byte asign = a.sign;
  1478.         final byte bsign = b.sign;

  1480.         a.sign = 1;
  1481.         b.sign = 1;

  1483.         /* The result will be signed like the arg with greatest magnitude */
  1484.         byte rsign = bsign;
  1485.         if (compare(a, b) > 0) {
  1486.             rsign = asign;
  1487.         }

  1489.         /* Handle special case when a or b is zero, by setting the exponent
  1490.        of the zero number equal to the other one.  This avoids an alignment
  1491.        which would cause catastropic loss of precision */
  1492.         if (b.mant[mant.length-1] == 0) {
  1493.             b.exp = a.exp;
  1494.         }

  1496.         if (a.mant[mant.length-1] == 0) {
  1497.             a.exp = b.exp;
  1498.         }

  1500.         /* align number with the smaller exponent */
  1501.         int aextradigit = 0;
  1502.         int bextradigit = 0;
  1503.         if (a.exp < b.exp) {
  1504.             aextradigit = a.align(b.exp);
  1505.         } else {
  1506.             bextradigit = b.align(a.exp);
  1507.         }

  1509.         /* complement the smaller of the two if the signs are different */
  1510.         if (asign != bsign) {
  1511.             if (asign == rsign) {
  1512.                 bextradigit = b.complement(bextradigit);
  1513.             } else {
  1514.                 aextradigit = a.complement(aextradigit);
  1515.             }
  1516.         }

  1518.         /* add the mantissas */
  1519.         int rh = 0; /* acts as a carry */
  1520.         for (int i = 0; i < mant.length; i++) {
  1521.             final int r = a.mant[i]+b.mant[i]+rh;
  1522.             rh = r / RADIX;
  1523.             result.mant[i] = r - rh * RADIX;
  1524.         }
  1525.         result.exp = a.exp;
  1526.         result.sign = rsign;

  1528.         /* handle overflow -- note, when asign!=bsign an overflow is
  1529.          * normal and should be ignored.  */

  1531.         if (rh != 0 && (asign == bsign)) {
  1532.             final int lostdigit = result.mant[0];
  1533.             result.shiftRight();
  1534.             result.mant[mant.length-1] = rh;
  1535.             final int excp = result.round(lostdigit);
  1536.             if (excp != 0) {
  1537.                 result = dotrap(excp, ADD_TRAP, x, result);
  1538.             }
  1539.         }

  1541.         /* normalize the result */
  1542.         for (int i = 0; i < mant.length; i++) {
  1543.             if (result.mant[mant.length-1] != 0) {
  1544.                 break;
  1545.             }
  1546.             result.shiftLeft();
  1547.             if (i == 0) {
  1548.                 result.mant[0] = aextradigit+bextradigit;
  1549.                 aextradigit = 0;
  1550.                 bextradigit = 0;
  1551.             }
  1552.         }

  1554.         /* result is zero if after normalization the most sig. digit is zero */
  1555.         if (result.mant[mant.length-1] == 0) {
  1556.             result.exp = 0;

  1558.             if (asign != bsign) {
  1559.                 // Unless adding 2 negative zeros, sign is positive
  1560.                 result.sign = 1;  // Per IEEE 854-1987 Section 6.3
  1561.             }
  1562.         }

  1564.         /* Call round to test for over/under flows */
  1565.         final int excp = result.round(aextradigit + bextradigit);
  1566.         if (excp != 0) {
  1567.             result = dotrap(excp, ADD_TRAP, x, result);
  1568.         }

  1570.         return result;
  1571.     }

  1573.     /** Returns a number that is this number with the sign bit reversed.
  1574.      * @return the opposite of this
  1575.      */
  1576.     @Override
  1577.     public Dfp negate() {
  1578.         Dfp result = newInstance(this);
  1579.         result.sign = (byte) - result.sign;
  1580.         return result;
  1581.     }

  1583.     /** Subtract x from this.
  1584.      * @param x number to subtract
  1585.      * @return difference of this and a
  1586.      */
  1587.     @Override
  1588.     public Dfp subtract(final Dfp x) {
  1589.         return add(x.negate());
  1590.     }

  1592.     /** Round this given the next digit n using the current rounding mode.
  1593.      * @param n ???
  1594.      * @return the IEEE flag if an exception occurred
  1595.      */
  1596.     protected int round(int n) {
  1597.         boolean inc;
  1598.         switch (field.getRoundingMode()) {
  1599.             case ROUND_DOWN:
  1600.                 inc = false;
  1601.                 break;

  1603.             case ROUND_UP:
  1604.                 inc = n != 0;       // round up if n!=0
  1605.                 break;

  1607.             case ROUND_HALF_UP:
  1608.                 inc = n >= 5000;  // round half up
  1609.                 break;

  1611.             case ROUND_HALF_DOWN:
  1612.                 inc = n > 5000;  // round half down
  1613.                 break;

  1615.             case ROUND_HALF_EVEN:
  1616.                 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1);  // round half-even
  1617.                 break;

  1619.             case ROUND_HALF_ODD:
  1620.                 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0);  // round half-odd
  1621.                 break;

  1623.             case ROUND_CEIL:
  1624.                 inc = sign == 1 && n != 0;  // round ceil
  1625.                 break;

  1627.             case ROUND_FLOOR:
  1628.             default:
  1629.                 inc = sign == -1 && n != 0;  // round floor
  1630.                 break;
  1631.         }

  1633.         if (inc) {
  1634.             // increment if necessary
  1635.             int rh = 1;
  1636.             for (int i = 0; i < mant.length; i++) {
  1637.                 final int r = mant[i] + rh;
  1638.                 rh = r / RADIX;
  1639.                 mant[i] = r - rh * RADIX;
  1640.             }

  1642.             if (rh != 0) {
  1643.                 shiftRight();
  1644.                 mant[mant.length-1] = rh;
  1645.             }
  1646.         }

  1648.         // check for exceptional cases and raise signals if necessary
  1649.         if (exp < MIN_EXP) {
  1650.             // Gradual Underflow
  1651.             field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
  1652.             return DfpField.FLAG_UNDERFLOW;
  1653.         }

  1655.         if (exp > MAX_EXP) {
  1656.             // Overflow
  1657.             field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
  1658.             return DfpField.FLAG_OVERFLOW;
  1659.         }

  1661.         if (n != 0) {
  1662.             // Inexact
  1663.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
  1664.             return DfpField.FLAG_INEXACT;
  1665.         }

  1667.         return 0;

  1669.     }

  1671.     /** Multiply this by x.
  1672.      * @param x multiplicand
  1673.      * @return product of this and x
  1674.      */
  1675.     @Override
  1676.     public Dfp multiply(final Dfp x) {

  1678.         // make sure we don't mix number with different precision
  1679.         if (field.getRadixDigits() != x.field.getRadixDigits()) {
  1680.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1681.             final Dfp result = newInstance(getZero());
  1682.             result.nans = QNAN;
  1683.             return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
  1684.         }

  1686.         Dfp result = newInstance(getZero());

  1688.         /* handle special cases */
  1689.         if (nans != FINITE || x.nans != FINITE) {
  1690.             if (isNaN()) {
  1691.                 return this;
  1692.             }

  1694.             if (x.isNaN()) {
  1695.                 return x;
  1696.             }

  1698.             if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] != 0) {
  1699.                 result = newInstance(this);
  1700.                 result.sign = (byte) (sign * x.sign);
  1701.                 return result;
  1702.             }

  1704.             if (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] != 0) {
  1705.                 result = newInstance(x);
  1706.                 result.sign = (byte) (sign * x.sign);
  1707.                 return result;
  1708.             }

  1710.             if (x.nans == INFINITE && nans == INFINITE) {
  1711.                 result = newInstance(this);
  1712.                 result.sign = (byte) (sign * x.sign);
  1713.                 return result;
  1714.             }

  1716.             if ( (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] == 0) ||
  1717.                     (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] == 0) ) {
  1718.                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1719.                 result = newInstance(getZero());
  1720.                 result.nans = QNAN;
  1721.                 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
  1722.                 return result;
  1723.             }
  1724.         }

  1726.         int[] product = new int[mant.length*2];  // Big enough to hold even the largest result

  1728.         for (int i = 0; i < mant.length; i++) {
  1729.             int rh = 0;  // acts as a carry
  1730.             for (int j=0; j<mant.length; j++) {
  1731.                 int r = mant[i] * x.mant[j];    // multiply the 2 digits
  1732.                 r += product[i+j] + rh;  // add to the product digit with carry in

  1734.                 rh = r / RADIX;
  1735.                 product[i+j] = r - rh * RADIX;
  1736.             }
  1737.             product[i+mant.length] = rh;
  1738.         }

  1740.         // Find the most sig digit
  1741.         int md = mant.length * 2 - 1;  // default, in case result is zero
  1742.         for (int i = mant.length * 2 - 1; i >= 0; i--) {
  1743.             if (product[i] != 0) {
  1744.                 md = i;
  1745.                 break;
  1746.             }
  1747.         }

  1749.         // Copy the digits into the result
  1750.         for (int i = 0; i < mant.length; i++) {
  1751.             result.mant[mant.length - i - 1] = product[md - i];
  1752.         }

  1754.         // Fixup the exponent.
  1755.         result.exp = exp + x.exp + md - 2 * mant.length + 1;
  1756.         result.sign = (byte)((sign == x.sign)?1:-1);

  1758.         if (result.mant[mant.length-1] == 0) {
  1759.             // if result is zero, set exp to zero
  1760.             result.exp = 0;
  1761.         }

  1763.         final int excp;
  1764.         if (md > (mant.length-1)) {
  1765.             excp = result.round(product[md-mant.length]);
  1766.         } else {
  1767.             excp = result.round(0); // has no effect except to check status
  1768.         }

  1770.         if (excp != 0) {
  1771.             result = dotrap(excp, MULTIPLY_TRAP, x, result);
  1772.         }

  1774.         return result;

  1776.     }

  1778.     /** Multiply this by a single digit x.
  1779.      * @param x multiplicand
  1780.      * @return product of this and x
  1781.      */
  1782.     @Override
  1783.     public Dfp multiply(final int x) {
  1784.         if (x >= 0 && x < RADIX) {
  1785.             return multiplyFast(x);
  1786.         } else {
  1787.             return multiply(newInstance(x));
  1788.         }
  1789.     }

  1791.     /** Multiply this by a single digit 0&lt;=x&lt;radix.
  1792.      * There are speed advantages in this special case.
  1793.      * @param x multiplicand
  1794.      * @return product of this and x
  1795.      */
  1796.     private Dfp multiplyFast(final int x) {
  1797.         Dfp result = newInstance(this);

  1799.         /* handle special cases */
  1800.         if (nans != FINITE) {
  1801.             if (isNaN()) {
  1802.                 return this;
  1803.             }

  1805.             if (nans == INFINITE && x != 0) {
  1806.                 result = newInstance(this);
  1807.                 return result;
  1808.             }

  1810.             if (nans == INFINITE && x == 0) {
  1811.                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1812.                 result = newInstance(getZero());
  1813.                 result.nans = QNAN;
  1814.                 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
  1815.                 return result;
  1816.             }
  1817.         }

  1819.         /* range check x */
  1820.         if (x < 0 || x >= RADIX) {
  1821.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1822.             result = newInstance(getZero());
  1823.             result.nans = QNAN;
  1824.             result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
  1825.             return result;
  1826.         }

  1828.         int rh = 0;
  1829.         for (int i = 0; i < mant.length; i++) {
  1830.             final int r = mant[i] * x + rh;
  1831.             rh = r / RADIX;
  1832.             result.mant[i] = r - rh * RADIX;
  1833.         }

  1835.         int lostdigit = 0;
  1836.         if (rh != 0) {
  1837.             lostdigit = result.mant[0];
  1838.             result.shiftRight();
  1839.             result.mant[mant.length-1] = rh;
  1840.         }

  1842.         if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
  1843.             result.exp = 0;
  1844.         }

  1846.         final int excp = result.round(lostdigit);
  1847.         if (excp != 0) {
  1848.             result = dotrap(excp, MULTIPLY_TRAP, result, result);
  1849.         }

  1851.         return result;
  1852.     }

  1854.     /** {@inheritDoc} */
  1855.     @Override
  1856.     public Dfp square() {
  1857.         return multiply(this);
  1858.     }

  1860.     /** Divide this by divisor.
  1861.      * @param divisor divisor
  1862.      * @return quotient of this by divisor
  1863.      */
  1864.     @Override
  1865.     public Dfp divide(Dfp divisor) {
  1866.         int[] dividend; // current status of the dividend
  1867.         int[] quotient; // quotient
  1868.         int[] remainder;// remainder
  1869.         int qd;         // current quotient digit we're working with
  1870.         int nsqd;       // number of significant quotient digits we have
  1871.         int trial=0;    // trial quotient digit
  1872.         int minadj;     // minimum adjustment
  1873.         boolean trialgood; // Flag to indicate a good trail digit
  1874.         int md;         // most sig digit in result
  1875.         int excp;       // exceptions

  1877.         // make sure we don't mix number with different precision
  1878.         if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
  1879.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1880.             final Dfp result = newInstance(getZero());
  1881.             result.nans = QNAN;
  1882.             return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
  1883.         }

  1885.         Dfp result = newInstance(getZero());

  1887.         /* handle special cases */
  1888.         if (nans != FINITE || divisor.nans != FINITE) {
  1889.             if (isNaN()) {
  1890.                 return this;
  1891.             }

  1893.             if (divisor.isNaN()) {
  1894.                 return divisor;
  1895.             }

  1897.             if (nans == INFINITE && divisor.nans == FINITE) {
  1898.                 result = newInstance(this);
  1899.                 result.sign = (byte) (sign * divisor.sign);
  1900.                 return result;
  1901.             }

  1903.             if (divisor.nans == INFINITE && nans == FINITE) {
  1904.                 result = newInstance(getZero());
  1905.                 result.sign = (byte) (sign * divisor.sign);
  1906.                 return result;
  1907.             }

  1909.             if (divisor.nans == INFINITE && nans == INFINITE) {
  1910.                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  1911.                 result = newInstance(getZero());
  1912.                 result.nans = QNAN;
  1913.                 result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
  1914.                 return result;
  1915.             }
  1916.         }

  1918.         /* Test for divide by zero */
  1919.         if (divisor.mant[mant.length-1] == 0) {
  1920.             field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
  1921.             result = newInstance(getZero());
  1922.             result.sign = (byte) (sign * divisor.sign);
  1923.             result.nans = INFINITE;
  1924.             result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
  1925.             return result;
  1926.         }

  1928.         dividend = new int[mant.length+1];  // one extra digit needed
  1929.         quotient = new int[mant.length+2];  // two extra digits needed 1 for overflow, 1 for rounding
  1930.         remainder = new int[mant.length+1]; // one extra digit needed

  1932.         /* Initialize our most significant digits to zero */

  1934.         dividend[mant.length] = 0;
  1935.         quotient[mant.length] = 0;
  1936.         quotient[mant.length+1] = 0;
  1937.         remainder[mant.length] = 0;

  1939.         /* copy our mantissa into the dividend, initialize the
  1940.        quotient while we are at it */

  1942.         for (int i = 0; i < mant.length; i++) {
  1943.             dividend[i] = mant[i];
  1944.             quotient[i] = 0;
  1945.             remainder[i] = 0;
  1946.         }

  1948.         /* outer loop.  Once per quotient digit */
  1949.         nsqd = 0;
  1950.         for (qd = mant.length+1; qd >= 0; qd--) {
  1951.             /* Determine outer limits of our quotient digit */

  1953.             // r =  most sig 2 digits of dividend
  1954.             final int divMsb = dividend[mant.length]*RADIX+dividend[mant.length-1];
  1955.             int min = divMsb       / (divisor.mant[mant.length-1]+1);
  1956.             int max = (divMsb + 1) / divisor.mant[mant.length-1];

  1958.             trialgood = false;
  1959.             while (!trialgood) {
  1960.                 // try the mean
  1961.                 trial = (min+max)/2;

  1963.                 /* Multiply by divisor and store as remainder */
  1964.                 int rh = 0;
  1965.                 for (int i = 0; i < mant.length + 1; i++) {
  1966.                     int dm = (i<mant.length)?divisor.mant[i]:0;
  1967.                     final int r = (dm * trial) + rh;
  1968.                     rh = r / RADIX;
  1969.                     remainder[i] = r - rh * RADIX;
  1970.                 }

  1972.                 /* subtract the remainder from the dividend */
  1973.                 rh = 1;  // carry in to aid the subtraction
  1974.                 for (int i = 0; i < mant.length + 1; i++) {
  1975.                     final int r = ((RADIX-1) - remainder[i]) + dividend[i] + rh;
  1976.                     rh = r / RADIX;
  1977.                     remainder[i] = r - rh * RADIX;
  1978.                 }

  1980.                 /* Lets analyze what we have here */
  1981.                 if (rh == 0) {
  1982.                     // trial is too big -- negative remainder
  1983.                     max = trial-1;
  1984.                     continue;
  1985.                 }

  1987.                 /* find out how far off the remainder is telling us we are */
  1988.                 minadj = (remainder[mant.length] * RADIX)+remainder[mant.length-1];
  1989.                 minadj /= divisor.mant[mant.length-1] + 1;

  1991.                 if (minadj >= 2) {
  1992.                     min = trial+minadj;  // update the minimum
  1993.                     continue;
  1994.                 }

  1996.                 /* May have a good one here, check more thoroughly.  Basically
  1997.            its a good one if it is less than the divisor */
  1998.                 trialgood = false;  // assume false
  1999.                 for (int i = mant.length - 1; i >= 0; i--) {
  2000.                     if (divisor.mant[i] > remainder[i]) {
  2001.                         trialgood = true;
  2002.                     }
  2003.                     if (divisor.mant[i] < remainder[i]) {
  2004.                         break;
  2005.                     }
  2006.                 }

  2008.                 if (remainder[mant.length] != 0) {
  2009.                     trialgood = false;
  2010.                 }

  2012.                 if (!trialgood) {
  2013.                     min = trial+1;
  2014.                 }
  2015.             }

  2017.             /* Great we have a digit! */
  2018.             quotient[qd] = trial;
  2019.             if (trial != 0 || nsqd != 0) {
  2020.                 nsqd++;
  2021.             }

  2023.             if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
  2024.                 // We have enough for this mode
  2025.                 break;
  2026.             }

  2028.             if (nsqd > mant.length) {
  2029.                 // We have enough digits
  2030.                 break;
  2031.             }

  2033.             /* move the remainder into the dividend while left shifting */
  2034.             dividend[0] = 0;
  2035.             System.arraycopy(remainder, 0, dividend, 1, mant.length);
  2036.         }

  2038.         /* Find the most sig digit */
  2039.         md = mant.length;  // default
  2040.         for (int i = mant.length + 1; i >= 0; i--) {
  2041.             if (quotient[i] != 0) {
  2042.                 md = i;
  2043.                 break;
  2044.             }
  2045.         }

  2047.         /* Copy the digits into the result */
  2048.         for (int i=0; i<mant.length; i++) {
  2049.             result.mant[mant.length-i-1] = quotient[md-i];
  2050.         }

  2052.         /* Fixup the exponent. */
  2053.         result.exp = exp - divisor.exp + md - mant.length;
  2054.         result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);

  2056.         if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
  2057.             result.exp = 0;
  2058.         }

  2060.         if (md > (mant.length-1)) {
  2061.             excp = result.round(quotient[md-mant.length]);
  2062.         } else {
  2063.             excp = result.round(0);
  2064.         }

  2066.         if (excp != 0) {
  2067.             result = dotrap(excp, DIVIDE_TRAP, divisor, result);
  2068.         }

  2070.         return result;
  2071.     }

  2073.     /** Divide by a single digit less than radix.
  2074.      *  Special case, so there are speed advantages. 0 &lt;= divisor &lt; radix
  2075.      * @param divisor divisor
  2076.      * @return quotient of this by divisor
  2077.      */
  2078.     public Dfp divide(int divisor) {

  2080.         // Handle special cases
  2081.         if (nans != FINITE) {
  2082.             if (isNaN()) {
  2083.                 return this;
  2084.             }

  2086.             if (nans == INFINITE) {
  2087.                 return newInstance(this);
  2088.             }
  2089.         }

  2091.         // Test for divide by zero
  2092.         if (divisor == 0) {
  2093.             field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
  2094.             Dfp result = newInstance(getZero());
  2095.             result.sign = sign;
  2096.             result.nans = INFINITE;
  2097.             result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
  2098.             return result;
  2099.         }

  2101.         // range check divisor
  2102.         if (divisor < 0 || divisor >= RADIX) {
  2103.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  2104.             Dfp result = newInstance(getZero());
  2105.             result.nans = QNAN;
  2106.             result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
  2107.             return result;
  2108.         }

  2110.         Dfp result = newInstance(this);

  2112.         int rl = 0;
  2113.         for (int i = mant.length-1; i >= 0; i--) {
  2114.             final int r = rl*RADIX + result.mant[i];
  2115.             final int rh = r / divisor;
  2116.             rl = r - rh * divisor;
  2117.             result.mant[i] = rh;
  2118.         }

  2120.         if (result.mant[mant.length-1] == 0) {
  2121.             // normalize
  2122.             result.shiftLeft();
  2123.             final int r = rl * RADIX;        // compute the next digit and put it in
  2124.             final int rh = r / divisor;
  2125.             rl = r - rh * divisor;
  2126.             result.mant[0] = rh;
  2127.         }

  2129.         final int excp = result.round(rl * RADIX / divisor);  // do the rounding
  2130.         if (excp != 0) {
  2131.             result = dotrap(excp, DIVIDE_TRAP, result, result);
  2132.         }

  2134.         return result;

  2136.     }

  2138.     /** {@inheritDoc} */
  2139.     @Override
  2140.     public Dfp reciprocal() {
  2141.         return field.getOne().divide(this);
  2142.     }

  2144.     /** Compute the square root.
  2145.      * @return square root of the instance
  2146.      */
  2147.     @Override
  2148.     public Dfp sqrt() {

  2150.         // check for unusual cases
  2151.         if (nans == FINITE && mant[mant.length-1] == 0) {
  2152.             // if zero
  2153.             return newInstance(this);
  2154.         }

  2156.         if (nans != FINITE) {
  2157.             if (nans == INFINITE && sign == 1) {
  2158.                 // if positive infinity
  2159.                 return newInstance(this);
  2160.             }

  2162.             if (nans == QNAN) {
  2163.                 return newInstance(this);
  2164.             }

  2166.             if (nans == SNAN) {
  2167.                 Dfp result;

  2169.                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  2170.                 result = newInstance(this);
  2171.                 result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
  2172.                 return result;
  2173.             }
  2174.         }

  2176.         if (sign == -1) {
  2177.             // if negative
  2178.             Dfp result;

  2180.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  2181.             result = newInstance(this);
  2182.             result.nans = QNAN;
  2183.             result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
  2184.             return result;
  2185.         }

  2187.         Dfp x = newInstance(this);

  2189.         /* Lets make a reasonable guess as to the size of the square root */
  2190.         if (x.exp < -1 || x.exp > 1) {
  2191.             x.exp = this.exp / 2;
  2192.         }

  2194.         /* Coarsely estimate the mantissa */
  2195.         switch (x.mant[mant.length-1] / 2000) {
  2196.             case 0:
  2197.                 x.mant[mant.length-1] = x.mant[mant.length-1]/2+1;
  2198.                 break;
  2199.             case 2:
  2200.                 x.mant[mant.length-1] = 1500;
  2201.                 break;
  2202.             case 3:
  2203.                 x.mant[mant.length-1] = 2200;
  2204.                 break;
  2205.             default:
  2206.                 x.mant[mant.length-1] = 3000;
  2207.                 break;
  2208.         }

  2210.         /* Now that we have the first pass estimate, compute the rest
  2211.        by the formula dx = (y - x*x) / (2x); */

  2213.         Dfp dx;
  2214.         Dfp px  = getZero();
  2215.         Dfp ppx;
  2216.         while (x.unequal(px)) {
  2217.             dx = newInstance(x);
  2218.             dx.sign = -1;
  2219.             dx = dx.add(this.divide(x));
  2220.             dx = dx.divide(2);
  2221.             ppx = px;
  2222.             px = x;
  2223.             x = x.add(dx);

  2225.             if (x.equals(ppx)) {
  2226.                 // alternating between two values
  2227.                 break;
  2228.             }

  2230.             // if dx is zero, break.  Note testing the most sig digit
  2231.             // is a sufficient test since dx is normalized
  2232.             if (dx.mant[mant.length-1] == 0) {
  2233.                 break;
  2234.             }
  2235.         }

  2237.         return x;

  2239.     }

  2241.     /** Get a string representation of the instance.
  2242.      * @return string representation of the instance
  2243.      */
  2244.     @Override
  2245.     public String toString() {
  2246.         if (nans != FINITE) {
  2247.             // if non-finite exceptional cases
  2248.             if (nans == INFINITE) {
  2249.                 return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
  2250.             } else {
  2251.                 return NAN_STRING;
  2252.             }
  2253.         }

  2255.         if (exp > mant.length || exp < -1) {
  2256.             return dfp2sci();
  2257.         }

  2259.         return dfp2string();

  2261.     }

  2263.     /** Convert an instance to a string using scientific notation.
  2264.      * @return string representation of the instance in scientific notation
  2265.      */
  2266.     protected String dfp2sci() {
  2267.         char[] rawdigits = new char[mant.length * 4];
  2268.         int p;
  2269.         int e;
  2270.         int ae;
  2271.         int shf;

  2273.         // Get all the digits
  2274.         p = 0;
  2275.         for (int i = mant.length - 1; i >= 0; i--) {
  2276.             rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
  2277.             rawdigits[p++] = (char) (((mant[i] / 100) %10) + '0');
  2278.             rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
  2279.             rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
  2280.         }

  2282.         // Find the first non-zero one
  2283.         for (p = 0; p < rawdigits.length; p++) {
  2284.             if (rawdigits[p] != '0') {
  2285.                 break;
  2286.             }
  2287.         }
  2288.         shf = p;

  2290.         // Now do the conversion
  2291.         StringBuilder builder = new StringBuilder();
  2292.         if (sign == -1) {
  2293.             builder.append('-');
  2294.         }

  2296.         if (p != rawdigits.length) {
  2297.             // there are non zero digits...
  2298.             builder.append(rawdigits[p++]).append('.');

  2300.             while (p<rawdigits.length) {
  2301.                 builder.append(rawdigits[p++]);
  2302.             }
  2303.         } else {
  2304.             builder.append("0.0e0");
  2305.             return builder.toString();
  2306.         }

  2308.         builder.append('e');

  2310.         // Find the msd of the exponent

  2312.         e = exp * 4 - shf - 1;
  2313.         ae = e;
  2314.         if (e < 0) {
  2315.             ae = -e;
  2316.         }

  2318.         // Find the largest p such that p < e
  2319.         for (p = 1000000000; p > ae; p /= 10) { // NOPMD - empty loop is normal here
  2320.             // nothing to do
  2321.         }

  2323.         if (e < 0) {
  2324.             builder.append('-');
  2325.         }

  2327.         while (p > 0) {
  2328.             builder.append((char)(ae / p + '0'));
  2329.             ae %= p;
  2330.             p /= 10;
  2331.         }

  2333.         return builder.toString();

  2335.     }

  2337.     /** Convert an instance to a string using normal notation.
  2338.      * @return string representation of the instance in normal notation
  2339.      */
  2340.     protected String dfp2string() {
  2341.         final String fourZero = "0000";
  2342.         int e = exp;
  2343.         boolean pointInserted = false;

  2345.         StringBuilder builder = new StringBuilder(23);

  2347.         if (e <= 0) {
  2348.             builder.append("0.");
  2349.             pointInserted = true;
  2350.         }

  2352.         while (e < 0) {
  2353.             builder.append(fourZero);
  2354.             e++;
  2355.         }

  2357.         for (int i = mant.length - 1; i >= 0; i--) {
  2358.             builder.
  2359.                 append((char) ((mant[i] / 1000) + '0')).
  2360.                 append((char) (((mant[i] / 100) % 10) + '0')).
  2361.                 append((char) (((mant[i] / 10) % 10) + '0')).
  2362.                 append((char) (((mant[i]) % 10) + '0'));
  2363.             --e;
  2364.             if (e == 0) {
  2365.                 builder.append('.');
  2366.                 pointInserted = true;
  2367.             }
  2368.         }

  2370.         while (e > 0) {
  2371.             builder.append(fourZero);
  2372.             e--;
  2373.         }

  2375.         if (!pointInserted) {
  2376.             // Ensure we have a radix point!
  2377.             builder.append('.');
  2378.         }

  2380.         // Suppress leading zeros
  2381.         while (builder.charAt(0) == '0') {
  2382.             builder.deleteCharAt(0);
  2383.         }
  2384.         if (builder.charAt(0) == '.') {
  2385.             builder.insert(0, '0');
  2386.         }

  2388.         // Suppress trailing zeros
  2389.         while (builder.charAt(builder.length() - 1) == '0') {
  2390.             builder.deleteCharAt(builder.length() - 1);
  2391.         }

  2393.         // Insert sign
  2394.         if (sign < 0) {
  2395.             builder.insert(0, '-');
  2396.         }

  2398.         return builder.toString();

  2400.     }

  2402.     /** Raises a trap.  This does not set the corresponding flag however.
  2403.      *  @param type the trap type
  2404.      *  @param what - name of routine trap occurred in
  2405.      *  @param oper - input operator to function
  2406.      *  @param result - the result computed prior to the trap
  2407.      *  @return The suggested return value from the trap handler
  2408.      */
  2409.     public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
  2410.         Dfp def = result;

  2412.         switch (type) {
  2413.             case DfpField.FLAG_INVALID:
  2414.                 def = newInstance(getZero());
  2415.                 def.sign = result.sign;
  2416.                 def.nans = QNAN;
  2417.                 break;

  2419.             case DfpField.FLAG_DIV_ZERO:
  2420.                 if (nans == FINITE && mant[mant.length-1] != 0) {
  2421.                     // normal case, we are finite, non-zero
  2422.                     def = newInstance(getZero());
  2423.                     def.sign = (byte)(sign*oper.sign);
  2424.                     def.nans = INFINITE;
  2425.                 }

  2427.                 if (nans == FINITE && mant[mant.length-1] == 0) {
  2428.                     //  0/0
  2429.                     def = newInstance(getZero());
  2430.                     def.nans = QNAN;
  2431.                 }

  2433.                 if (nans == INFINITE || nans == QNAN) {
  2434.                     def = newInstance(getZero());
  2435.                     def.nans = QNAN;
  2436.                 }

  2438.                 if (nans == INFINITE || nans == SNAN) {
  2439.                     def = newInstance(getZero());
  2440.                     def.nans = QNAN;
  2441.                 }
  2442.                 break;

  2444.             case DfpField.FLAG_UNDERFLOW:
  2445.                 if ( (result.exp+mant.length) < MIN_EXP) {
  2446.                     def = newInstance(getZero());
  2447.                     def.sign = result.sign;
  2448.                 } else {
  2449.                     def = newInstance(result);  // gradual underflow
  2450.                 }
  2451.                 result.exp += ERR_SCALE;
  2452.                 break;

  2454.             case DfpField.FLAG_OVERFLOW:
  2455.                 result.exp -= ERR_SCALE;
  2456.                 def = newInstance(getZero());
  2457.                 def.sign = result.sign;
  2458.                 def.nans = INFINITE;
  2459.                 break;

  2461.             default: def = result; break;
  2462.         }

  2464.         return trap(type, what, oper, def, result);

  2466.     }

  2468.     /** Trap handler.  Subclasses may override this to provide trap
  2469.      *  functionality per IEEE 854-1987.
  2470.      *
  2471.      *  @param type  The exception type - e.g. FLAG_OVERFLOW
  2472.      *  @param what  The name of the routine we were in e.g. divide()
  2473.      *  @param oper  An operand to this function if any
  2474.      *  @param def   The default return value if trap not enabled
  2475.      *  @param result    The result that is specified to be delivered per
  2476.      *                   IEEE 854, if any
  2477.      *  @return the value that should be return by the operation triggering the trap
  2478.      */
  2479.     protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
  2480.         return def;
  2481.     }

  2483.     /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN.
  2484.      * @return type of the number
  2485.      */
  2486.     public int classify() {
  2487.         return nans;
  2488.     }

  2490.     /** Creates an instance that is the same as x except that it has the sign of y.
  2491.      * abs(x) = dfp.copysign(x, dfp.one)
  2492.      * @param x number to get the value from
  2493.      * @param y number to get the sign from
  2494.      * @return a number with the value of x and the sign of y
  2495.      */
  2496.     public static Dfp copysign(final Dfp x, final Dfp y) {
  2497.         Dfp result = x.newInstance(x);
  2498.         result.sign = y.sign;
  2499.         return result;
  2500.     }

  2502.     /** Returns the next number greater than this one in the direction of x.
  2503.      * If this==x then simply returns this.
  2504.      * @param x direction where to look at
  2505.      * @return closest number next to instance in the direction of x
  2506.      */
  2507.     public Dfp nextAfter(final Dfp x) {

  2509.         // make sure we don't mix number with different precision
  2510.         if (field.getRadixDigits() != x.field.getRadixDigits()) {
  2511.             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
  2512.             final Dfp result = newInstance(getZero());
  2513.             result.nans = QNAN;
  2514.             return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
  2515.         }

  2517.         // if this is greater than x
  2518.         boolean up = false;
  2519.         if (this.lessThan(x)) {
  2520.             up = true;
  2521.         }

  2523.         if (compare(this, x) == 0) {
  2524.             return newInstance(x);
  2525.         }

  2527.         if (lessThan(getZero())) {
  2528.             up = !up;
  2529.         }

  2531.         final Dfp inc;
  2532.         Dfp result;
  2533.         if (up) {
  2534.             inc = newInstance(getOne());
  2535.             inc.exp = this.exp-mant.length+1;
  2536.             inc.sign = this.sign;

  2538.             if (this.equals(getZero())) {
  2539.                 inc.exp = MIN_EXP-mant.length;
  2540.             }

  2542.             result = add(inc);
  2543.         } else {
  2544.             inc = newInstance(getOne());
  2545.             inc.exp = this.exp;
  2546.             inc.sign = this.sign;

  2548.             if (this.equals(inc)) {
  2549.                 inc.exp = this.exp-mant.length;
  2550.             } else {
  2551.                 inc.exp = this.exp-mant.length+1;
  2552.             }

  2554.             if (this.equals(getZero())) {
  2555.                 inc.exp = MIN_EXP-mant.length;
  2556.             }

  2558.             result = this.subtract(inc);
  2559.         }

  2561.         if (result.classify() == INFINITE && this.classify() != INFINITE) {
  2562.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
  2563.             result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
  2564.         }

  2566.         if (result.equals(getZero()) && !this.equals(getZero())) {
  2567.             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
  2568.             result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
  2569.         }

  2571.         return result;

  2573.     }

  2575.     /** Convert the instance into a double.
  2576.      * @return a double approximating the instance
  2577.      * @see #toSplitDouble()
  2578.      */
  2579.     public double toDouble() {

  2581.         if (isInfinite()) {
  2582.             if (lessThan(getZero())) {
  2583.                 return Double.NEGATIVE_INFINITY;
  2584.             } else {
  2585.                 return Double.POSITIVE_INFINITY;
  2586.             }
  2587.         }

  2589.         if (isNaN()) {
  2590.             return Double.NaN;
  2591.         }

  2593.         Dfp y = this;
  2594.         boolean negate = false;
  2595.         int cmp0 = compare(this, getZero());
  2596.         if (cmp0 == 0) {
  2597.             return sign < 0 ? -0.0 : +0.0;
  2598.         } else if (cmp0 < 0) {
  2599.             y = negate();
  2600.             negate = true;
  2601.         }

  2603.         /* Find the exponent, first estimate by integer log10, then adjust.
  2604.          Should be faster than doing a natural logarithm.  */
  2605.         int exponent = (int)(y.intLog10() * 3.32);
  2606.         if (exponent < 0) {
  2607.             exponent--;
  2608.         }

  2610.         Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
  2611.         while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
  2612.             tempDfp = tempDfp.multiply(2);
  2613.             exponent++;
  2614.         }
  2615.         exponent--;

  2617.         /* We have the exponent, now work on the mantissa */

  2619.         y = y.divide(DfpMath.pow(getTwo(), exponent));
  2620.         if (exponent > -1023) {
  2621.             y = y.subtract(getOne());
  2622.         }

  2624.         if (exponent < -1074) {
  2625.             return 0;
  2626.         }

  2628.         if (exponent > 1023) {
  2629.             return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
  2630.         }


  2633.         y = y.multiply(newInstance(4503599627370496L)).rint();
  2634.         String str = y.toString();
  2635.         str = str.substring(0, str.length()-1);
  2636.         long mantissa = Long.parseLong(str);

  2638.         if (mantissa == 4503599627370496L) {
  2639.             // Handle special case where we round up to next power of two
  2640.             mantissa = 0;
  2641.             exponent++;
  2642.         }

  2644.         /* Its going to be subnormal, so make adjustments */
  2645.         if (exponent <= -1023) {
  2646.             exponent--;
  2647.         }

  2649.         while (exponent < -1023) {
  2650.             exponent++;
  2651.             mantissa >>>= 1;
  2652.         }

  2654.         long bits = mantissa | ((exponent + 1023L) << 52);
  2655.         double x = Double.longBitsToDouble(bits);

  2657.         if (negate) {
  2658.             x = -x;
  2659.         }

  2661.         return x;

  2663.     }

  2665.     /** Convert the instance into a split double.
  2666.      * @return an array of two doubles which sum represent the instance
  2667.      * @see #toDouble()
  2668.      */
  2669.     public double[] toSplitDouble() {
  2670.         double[] split = new double[2];
  2671.         long mask = 0xffffffffc0000000L;

  2673.         split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
  2674.         split[1] = subtract(newInstance(split[0])).toDouble();

  2676.         return split;
  2677.     }

  2679.     /** {@inheritDoc}
  2680.      */
  2681.     @Override
  2682.     public double getReal() {
  2683.         return toDouble();
  2684.     }

  2686.     /** {@inheritDoc} */
  2687.     @Override
  2688.     public Dfp getAddendum() {
  2689.         return isFinite() ? subtract(getReal()) : getZero();
  2690.     }

  2692.     /** {@inheritDoc}
  2693.      */
  2694.     @Override
  2695.     public Dfp remainder(final double a) {
  2696.         return remainder(newInstance(a));
  2697.     }

  2699.     /** {@inheritDoc}
  2700.      */
  2701.     @Override
  2702.     public Dfp sign() {
  2703.         if (isNaN() || isZero()) {
  2704.             return this;
  2705.         } else {
  2706.             return newInstance(sign > 0 ? +1 : -1);
  2707.         }
  2708.     }

  2710.     /** {@inheritDoc}
  2711.      */
  2712.     @Override
  2713.     public Dfp copySign(final Dfp s) {
  2714.         if ((sign >= 0 && s.sign >= 0) || (sign < 0 && s.sign < 0)) { // Sign is currently OK
  2715.             return this;
  2716.         }
  2717.         return negate(); // flip sign
  2718.     }

  2720.     /** {@inheritDoc}
  2721.      */
  2722.     @Override
  2723.     public Dfp copySign(final double s) {
  2724.         long sb = Double.doubleToLongBits(s);
  2725.         if ((sign >= 0 && sb >= 0) || (sign < 0 && sb < 0)) { // Sign is currently OK
  2726.             return this;
  2727.         }
  2728.         return negate(); // flip sign
  2729.     }

  2731.     /** {@inheritDoc}
  2732.      */
  2733.     @Override
  2734.     public int getExponent() {

  2736.         if (nans != FINITE) {
  2737.             // 2⁴³⁵⁴¹¹ < 10000³²⁷⁶⁸ < 2⁴³⁵⁴¹²
  2738.             return 435411;
  2739.         }
  2740.         if (isZero()) {
  2741.             return -435412;
  2742.         }

  2744.         final Dfp abs = abs();

  2746.         // estimate a lower bound for binary exponent
  2747.         // 13301/1001 is a continued fraction approximation of ln(10000)/ln(2)
  2748.         int p = FastMath.max(13301 * exp / 1001 - 15, -435411);
  2749.         Dfp twoP = DfpMath.pow(getTwo(), p);
  2750.         while (compare(abs, twoP) >= 0) {
  2751.             twoP = twoP.add(twoP);
  2752.             ++p;
  2753.         }

  2755.         return p - 1;

  2757.     }

  2759.     /** {@inheritDoc}
  2760.      */
  2761.     @Override
  2762.     public Dfp scalb(final int n) {
  2763.         return multiply(DfpMath.pow(getTwo(), n));
  2764.     }

  2766.     /** {@inheritDoc}
  2767.      */
  2768.     @Override
  2769.     public Dfp ulp() {
  2770.         final Dfp result = new Dfp(field);
  2771.         result.mant[result.mant.length - 1] = 1;
  2772.         result.exp = exp - (result.mant.length - 1);
  2773.         return result;
  2774.     }

  2776.     /** {@inheritDoc}
  2777.      */
  2778.     @Override
  2779.     public Dfp hypot(final Dfp y) {

  2781.         if (isInfinite() || y.isInfinite()) {
  2782.             return field.newDfp(Double.POSITIVE_INFINITY);
  2783.         } else if (isNaN() || y.isNaN()) {
  2784.             return field.newDfp(Double.NaN);
  2785.         } else {
  2786.             // find scaling to avoid both overflow and underflow
  2787.             final int scalingExp = (exp + y.exp) / 2;

  2789.             // scale both operands
  2790.             final Dfp scaledX = new Dfp(this);
  2791.             scaledX.exp -= scalingExp;
  2792.             final Dfp scaledY = new Dfp(y);
  2793.             scaledY.exp -= scalingExp;

  2795.             // compute scaled hypothenuse
  2796.             final Dfp h = scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();

  2798.             // scale result
  2799.             h.exp += scalingExp;

  2801.             return h;
  2802.         }

  2804.     }

  2806.     /** {@inheritDoc}
  2807.      */
  2808.     @Override
  2809.     public Dfp rootN(final int n) {
  2810.         return (sign >= 0) ?
  2811.                DfpMath.pow(this, getOne().divide(n)) :
  2812.                DfpMath.pow(negate(), getOne().divide(n)).negate();
  2813.     }

  2815.     /** {@inheritDoc}
  2816.      */
  2817.     @Override
  2818.     public Dfp pow(final double p) {
  2819.         return DfpMath.pow(this, newInstance(p));
  2820.     }

  2822.     /** {@inheritDoc}
  2823.      */
  2824.     @Override
  2825.     public Dfp pow(final int n) {
  2826.         return DfpMath.pow(this, n);
  2827.     }

  2829.     /** {@inheritDoc}
  2830.      */
  2831.     @Override
  2832.     public Dfp pow(final Dfp e) {
  2833.         return DfpMath.pow(this, e);
  2834.     }

  2836.     /** {@inheritDoc}
  2837.      */
  2838.     @Override
  2839.     public Dfp exp() {
  2840.         return DfpMath.exp(this);
  2841.     }

  2843.     /** {@inheritDoc}
  2844.      */
  2845.     @Override
  2846.     public Dfp expm1() {
  2847.         return DfpMath.exp(this).subtract(getOne());
  2848.     }

  2850.     /** {@inheritDoc}
  2851.      */
  2852.     @Override
  2853.     public Dfp log() {
  2854.         return DfpMath.log(this);
  2855.     }

  2857.     /** {@inheritDoc}
  2858.      */
  2859.     @Override
  2860.     public Dfp log1p() {
  2861.         return DfpMath.log(this.add(getOne()));
  2862.     }

  2864.     /** {@inheritDoc}
  2865.      */
  2866.     @Override
  2867.     public Dfp log10() {
  2868.         return DfpMath.log(this).divide(DfpMath.log(newInstance(10)));
  2869.     }

  2871.     /** {@inheritDoc}
  2872.      */
  2873.     @Override
  2874.     public Dfp cos() {
  2875.         return DfpMath.cos(this);
  2876.     }

  2878.     /** {@inheritDoc}
  2879.      */
  2880.     @Override
  2881.     public Dfp sin() {
  2882.         return DfpMath.sin(this);
  2883.     }

  2885.     /** {@inheritDoc}
  2886.      */
  2887.     @Override
  2888.     public Dfp tan() {
  2889.         return DfpMath.tan(this);
  2890.     }

  2892.     /** {@inheritDoc}
  2893.      */
  2894.     @Override
  2895.     public Dfp acos() {
  2896.         return DfpMath.acos(this);
  2897.     }

  2899.     /** {@inheritDoc}
  2900.      */
  2901.     @Override
  2902.     public Dfp asin() {
  2903.         return DfpMath.asin(this);
  2904.     }

  2906.     /** {@inheritDoc}
  2907.      */
  2908.     @Override
  2909.     public Dfp atan() {
  2910.         return DfpMath.atan(this);
  2911.     }

  2913.     /** {@inheritDoc}
  2914.      */
  2915.     @Override
  2916.     public Dfp atan2(final Dfp x)
  2917.         throws MathIllegalArgumentException {

  2919.         // compute r = sqrt(x^2+y^2)
  2920.         final Dfp r = x.square().add(multiply(this)).sqrt();
  2921.         if (r.isZero()) {
  2922.             // special cases handling
  2923.             if (x.sign >= 0) {
  2924.                 return this; // ±0.0
  2925.             } else {
  2926.                 return newInstance((sign <= 0) ? -FastMath.PI : FastMath.PI); // ±π
  2927.             }
  2928.         }

  2930.         if (x.sign >= 0) {

  2932.             // compute atan2(y, x) = 2 atan(y / (r + x))
  2933.             return getTwo().multiply(divide(r.add(x)).atan());

  2935.         } else {

  2937.             // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
  2938.             final Dfp tmp = getTwo().multiply(divide(r.subtract(x)).atan());
  2939.             final Dfp pmPi = newInstance((tmp.sign <= 0) ? -FastMath.PI : FastMath.PI);
  2940.             return pmPi.subtract(tmp);

  2942.         }

  2944.     }

  2946.     /** {@inheritDoc}
  2947.      */
  2948.     @Override
  2949.     public Dfp cosh() {
  2950.         return DfpMath.exp(this).add(DfpMath.exp(negate())).multiply(0.5);
  2951.     }

  2953.     /** {@inheritDoc}
  2954.      */
  2955.     @Override
  2956.     public Dfp sinh() {
  2957.         return DfpMath.exp(this).subtract(DfpMath.exp(negate())).multiply(0.5);
  2958.     }

  2960.     /** {@inheritDoc}
  2961.      */
  2962.     @Override
  2963.     public FieldSinhCosh<Dfp> sinhCosh() {
  2964.         final Dfp p = DfpMath.exp(this);
  2965.         final Dfp m = DfpMath.exp(negate());
  2966.         return new FieldSinhCosh<>(p.subtract(m).multiply(0.5), p.add(m).multiply(0.5));
  2967.     }

  2969.     /** {@inheritDoc}
  2970.      */
  2971.     @Override
  2972.     public Dfp tanh() {
  2973.         final Dfp ePlus  = DfpMath.exp(this);
  2974.         final Dfp eMinus = DfpMath.exp(negate());
  2975.         return ePlus.subtract(eMinus).divide(ePlus.add(eMinus));
  2976.     }

  2978.     /** {@inheritDoc}
  2979.      */
  2980.     @Override
  2981.     public Dfp acosh() {
  2982.         return square().subtract(getOne()).sqrt().add(this).log();
  2983.     }

  2985.     /** {@inheritDoc}
  2986.      */
  2987.     @Override
  2988.     public Dfp asinh() {
  2989.         return square().add(getOne()).sqrt().add(this).log();
  2990.     }

  2992.     /** {@inheritDoc}
  2993.      */
  2994.     @Override
  2995.     public Dfp atanh() {
  2996.         return getOne().add(this).divide(getOne().subtract(this)).log().divide(2);
  2997.     }

  2999.     /** {@inheritDoc} */
  3000.     @Override
  3001.     public Dfp toDegrees() {
  3002.         return multiply(field.getRadToDeg());
  3003.     }

  3005.     /** {@inheritDoc} */
  3006.     @Override
  3007.     public Dfp toRadians() {
  3008.         return multiply(field.getDegToRad());
  3009.     }

  3011.     /** {@inheritDoc}
  3012.      */
  3013.     @Override
  3014.     public Dfp linearCombination(final Dfp[] a, final Dfp[] b)
  3015.         throws MathIllegalArgumentException {

  3017.         MathUtils.checkDimension(a.length, b.length);

  3019.         // compute in extended accuracy
  3020.         final DfpField extendedField = a[0].field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3021.         Dfp r = extendedField.getZero();
  3022.         for (int i = 0; i < a.length; ++i) {
  3023.             final Dfp aiExt = a[i].newInstance(extendedField, null);
  3024.             final Dfp biExt = b[i].newInstance(extendedField, null);
  3025.             r = r.add(aiExt.multiply(biExt));
  3026.         }

  3028.         // back to normal accuracy
  3029.         return r.newInstance(a[0].field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3031.     }

  3033.     /** {@inheritDoc}
  3034.      */
  3035.     @Override
  3036.     public Dfp linearCombination(final double[] a, final Dfp[] b)
  3037.         throws MathIllegalArgumentException {

  3039.         MathUtils.checkDimension(a.length, b.length);

  3041.         // compute in extended accuracy
  3042.         final DfpField extendedField = b[0].field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3043.         Dfp r = extendedField.getZero();
  3044.         for (int i = 0; i < a.length; ++i) {
  3045.             final Dfp biExt = b[i].newInstance(extendedField, null);
  3046.             r = r.add(biExt.multiply(a[i]));
  3047.         }

  3049.         // back to normal accuracy
  3050.         return r.newInstance(b[0].field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3052.     }

  3054.     /** {@inheritDoc}
  3055.      */
  3056.     @Override
  3057.     public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2) {

  3059.         // switch to extended accuracy
  3060.         final DfpField extendedField = a1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3061.         final Dfp a1Ext = a1.newInstance(extendedField, null);
  3062.         final Dfp b1Ext = b1.newInstance(extendedField, null);
  3063.         final Dfp a2Ext = a2.newInstance(extendedField, null);
  3064.         final Dfp b2Ext = b2.newInstance(extendedField, null);

  3066.         // compute linear combination in extended accuracy
  3067.         final Dfp resultExt = a1Ext.multiply(b1Ext).
  3068.                           add(a2Ext.multiply(b2Ext));

  3070.         // back to normal accuracy
  3071.         return resultExt.newInstance(a1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3073.     }

  3075.     /** {@inheritDoc}
  3076.      */
  3077.     @Override
  3078.     public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2) {

  3080.         // switch to extended accuracy
  3081.         final DfpField extendedField = b1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3082.         final Dfp b1Ext = b1.newInstance(extendedField, null);
  3083.         final Dfp b2Ext = b2.newInstance(extendedField, null);

  3085.         // compute linear combination in extended accuracy
  3086.         final Dfp resultExt = b1Ext.multiply(a1).
  3087.                           add(b2Ext.multiply(a2));

  3089.         // back to normal accuracy
  3090.         return resultExt.newInstance(b1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3092.     }

  3094.     /** {@inheritDoc}
  3095.      */
  3096.     @Override
  3097.     public Dfp linearCombination(final Dfp a1, final Dfp b1,
  3098.                                  final Dfp a2, final Dfp b2,
  3099.                                  final Dfp a3, final Dfp b3) {

  3101.         // switch to extended accuracy
  3102.         final DfpField extendedField = a1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3103.         final Dfp a1Ext = a1.newInstance(extendedField, null);
  3104.         final Dfp b1Ext = b1.newInstance(extendedField, null);
  3105.         final Dfp a2Ext = a2.newInstance(extendedField, null);
  3106.         final Dfp b2Ext = b2.newInstance(extendedField, null);
  3107.         final Dfp a3Ext = a3.newInstance(extendedField, null);
  3108.         final Dfp b3Ext = b3.newInstance(extendedField, null);

  3110.         // compute linear combination in extended accuracy
  3111.         final Dfp resultExt = a1Ext.multiply(b1Ext).
  3112.                           add(a2Ext.multiply(b2Ext)).
  3113.                           add(a3Ext.multiply(b3Ext));

  3115.         // back to normal accuracy
  3116.         return resultExt.newInstance(a1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3118.     }

  3120.     /** {@inheritDoc}
  3121.      */
  3122.     @Override
  3123.     public Dfp linearCombination(final double a1, final Dfp b1,
  3124.                                  final double a2, final Dfp b2,
  3125.                                  final double a3, final Dfp b3) {

  3127.         // switch to extended accuracy
  3128.         final DfpField extendedField = b1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3129.         final Dfp b1Ext = b1.newInstance(extendedField, null);
  3130.         final Dfp b2Ext = b2.newInstance(extendedField, null);
  3131.         final Dfp b3Ext = b3.newInstance(extendedField, null);

  3133.         // compute linear combination in extended accuracy
  3134.         final Dfp resultExt = b1Ext.multiply(a1).
  3135.                           add(b2Ext.multiply(a2)).
  3136.                           add(b3Ext.multiply(a3));

  3138.         // back to normal accuracy
  3139.         return resultExt.newInstance(b1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3141.     }

  3143.     /** {@inheritDoc}
  3144.      */
  3145.     @Override
  3146.     public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2,
  3147.                                  final Dfp a3, final Dfp b3, final Dfp a4, final Dfp b4) {

  3149.         // switch to extended accuracy
  3150.         final DfpField extendedField = a1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3151.         final Dfp a1Ext = a1.newInstance(extendedField, null);
  3152.         final Dfp b1Ext = b1.newInstance(extendedField, null);
  3153.         final Dfp a2Ext = a2.newInstance(extendedField, null);
  3154.         final Dfp b2Ext = b2.newInstance(extendedField, null);
  3155.         final Dfp a3Ext = a3.newInstance(extendedField, null);
  3156.         final Dfp b3Ext = b3.newInstance(extendedField, null);
  3157.         final Dfp a4Ext = a4.newInstance(extendedField, null);
  3158.         final Dfp b4Ext = b4.newInstance(extendedField, null);

  3160.         // compute linear combination in extended accuracy
  3161.         final Dfp resultExt = a1Ext.multiply(b1Ext).
  3162.                           add(a2Ext.multiply(b2Ext)).
  3163.                           add(a3Ext.multiply(b3Ext)).
  3164.                           add(a4Ext.multiply(b4Ext));

  3166.         // back to normal accuracy
  3167.         return resultExt.newInstance(a1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3169.     }

  3171.     /** {@inheritDoc}
  3172.      */
  3173.     @Override
  3174.     public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2,
  3175.                                  final double a3, final Dfp b3, final double a4, final Dfp b4) {

  3177.         // switch to extended accuracy
  3178.         final DfpField extendedField = b1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
  3179.         final Dfp b1Ext = b1.newInstance(extendedField, null);
  3180.         final Dfp b2Ext = b2.newInstance(extendedField, null);
  3181.         final Dfp b3Ext = b3.newInstance(extendedField, null);
  3182.         final Dfp b4Ext = b4.newInstance(extendedField, null);

  3184.         // compute linear combination in extended accuracy
  3185.         final Dfp resultExt = b1Ext.multiply(a1).
  3186.                           add(b2Ext.multiply(a2)).
  3187.                           add(b3Ext.multiply(a3)).
  3188.                           add(b4Ext.multiply(a4));

  3190.         // back to normal accuracy
  3191.         return resultExt.newInstance(b1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);

  3193.     }

  3195.     /** {@inheritDoc} */
  3196.     @Override
  3197.     public Dfp getPi() {
  3198.         return field.getPi();
  3199.     }

  3201. }
Created with JaCoCo 0.8.12.202403310830