/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      https://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  
  /*
   * This is not the original file distributed by the Apache Software Foundation
   * It has been modified by the Hipparchus project
   */
  
  package org.hipparchus.dfp;
  
  import java.util.Arrays;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.exception.MathIllegalArgumentException;
  import org.hipparchus.exception.MathRuntimeException;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.FieldSinhCosh;
  import org.hipparchus.util.MathUtils;
  
  /**
   *  Decimal floating point library for Java
   *
   *  <p>Another floating point class.  This one is built using radix 10000
   *  which is 10<sup>4</sup>, so its almost decimal.</p>
   *
   *  <p>The design goals here are:</p>
   *  <ol>
   *    <li>Decimal math, or close to it</li>
   *    <li>Settable precision (but no mix between numbers using different settings)</li>
   *    <li>Portability.  Code should be kept as portable as possible.</li>
   *    <li>Performance</li>
   *    <li>Accuracy  - Results should always be +/- 1 ULP for basic
   *         algebraic operation</li>
   *    <li>Comply with IEEE 854-1987 as much as possible.
   *         (See IEEE 854-1987 notes below)</li>
   *  </ol>
   *
   *  <p>Trade offs:</p>
   *  <ol>
   *    <li>Memory foot print.  I'm using more memory than necessary to
   *         represent numbers to get better performance.</li>
   *    <li>Digits are bigger, so rounding is a greater loss.  So, if you
   *         really need 12 decimal digits, better use 4 base 10000 digits
   *         there can be one partially filled.</li>
   *  </ol>
   *
   *  <p>Numbers are represented  in the following form:
   *  \[
   *  n  =  \mathrm{sign} \times \mathrm{mant} \times \mathrm{radix}^\mathrm{exp}
   *  \]
   *  where sign is &plusmn;1, mantissa represents a fractional number between
   *  zero and one.  mant[0] is the least significant digit.
   *  exp is in the range of -32767 to 32768</p>
   *
   *  <p>IEEE 854-1987  Notes and differences</p>
   *
   *  <p>IEEE 854 requires the radix to be either 2 or 10.  The radix here is
   *  10000, so that requirement is not met, but  it is possible that a
   *  subclassed can be made to make it behave as a radix 10
   *  number.  It is my opinion that if it looks and behaves as a radix
   *  10 number then it is one and that requirement would be met.</p>
   *
   *  <p>The radix of 10000 was chosen because it should be faster to operate
   *  on 4 decimal digits at once instead of one at a time.  Radix 10 behavior
   *  can be realized by adding an additional rounding step to ensure that
   *  the number of decimal digits represented is constant.</p>
   *
   *  <p>The IEEE standard specifically leaves out internal data encoding,
   *  so it is reasonable to conclude that such a subclass of this radix
   *  10000 system is merely an encoding of a radix 10 system.</p>
   *
   *  <p>IEEE 854 also specifies the existence of "sub-normal" numbers.  This
   *  class does not contain any such entities.  The most significant radix
   *  10000 digit is always non-zero.  Instead, we support "gradual underflow"
   *  by raising the underflow flag for numbers less with exponent less than
   *  expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
   *  Thus the smallest number we can represent would be:
   *  1E(-(MIN_EXP-digits-1)*4),  eg, for digits=5, MIN_EXP=-32767, that would
   *  be 1e-131092.</p>
   *
   *  <p>IEEE 854 defines that the implied radix point lies just to the right
   *  of the most significant digit and to the left of the remaining digits.
   *  This implementation puts the implied radix point to the left of all
   *  digits including the most significant one.  The most significant digit
   *  here is the one just to the right of the radix point.  This is a fine
  *  detail and is really only a matter of definition.  Any side effects of
  *  this can be rendered invisible by a subclass.</p>
  * @see DfpField
  */
 public class Dfp implements CalculusFieldElement<Dfp> {
 
     /** The radix, or base of this system.  Set to 10000 */
     public static final int RADIX = 10000;
 
     /** The minimum exponent before underflow is signaled.  Flush to zero
      *  occurs at minExp-DIGITS */
     public static final int MIN_EXP = -32767;
 
     /** The maximum exponent before overflow is signaled and results flushed
      *  to infinity */
     public static final int MAX_EXP =  32768;
 
     /** The amount under/overflows are scaled by before going to trap handler */
     public static final int ERR_SCALE = 32760;
 
     /** Indicator value for normal finite numbers. */
     public static final byte FINITE = 0;
 
     /** Indicator value for Infinity. */
     public static final byte INFINITE = 1;
 
     /** Indicator value for signaling NaN. */
     public static final byte SNAN = 2;
 
     /** Indicator value for quiet NaN. */
     public static final byte QNAN = 3;
 
     /** String for NaN representation. */
     private static final String NAN_STRING = "NaN";
 
     /** String for positive infinity representation. */
     private static final String POS_INFINITY_STRING = "Infinity";
 
     /** String for negative infinity representation. */
     private static final String NEG_INFINITY_STRING = "-Infinity";
 
     /** Name for traps triggered by addition. */
     private static final String ADD_TRAP = "add";
 
     /** Name for traps triggered by multiplication. */
     private static final String MULTIPLY_TRAP = "multiply";
 
     /** Name for traps triggered by division. */
     private static final String DIVIDE_TRAP = "divide";
 
     /** Name for traps triggered by square root. */
     private static final String SQRT_TRAP = "sqrt";
 
     /** Name for traps triggered by alignment. */
     private static final String ALIGN_TRAP = "align";
 
     /** Name for traps triggered by truncation. */
     private static final String TRUNC_TRAP = "trunc";
 
     /** Name for traps triggered by nextAfter. */
     private static final String NEXT_AFTER_TRAP = "nextAfter";
 
     /** Name for traps triggered by lessThan. */
     private static final String LESS_THAN_TRAP = "lessThan";
 
     /** Name for traps triggered by greaterThan. */
     private static final String GREATER_THAN_TRAP = "greaterThan";
 
     /** Name for traps triggered by newInstance. */
     private static final String NEW_INSTANCE_TRAP = "newInstance";
 
     /** Multiplication factor for number of digits used to compute linear combinations. */
     private static final int LINEAR_COMBINATION_DIGITS_FACTOR = 2;
 
     /** Mantissa. */
     protected int[] mant;
 
     /** Sign bit: 1 for positive, -1 for negative. */
     protected byte sign;
 
     /** Exponent. */
     protected int exp;
 
     /** Indicator for non-finite / non-number values. */
     protected byte nans;
 
     /** Factory building similar Dfp's. */
     private final DfpField field;
 
     /** Makes an instance with a value of zero.
      * @param field field to which this instance belongs
      */
     protected Dfp(final DfpField field) {
         mant = new int[field.getRadixDigits()];
         sign = 1;
         exp = 0;
         nans = FINITE;
         this.field = field;
     }
 
     /** Create an instance from a byte value.
      * @param field field to which this instance belongs
      * @param x value to convert to an instance
      */
     protected Dfp(final DfpField field, byte x) {
         this(field, (long) x);
     }
 
     /** Create an instance from an int value.
      * @param field field to which this instance belongs
      * @param x value to convert to an instance
      */
     protected Dfp(final DfpField field, int x) {
         this(field, (long) x);
     }
 
     /** Create an instance from a long value.
      * @param field field to which this instance belongs
      * @param x value to convert to an instance
      */
     protected Dfp(final DfpField field, long x) {
 
         // initialize as if 0
         mant = new int[field.getRadixDigits()];
         nans = FINITE;
         this.field = field;
 
         boolean isLongMin = false;
         if (x == Long.MIN_VALUE) {
             // special case for Long.MIN_VALUE (-9223372036854775808)
             // we must shift it before taking its absolute value
             isLongMin = true;
             ++x;
         }
 
         // set the sign
         if (x < 0) {
             sign = -1;
             x = -x;
         } else {
             sign = 1;
         }
 
         exp = 0;
         while (x != 0) {
             System.arraycopy(mant, mant.length - exp, mant, mant.length - 1 - exp, exp);
             mant[mant.length - 1] = (int) (x % RADIX);
             x /= RADIX;
             exp++;
         }
 
         if (isLongMin) {
             // remove the shift added for Long.MIN_VALUE
             // we know in this case that fixing the last digit is sufficient
             for (int i = 0; i < mant.length - 1; i++) {
                 if (mant[i] != 0) {
                     mant[i]++;
                     break;
                 }
             }
         }
     }
 
     /** Create an instance from a double value.
      * @param field field to which this instance belongs
      * @param x value to convert to an instance
      */
     protected Dfp(final DfpField field, double x) {
 
         // initialize as if 0
         mant = new int[field.getRadixDigits()];
         this.field = field;
 
         long bits = Double.doubleToLongBits(x);
         long mantissa = bits & 0x000fffffffffffffL;
         int exponent = (int) ((bits & 0x7ff0000000000000L) >> 52) - 1023;
 
         if (exponent == -1023) {
             // Zero or sub-normal
             if (x == 0) {
                 // make sure 0 has the right sign
                 if ((bits & 0x8000000000000000L) != 0) {
                     sign = -1;
                 } else {
                     sign = 1;
                 }
                 return;
             }
 
             exponent++;
 
             // Normalize the subnormal number
             while ( (mantissa & 0x0010000000000000L) == 0) {
                 exponent--;
                 mantissa <<= 1;
             }
             mantissa &= 0x000fffffffffffffL;
         }
 
         if (exponent == 1024) {
             // infinity or NAN
             if (x != x) {
                 sign = (byte) 1;
                 nans = QNAN;
             } else if (x < 0) {
                 sign = (byte) -1;
                 nans = INFINITE;
             } else {
                 sign = (byte) 1;
                 nans = INFINITE;
             }
             return;
         }
 
         Dfp xdfp = new Dfp(field, mantissa);
         xdfp = xdfp.divide(new Dfp(field, 4503599627370496L)).add(field.getOne());  // Divide by 2^52, then add one
         xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));
 
         if ((bits & 0x8000000000000000L) != 0) {
             xdfp = xdfp.negate();
         }
 
         System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
         sign = xdfp.sign;
         exp  = xdfp.exp;
         nans = xdfp.nans;
 
     }
 
     /** Copy constructor.
      * @param d instance to copy
      */
     public Dfp(final Dfp d) {
         mant  = d.mant.clone();
         sign  = d.sign;
         exp   = d.exp;
         nans  = d.nans;
         field = d.field;
     }
 
     /** Create an instance from a String representation.
      * @param field field to which this instance belongs
      * @param s string representation of the instance
      */
     protected Dfp(final DfpField field, final String s) {
 
         // initialize as if 0
         mant = new int[field.getRadixDigits()];
         sign = 1;
         nans = FINITE;
         this.field = field;
 
         boolean decimalFound = false;
         final int rsize = 4;   // size of radix in decimal digits
         final int offset = 4;  // Starting offset into Striped
         final char[] striped = new char[getRadixDigits() * rsize + offset * 2];
 
         // Check some special cases
         if (POS_INFINITY_STRING.equals(s)) {
             sign = (byte) 1;
             nans = INFINITE;
             return;
         }
 
         if (NEG_INFINITY_STRING.equals(s)) {
             sign = (byte) -1;
             nans = INFINITE;
             return;
         }
 
         if (NAN_STRING.equals(s)) {
             sign = (byte) 1;
             nans = QNAN;
             return;
         }
 
         // Check for scientific notation
         int p = s.indexOf('e');
         if (p == -1) { // try upper case?
             p = s.indexOf('E');
         }
 
         final String fpdecimal;
         int sciexp = 0;
         if (p != -1) {
             // scientific notation
             fpdecimal = s.substring(0, p);
             String fpexp = s.substring(p+1);
             boolean negative = false;
 
             for (int i=0; i<fpexp.length(); i++)
             {
                 if (fpexp.charAt(i) == '-')
                 {
                     negative = true;
                     continue;
                 }
                 if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') {
                     sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
                 }
             }
 
             if (negative) {
                 sciexp = -sciexp;
             }
         } else {
             // normal case
             fpdecimal = s;
         }
 
         // If there is a minus sign in the number then it is negative
         if (fpdecimal.indexOf('-') !=  -1) {
             sign = -1;
         }
 
         // First off, find all of the leading zeros, trailing zeros, and significant digits
         p = 0;
 
         // Move p to first significant digit
         int decimalPos = 0;
         for (;;) {
             if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
                 break;
             }
 
             if (decimalFound && fpdecimal.charAt(p) == '0') {
                 decimalPos--;
             }
 
             if (fpdecimal.charAt(p) == '.') {
                 decimalFound = true;
             }
 
             p++;
 
             if (p == fpdecimal.length()) {
                 break;
             }
         }
 
         // Copy the string onto Stripped
         int q = offset;
         striped[0] = '0';
         striped[1] = '0';
         striped[2] = '0';
         striped[3] = '0';
         int significantDigits=0;
         for(;;) {
             if (p == (fpdecimal.length())) {
                 break;
             }
 
             // Don't want to run pass the end of the array
             if (q == mant.length*rsize+offset+1) {
                 break;
             }
 
             if (fpdecimal.charAt(p) == '.') {
                 decimalFound = true;
                 decimalPos = significantDigits;
                 p++;
                 continue;
             }
 
             if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
                 p++;
                 continue;
             }
 
             striped[q] = fpdecimal.charAt(p);
             q++;
             p++;
             significantDigits++;
         }
 
 
         // If the decimal point has been found then get rid of trailing zeros.
         if (decimalFound && q != offset) {
             for (;;) {
                 q--;
                 if (q == offset) {
                     break;
                 }
                 if (striped[q] == '0') {
                     significantDigits--;
                 } else {
                     break;
                 }
             }
         }
 
         // special case of numbers like "0.00000"
         if (decimalFound && significantDigits == 0) {
             decimalPos = 0;
         }
 
         // Implicit decimal point at end of number if not present
         if (!decimalFound) {
             decimalPos = q-offset;
         }
 
         // Find the number of significant trailing zeros
         q = offset;  // set q to point to first sig digit
         p = significantDigits-1+offset;
 
         while (p > q) {
             if (striped[p] != '0') {
                 break;
             }
             p--;
         }
 
         // Make sure the decimal is on a mod 10000 boundary
         int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
         q -= i;
         decimalPos += i;
 
         // Make the mantissa length right by adding zeros at the end if necessary
         while ((p - q) < (mant.length * rsize)) {
             for (i = 0; i < rsize; i++) {
                 striped[++p] = '0';
             }
         }
 
         // Ok, now we know how many trailing zeros there are,
         // and where the least significant digit is
         for (i = mant.length - 1; i >= 0; i--) {
             mant[i] = (striped[q]   - '0') * 1000 +
                       (striped[q+1] - '0') * 100  +
                       (striped[q+2] - '0') * 10   +
                       (striped[q+3] - '0');
             q += 4;
         }
 
         exp = (decimalPos+sciexp) / rsize;
 
         if (q < striped.length) {
             // Is there possible another digit?
             round((striped[q] - '0')*1000);
         }
 
     }
 
     /** Creates an instance with a non-finite value.
      * @param field field to which this instance belongs
      * @param sign sign of the Dfp to create
      * @param nans code of the value, must be one of {@link #INFINITE},
      * {@link #SNAN},  {@link #QNAN}
      */
     protected Dfp(final DfpField field, final byte sign, final byte nans) {
         this.field = field;
         this.mant    = new int[field.getRadixDigits()];
         this.sign    = sign;
         this.exp     = 0;
         this.nans    = nans;
     }
 
     /** Create an instance with a value of 0.
      * Use this internally in preference to constructors to facilitate subclasses
      * @return a new instance with a value of 0
      */
     public Dfp newInstance() {
         return new Dfp(getField());
     }
 
     /** Create an instance from a byte value.
      * @param x value to convert to an instance
      * @return a new instance with value x
      */
     public Dfp newInstance(final byte x) {
         return new Dfp(getField(), x);
     }
 
     /** Create an instance from an int value.
      * @param x value to convert to an instance
      * @return a new instance with value x
      */
     public Dfp newInstance(final int x) {
         return new Dfp(getField(), x);
     }
 
     /** Create an instance from a long value.
      * @param x value to convert to an instance
      * @return a new instance with value x
      */
     public Dfp newInstance(final long x) {
         return new Dfp(getField(), x);
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp newInstance(final double x) {
         return new Dfp(getField(), x);
     }
 
     /** Create an instance by copying an existing one.
      * Use this internally in preference to constructors to facilitate subclasses.
      * @param d instance to copy
      * @return a new instance with the same value as d
      */
     public Dfp newInstance(final Dfp d) {
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != d.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
         }
 
         return new Dfp(d);
 
     }
 
     /** Create an instance from a String representation.
      * Use this internally in preference to constructors to facilitate subclasses.
      * @param s string representation of the instance
      * @return a new instance parsed from specified string
      */
     public Dfp newInstance(final String s) {
         return new Dfp(field, s);
     }
 
     /** Creates an instance with a non-finite value.
      * @param sig sign of the Dfp to create
      * @param code code of the value, must be one of {@link #INFINITE},
      * {@link #SNAN},  {@link #QNAN}
      * @return a new instance with a non-finite value
      */
     public Dfp newInstance(final byte sig, final byte code) {
         return field.newDfp(sig, code);
     }
 
     /** Creates an instance by converting the instance to a different field (i.e. different accuracy).
      * <p>
      * If the target field as a greater number of digits, the extra least significant digits
      * will be set to zero.
      * </p>
      * @param targetField field to convert the instance to
      * @param rmode rounding mode to use if target field as less digits than the instance, can be null otherwise
      * @return converted instance (or the instance itself if it already has the required number of digits)
      * @see DfpField#getExtendedField(int, boolean)
      * @since 1.7
      */
     public Dfp newInstance(final DfpField targetField, final DfpField.RoundingMode rmode) {
         final int deltaLength = targetField.getRadixDigits() - field.getRadixDigits();
         if (deltaLength == 0) {
 
             // no conversion, we return the instance itself
             return this;
 
         } else {
 
             // create an instance (initially set to 0) with the expected number of digits
             Dfp result = new Dfp(targetField);
             result.sign = sign;
             result.exp  = exp;
             result.nans = nans;
             if (nans == 0) {
 
                 if (deltaLength < 0) {
 
                     // copy only the most significant digits, dropping the least significant ones
                     // the result corresponds to pure truncation, proper rounding will follow
                     System.arraycopy(mant, -deltaLength, result.mant, 0, result.mant.length);
 
                     // check if we have dropped any non-zero digits in the low part
                     // (not counting the last dropped digit which will be handled specially)
                     final int last = -(deltaLength + 1);
                     boolean zeroLSB = true;
                     for (int i = 0; i < last; ++i) {
                         zeroLSB &= mant[i] == 0;
                     }
 
                     if (!(zeroLSB && mant[last] == 0)) {
                         // there are some non-zero digits that have been discarded, perform rounding
 
                         if (shouldIncrement(rmode, zeroLSB, mant[last], result.mant[0], sign)) {
                             // rounding requires incrementing the mantissa
                             result.incrementMantissa();
                         }
 
                         targetField.setIEEEFlagsBits(DfpField.FLAG_INEXACT);  // signal inexact
                         result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
 
                     }
 
                 } else {
                     // copy all digits as the new most significant ones, leaving the least significant digits to zero
                     System.arraycopy(mant, 0, result.mant, deltaLength, mant.length);
                 }
 
             }
 
             return result;
 
         }
     }
 
     /** Check if mantissa of a truncated number must be incremented.
      * <p>
      * This method must be called <em>only</em> when some non-zero digits have been
      * discarded (i.e. when either {@code zeroLSB} is false or {@code lastDiscarded} is non-zero),
      * otherwise it would generate false positive
      * </p>
      * @param rmode rounding mode to use if target field as less digits than the instance, can be null otherwise
      * @param zeroLSB true is least significant discarded digits (except last) are all zero
      * @param lastDiscarded last discarded digit
      * @param firstNonDiscarded first non-discarded digit
      * @param sign of the number
      * @return true if the already truncated mantissa should be incremented to achieve correct rounding
      * @since 1.7
      */
     private static boolean shouldIncrement(final DfpField.RoundingMode rmode,
                                            final boolean zeroLSB, final int lastDiscarded,
                                            final int firstNonDiscarded, final int sign) {
         switch (rmode) {
             case ROUND_DOWN :
                 return false;
 
             case ROUND_UP :
                 return true;
 
             case ROUND_HALF_UP :
                 return lastDiscarded >= 5000;
 
             case ROUND_HALF_DOWN :
                 return isAboveHalfWay(zeroLSB, lastDiscarded);
 
             case ROUND_HALF_EVEN :
                 return (isHalfWay(zeroLSB, lastDiscarded) && (firstNonDiscarded & 0x1) == 0x1) ||
                        isAboveHalfWay(zeroLSB, lastDiscarded);
 
             case ROUND_HALF_ODD :
                 return (isHalfWay(zeroLSB, lastDiscarded) && (firstNonDiscarded & 0x1) == 0x0) ||
                        isAboveHalfWay(zeroLSB, lastDiscarded);
 
             case ROUND_CEIL :
                 return sign > 0;
 
             case ROUND_FLOOR :
                 return sign < 0;
 
             default :
                 // this should never happen
                 throw MathRuntimeException.createInternalError();
         }
     }
 
     /** Increment the mantissa of the instance
      * @since 1.7
      */
     private void incrementMantissa() {
         boolean carry = true;
         for (int i = 0; carry && i < mant.length; ++i) {
             ++mant[i];
             if (mant[i] >= RADIX) {
                 mant[i] -= RADIX;
             } else {
                 carry = false;
             }
         }
         if (carry) {
             // we have exceeded capacity, we need to drop one digit
             for (int i = 0; i < mant.length - 1; i++) {
                 mant[i] = mant[i+1];
             }
             mant[mant.length - 1] = 1;
             exp++;
         }
     }
 
     /** Check if discarded digits are exactly halfway between two rounder numbers.
      * @param zeroLSB true is least significant discarded digits (except last) are all zero
      * @param lastDiscarded last discarded digit
      * @return true if discarded digits correspond to a number exactly halfway between two rounded numbers
      * @since 1.7
      */
     private static boolean isHalfWay(final boolean zeroLSB, final int lastDiscarded) {
         return lastDiscarded == 5000 && zeroLSB;
     }
 
     /** Check if discarded digits are strictly above halfway between two rounder numbers.
      * @param zeroLSB true is least significant discarded digits (except last) are all zero
      * @param lastDiscarded last discarded digit
      * @return true if discarded digits correspond to a number strictly above halfway between two rounded numbers
      * @since 1.7
      */
     private static boolean isAboveHalfWay(final boolean zeroLSB, final int lastDiscarded) {
         return (lastDiscarded > 5000) || (lastDiscarded == 5000 && !zeroLSB);
     }
 
     /** Get the {@link org.hipparchus.Field Field} (really a {@link DfpField}) to which the instance belongs.
      * <p>
      * The field is linked to the number of digits and acts as a factory
      * for {@link Dfp} instances.
      * </p>
      * @return {@link org.hipparchus.Field Field} (really a {@link DfpField}) to which the instance belongs
      */
     @Override
     public DfpField getField() {
         return field;
     }
 
     /** Get the number of radix digits of the instance.
      * @return number of radix digits
      */
     public int getRadixDigits() {
         return field.getRadixDigits();
     }
 
     /** Get the constant 0.
      * @return a Dfp with value zero
      */
     public Dfp getZero() {
         return field.getZero();
     }
 
     /** Get the constant 1.
      * @return a Dfp with value one
      */
     public Dfp getOne() {
         return field.getOne();
     }
 
     /** Get the constant 2.
      * @return a Dfp with value two
      */
     public Dfp getTwo() {
         return field.getTwo();
     }
 
     /** Shift the mantissa left, and adjust the exponent to compensate.
      */
     protected void shiftLeft() {
         for (int i = mant.length - 1; i > 0; i--) {
             mant[i] = mant[i-1];
         }
         mant[0] = 0;
         exp--;
     }
 
     /* Note that shiftRight() does not call round() as that round() itself
      uses shiftRight() */
     /** Shift the mantissa right, and adjust the exponent to compensate.
      */
     protected void shiftRight() {
         for (int i = 0; i < mant.length - 1; i++) {
             mant[i] = mant[i+1];
         }
         mant[mant.length - 1] = 0;
         exp++;
     }
 
     /** Make our exp equal to the supplied one, this may cause rounding.
      *  Also causes de-normalized numbers.  These numbers are generally
      *  dangerous because most routines assume normalized numbers.
      *  Align doesn't round, so it will return the last digit destroyed
      *  by shifting right.
      *  @param e desired exponent
      *  @return last digit destroyed by shifting right
      */
     protected int align(int e) {
         int lostdigit = 0;
         boolean inexact = false;
 
         int diff = exp - e;
 
         int adiff = diff;
         if (adiff < 0) {
             adiff = -adiff;
         }
 
         if (diff == 0) {
             return 0;
         }
 
         if (adiff > (mant.length + 1)) {
             // Special case
             Arrays.fill(mant, 0);
             exp = e;
 
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
             dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
 
             return 0;
         }
 
         for (int i = 0; i < adiff; i++) {
             if (diff < 0) {
                 /* Keep track of loss -- only signal inexact after losing 2 digits.
                  * the first lost digit is returned to add() and may be incorporated
                  * into the result.
                  */
                 if (lostdigit != 0) {
                     inexact = true;
                 }
 
                 lostdigit = mant[0];
 
                 shiftRight();
             } else {
                 shiftLeft();
             }
         }
 
         if (inexact) {
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
             dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
         }
 
         return lostdigit;
 
     }
 
     /** Check if instance is less than x.
      * @param x number to check instance against
      * @return true if instance is less than x and neither are NaN, false otherwise
      */
     public boolean lessThan(final Dfp x) {
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != x.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
             return false;
         }
 
         /* if a nan is involved, signal invalid and return false */
         if (isNaN() || x.isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
             return false;
         }
 
         return compare(this, x) < 0;
     }
 
     /** Check if instance is greater than x.
      * @param x number to check instance against
      * @return true if instance is greater than x and neither are NaN, false otherwise
      */
     public boolean greaterThan(final Dfp x) {
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != x.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
             return false;
         }
 
         /* if a nan is involved, signal invalid and return false */
         if (isNaN() || x.isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
             return false;
         }
 
         return compare(this, x) > 0;
     }
 
     /** Check if instance is less than or equal to 0.
      * @return true if instance is not NaN and less than or equal to 0, false otherwise
      */
     public boolean negativeOrNull() {
 
         if (isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
             return false;
         }
 
         return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
 
     }
 
     /** Check if instance is strictly less than 0.
      * @return true if instance is not NaN and less than or equal to 0, false otherwise
      */
     public boolean strictlyNegative() {
 
         if (isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
             return false;
         }
 
         return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite());
 
     }
 
     /** Check if instance is greater than or equal to 0.
      * @return true if instance is not NaN and greater than or equal to 0, false otherwise
      */
     public boolean positiveOrNull() {
 
         if (isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
             return false;
         }
 
         return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
 
     }
 
     /** Check if instance is strictly greater than 0.
      * @return true if instance is not NaN and greater than or equal to 0, false otherwise
      */
     public boolean strictlyPositive() {
 
         if (isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
             return false;
         }
 
         return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite());
 
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp abs() {
         Dfp result = newInstance(this);
         result.sign = 1;
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public boolean isInfinite() {
         return nans == INFINITE;
     }
 
     /** {@inheritDoc} */
     @Override
     public boolean isNaN() {
         return (nans == QNAN) || (nans == SNAN);
     }
 
     /** Check if instance is equal to zero.
      * @return true if instance is equal to zero
      */
     @Override
     public boolean isZero() {
 
         if (isNaN()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
             return false;
         }
 
         return (mant[mant.length - 1] == 0) && !isInfinite();
 
     }
 
     /** Check if instance is equal to x.
      * @param other object to check instance against
      * @return true if instance is equal to x and neither are NaN, false otherwise
      */
     @Override
     public boolean equals(final Object other) {
 
         if (other instanceof Dfp) {
             final Dfp x = (Dfp) other;
             if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
                 return false;
             }
 
             return compare(this, x) == 0;
         }
 
         return false;
 
     }
 
     /**
      * Gets a hashCode for the instance.
      * @return a hash code value for this object
      */
     @Override
     public int hashCode() {
         return 17 + (isZero() ? 0 : (sign << 8)) + (nans << 16) + exp + Arrays.hashCode(mant);
     }
 
     /** Check if instance is not equal to x.
      * @param x number to check instance against
      * @return true if instance is not equal to x and neither are NaN, false otherwise
      */
     public boolean unequal(final Dfp x) {
         if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
             return false;
         }
 
         return greaterThan(x) || lessThan(x);
     }
 
     /** Compare two instances.
      * @param a first instance in comparison
      * @param b second instance in comparison
      * @return -1 if a<b, 1 if a>b and 0 if a==b
      *  Note this method does not properly handle NaNs or numbers with different precision.
      */
     private static int compare(final Dfp a, final Dfp b) {
         // Ignore the sign of zero
         if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
             a.nans == FINITE && b.nans == FINITE) {
             return 0;
         }
 
         if (a.sign != b.sign) {
             if (a.sign == -1) {
                 return -1;
             } else {
                 return 1;
             }
         }
 
         // deal with the infinities
         if (a.nans == INFINITE && b.nans == FINITE) {
             return a.sign;
         }
 
         if (a.nans == FINITE && b.nans == INFINITE) {
             return -b.sign;
         }
 
         if (a.nans == INFINITE && b.nans == INFINITE) {
             return 0;
         }
 
         // Handle special case when a or b is zero, by ignoring the exponents
         if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) {
             if (a.exp < b.exp) {
                 return -a.sign;
             }
 
             if (a.exp > b.exp) {
                 return a.sign;
             }
         }
 
         // compare the mantissas
         for (int i = a.mant.length - 1; i >= 0; i--) {
             if (a.mant[i] > b.mant[i]) {
                 return a.sign;
             }
 
             if (a.mant[i] < b.mant[i]) {
                 return -a.sign;
             }
         }
 
         return 0;
 
     }
 
     /** Round to nearest integer using the round-half-even method.
      *  That is round to nearest integer unless both are equidistant.
      *  In which case round to the even one.
      *  @return rounded value
      */
     @Override
     public Dfp rint() {
         return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
     }
 
     /** Round to an integer using the round floor mode.
      * That is, round toward -Infinity
      *  @return rounded value
      */
     @Override
     public Dfp floor() {
         return trunc(DfpField.RoundingMode.ROUND_FLOOR);
     }
 
     /** Round to an integer using the round ceil mode.
      * That is, round toward +Infinity
      *  @return rounded value
      */
     @Override
     public Dfp ceil() {
         return trunc(DfpField.RoundingMode.ROUND_CEIL);
     }
 
     /** Returns the IEEE remainder.
      * @param d divisor
      * @return this less n &times; d, where n is the integer closest to this/d
      */
     @Override
     public Dfp remainder(final Dfp d) {
 
         final Dfp result = this.subtract(this.divide(d).rint().multiply(d));
 
         // IEEE 854-1987 says that if the result is zero, then it carries the sign of this
         if (result.mant[mant.length-1] == 0) {
             result.sign = sign;
         }
 
         return result;
 
     }
 
     /** Does the integer conversions with the specified rounding.
      * @param rmode rounding mode to use
      * @return truncated value
      */
     protected Dfp trunc(final DfpField.RoundingMode rmode) {
         boolean changed = false;
 
         if (isNaN()) {
             return newInstance(this);
         }
 
         if (nans == INFINITE) {
             return newInstance(this);
         }
 
         if (mant[mant.length-1] == 0) {
             // a is zero
             return newInstance(this);
         }
 
         /* If the exponent is less than zero then we can certainly
          * return -1, 0 or +1 depending on sign and rounding mode */
         if (exp < 0) {
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
             final Dfp result;
             if (sign == -1 && rmode == DfpField.RoundingMode.ROUND_FLOOR) {
                 result = newInstance(-1);
             } else if (sign == +1 && rmode == DfpField.RoundingMode.ROUND_CEIL) {
                 result = newInstance(+1);
             } else {
                 // for all other combinations of sign and mode, zero is the correct rounding
                 result = newInstance(0);
             }
             return dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
         }
 
         /* If the exponent is greater than or equal to digits, then it
          * must already be an integer since there is no precision left
          * for any fractional part */
 
         if (exp >= mant.length) {
             return newInstance(this);
         }
 
         /* General case:  create another dfp, result, that contains the
          * a with the fractional part lopped off.  */
 
         Dfp result = newInstance(this);
         for (int i = 0; i < mant.length-result.exp; i++) {
             changed |= result.mant[i] != 0;
             result.mant[i] = 0;
         }
 
         if (changed) {
             switch (rmode) {
                 case ROUND_FLOOR:
                     if (result.sign == -1) {
                         // then we must increment the mantissa by one
                         result = result.add(newInstance(-1));
                     }
                     break;
 
                 case ROUND_CEIL:
                     if (result.sign == 1) {
                         // then we must increment the mantissa by one
                         result = result.add(getOne());
                     }
                     break;
 
                 case ROUND_HALF_EVEN:
                 default:
                     final Dfp half = newInstance("0.5");
                     Dfp a = subtract(result);  // difference between this and result
                     a.sign = 1;            // force positive (take abs)
                     if (a.greaterThan(half)) {
                         a = newInstance(getOne());
                         a.sign = sign;
                         result = result.add(a);
                     }
 
                     // If exactly equal to 1/2 and odd then increment
                     if (a.equals(half) && result.exp > 0 && (result.mant[mant.length-result.exp]&1) != 0) {
                         a = newInstance(getOne());
                         a.sign = sign;
                         result = result.add(a);
                     }
                     break;
             }
 
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);  // signal inexact
             result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
             return result;
         }
 
         return result;
     }
 
     /** Convert this to an integer.
      * If greater than 2147483647, it returns 2147483647. If less than -2147483648 it returns -2147483648.
      * @return converted number
      */
     public int intValue() {
         Dfp rounded;
         int result = 0;
 
         rounded = rint();
 
         if (rounded.greaterThan(newInstance(2147483647))) {
             return 2147483647;
         }
 
         if (rounded.lessThan(newInstance(-2147483648))) {
             return -2147483648;
         }
 
         for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
             result = result * RADIX + rounded.mant[i];
         }
 
         if (rounded.sign == -1) {
             result = -result;
         }
 
         return result;
     }
 
     /** Get the exponent of the greatest power of 10000 that is
      *  less than or equal to the absolute value of this.  I.E.  if
      *  this is 10<sup>6</sup> then log10K would return 1.
      *  @return integer base 10000 logarithm
      */
     public int log10K() {
         return exp - 1;
     }
 
     /** Get the specified  power of 10000.
      * @param e desired power
      * @return 10000<sup>e</sup>
      */
     public Dfp power10K(final int e) {
         Dfp d = newInstance(getOne());
         d.exp = e + 1;
         return d;
     }
 
     /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
      *  @return integer base 10 logarithm
      */
     public int intLog10()  {
         if (mant[mant.length-1] > 1000) {
             return exp * 4 - 1;
         }
         if (mant[mant.length-1] > 100) {
             return exp * 4 - 2;
         }
         if (mant[mant.length-1] > 10) {
             return exp * 4 - 3;
         }
         return exp * 4 - 4;
     }
 
     /** Return the specified  power of 10.
      * @param e desired power
      * @return 10<sup>e</sup>
      */
     public Dfp power10(final int e) {
         Dfp d = newInstance(getOne());
 
         if (e >= 0) {
             d.exp = e / 4 + 1;
         } else {
             d.exp = (e + 1) / 4;
         }
 
         switch ((e % 4 + 4) % 4) {
             case 0:
                 break;
             case 1:
                 d = d.multiply(10);
                 break;
             case 2:
                 d = d.multiply(100);
                 break;
             default:
                 d = d.multiply(1000);
                 break;
         }
 
         return d;
     }
 
     /** Negate the mantissa of this by computing the complement.
      *  Leaves the sign bit unchanged, used internally by add.
      *  Denormalized numbers are handled properly here.
      *  @param extra ???
      *  @return ???
      */
     protected int complement(int extra) {
 
         extra = RADIX-extra;
         for (int i = 0; i < mant.length; i++) {
             mant[i] = RADIX-mant[i]-1;
         }
 
         int rh = extra / RADIX;
         extra -= rh * RADIX;
         for (int i = 0; i < mant.length; i++) {
             final int r = mant[i] + rh;
             rh = r / RADIX;
             mant[i] = r - rh * RADIX;
         }
 
         return extra;
     }
 
     /** Add x to this.
      * @param x number to add
      * @return sum of this and x
      */
     @Override
     public Dfp add(final Dfp x) {
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != x.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
         }
 
         /* handle special cases */
         if (nans != FINITE || x.nans != FINITE) {
             if (isNaN()) {
                 return this;
             }
 
             if (x.isNaN()) {
                 return x;
             }
 
             if (nans == INFINITE && x.nans == FINITE) {
                 return this;
             }
 
             if (x.nans == INFINITE && nans == FINITE) {
                 return x;
             }
 
             if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
                 return x;
             }
 
             if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
                 Dfp result = newInstance(getZero());
                 result.nans = QNAN;
                 result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
                 return result;
             }
         }
 
         /* copy this and the arg */
         Dfp a = newInstance(this);
         Dfp b = newInstance(x);
 
         /* initialize the result object */
         Dfp result = newInstance(getZero());
 
         /* Make all numbers positive, but remember their sign */
         final byte asign = a.sign;
         final byte bsign = b.sign;
 
         a.sign = 1;
         b.sign = 1;
 
         /* The result will be signed like the arg with greatest magnitude */
         byte rsign = bsign;
         if (compare(a, b) > 0) {
             rsign = asign;
         }
 
         /* Handle special case when a or b is zero, by setting the exponent
        of the zero number equal to the other one.  This avoids an alignment
        which would cause catastropic loss of precision */
         if (b.mant[mant.length-1] == 0) {
             b.exp = a.exp;
         }
 
         if (a.mant[mant.length-1] == 0) {
             a.exp = b.exp;
         }
 
         /* align number with the smaller exponent */
         int aextradigit = 0;
         int bextradigit = 0;
         if (a.exp < b.exp) {
             aextradigit = a.align(b.exp);
         } else {
             bextradigit = b.align(a.exp);
         }
 
         /* complement the smaller of the two if the signs are different */
         if (asign != bsign) {
             if (asign == rsign) {
                 bextradigit = b.complement(bextradigit);
             } else {
                 aextradigit = a.complement(aextradigit);
             }
         }
 
         /* add the mantissas */
         int rh = 0; /* acts as a carry */
         for (int i = 0; i < mant.length; i++) {
             final int r = a.mant[i]+b.mant[i]+rh;
             rh = r / RADIX;
             result.mant[i] = r - rh * RADIX;
         }
         result.exp = a.exp;
         result.sign = rsign;
 
         /* handle overflow -- note, when asign!=bsign an overflow is
          * normal and should be ignored.  */
 
         if (rh != 0 && (asign == bsign)) {
             final int lostdigit = result.mant[0];
             result.shiftRight();
             result.mant[mant.length-1] = rh;
             final int excp = result.round(lostdigit);
             if (excp != 0) {
                 result = dotrap(excp, ADD_TRAP, x, result);
             }
         }
 
         /* normalize the result */
         for (int i = 0; i < mant.length; i++) {
             if (result.mant[mant.length-1] != 0) {
                 break;
             }
             result.shiftLeft();
             if (i == 0) {
                 result.mant[0] = aextradigit+bextradigit;
                 aextradigit = 0;
                 bextradigit = 0;
             }
         }
 
         /* result is zero if after normalization the most sig. digit is zero */
         if (result.mant[mant.length-1] == 0) {
             result.exp = 0;
 
             if (asign != bsign) {
                 // Unless adding 2 negative zeros, sign is positive
                 result.sign = 1;  // Per IEEE 854-1987 Section 6.3
             }
         }
 
         /* Call round to test for over/under flows */
         final int excp = result.round(aextradigit + bextradigit);
         if (excp != 0) {
             result = dotrap(excp, ADD_TRAP, x, result);
         }
 
         return result;
     }
 
     /** Returns a number that is this number with the sign bit reversed.
      * @return the opposite of this
      */
     @Override
     public Dfp negate() {
         Dfp result = newInstance(this);
         result.sign = (byte) - result.sign;
         return result;
     }
 
     /** Subtract x from this.
      * @param x number to subtract
      * @return difference of this and a
      */
     @Override
     public Dfp subtract(final Dfp x) {
         return add(x.negate());
     }
 
     /** Round this given the next digit n using the current rounding mode.
      * @param n ???
      * @return the IEEE flag if an exception occurred
      */
     protected int round(int n) {
         boolean inc;
         switch (field.getRoundingMode()) {
             case ROUND_DOWN:
                 inc = false;
                 break;
 
             case ROUND_UP:
                 inc = n != 0;       // round up if n!=0
                 break;
 
             case ROUND_HALF_UP:
                 inc = n >= 5000;  // round half up
                 break;
 
             case ROUND_HALF_DOWN:
                 inc = n > 5000;  // round half down
                 break;
 
             case ROUND_HALF_EVEN:
                 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1);  // round half-even
                 break;
 
             case ROUND_HALF_ODD:
                 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0);  // round half-odd
                 break;
 
             case ROUND_CEIL:
                 inc = sign == 1 && n != 0;  // round ceil
                 break;
 
             case ROUND_FLOOR:
             default:
                 inc = sign == -1 && n != 0;  // round floor
                 break;
         }
 
         if (inc) {
             // increment if necessary
             int rh = 1;
             for (int i = 0; i < mant.length; i++) {
                 final int r = mant[i] + rh;
                 rh = r / RADIX;
                 mant[i] = r - rh * RADIX;
             }
 
             if (rh != 0) {
                 shiftRight();
                 mant[mant.length-1] = rh;
             }
         }
 
         // check for exceptional cases and raise signals if necessary
         if (exp < MIN_EXP) {
             // Gradual Underflow
             field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
             return DfpField.FLAG_UNDERFLOW;
         }
 
         if (exp > MAX_EXP) {
             // Overflow
             field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
             return DfpField.FLAG_OVERFLOW;
         }
 
         if (n != 0) {
             // Inexact
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
             return DfpField.FLAG_INEXACT;
         }
 
         return 0;
 
     }
 
     /** Multiply this by x.
      * @param x multiplicand
      * @return product of this and x
      */
     @Override
     public Dfp multiply(final Dfp x) {
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != x.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
         }
 
         Dfp result = newInstance(getZero());
 
         /* handle special cases */
         if (nans != FINITE || x.nans != FINITE) {
             if (isNaN()) {
                 return this;
             }
 
             if (x.isNaN()) {
                 return x;
             }
 
             if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] != 0) {
                 result = newInstance(this);
                 result.sign = (byte) (sign * x.sign);
                 return result;
             }
 
             if (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] != 0) {
                 result = newInstance(x);
                 result.sign = (byte) (sign * x.sign);
                 return result;
             }
 
             if (x.nans == INFINITE && nans == INFINITE) {
                 result = newInstance(this);
                 result.sign = (byte) (sign * x.sign);
                 return result;
             }
 
             if ( (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] == 0) ||
                     (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] == 0) ) {
                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
                 result = newInstance(getZero());
                 result.nans = QNAN;
                 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
                 return result;
             }
         }
 
         int[] product = new int[mant.length*2];  // Big enough to hold even the largest result
 
         for (int i = 0; i < mant.length; i++) {
             int rh = 0;  // acts as a carry
             for (int j=0; j<mant.length; j++) {
                 int r = mant[i] * x.mant[j];    // multiply the 2 digits
                 r += product[i+j] + rh;  // add to the product digit with carry in
 
                 rh = r / RADIX;
                 product[i+j] = r - rh * RADIX;
             }
             product[i+mant.length] = rh;
         }
 
         // Find the most sig digit
         int md = mant.length * 2 - 1;  // default, in case result is zero
         for (int i = mant.length * 2 - 1; i >= 0; i--) {
             if (product[i] != 0) {
                 md = i;
                 break;
             }
         }
 
         // Copy the digits into the result
         for (int i = 0; i < mant.length; i++) {
             result.mant[mant.length - i - 1] = product[md - i];
         }
 
         // Fixup the exponent.
         result.exp = exp + x.exp + md - 2 * mant.length + 1;
         result.sign = (byte)((sign == x.sign)?1:-1);
 
         if (result.mant[mant.length-1] == 0) {
             // if result is zero, set exp to zero
             result.exp = 0;
         }
 
         final int excp;
         if (md > (mant.length-1)) {
             excp = result.round(product[md-mant.length]);
         } else {
             excp = result.round(0); // has no effect except to check status
         }
 
         if (excp != 0) {
             result = dotrap(excp, MULTIPLY_TRAP, x, result);
         }
 
         return result;
 
     }
 
     /** Multiply this by a single digit x.
      * @param x multiplicand
      * @return product of this and x
      */
     @Override
     public Dfp multiply(final int x) {
         if (x >= 0 && x < RADIX) {
             return multiplyFast(x);
         } else {
             return multiply(newInstance(x));
         }
     }
 
     /** Multiply this by a single digit 0&lt;=x&lt;radix.
      * There are speed advantages in this special case.
      * @param x multiplicand
      * @return product of this and x
      */
     private Dfp multiplyFast(final int x) {
         Dfp result = newInstance(this);
 
         /* handle special cases */
         if (nans != FINITE) {
             if (isNaN()) {
                 return this;
             }
 
             if (nans == INFINITE && x != 0) {
                 result = newInstance(this);
                 return result;
             }
 
             if (nans == INFINITE && x == 0) {
                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
                 result = newInstance(getZero());
                 result.nans = QNAN;
                 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
                 return result;
             }
         }
 
         /* range check x */
         if (x < 0 || x >= RADIX) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             result = newInstance(getZero());
             result.nans = QNAN;
             result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
             return result;
         }
 
         int rh = 0;
         for (int i = 0; i < mant.length; i++) {
             final int r = mant[i] * x + rh;
             rh = r / RADIX;
             result.mant[i] = r - rh * RADIX;
         }
 
         int lostdigit = 0;
         if (rh != 0) {
             lostdigit = result.mant[0];
             result.shiftRight();
             result.mant[mant.length-1] = rh;
         }
 
         if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
             result.exp = 0;
         }
 
         final int excp = result.round(lostdigit);
         if (excp != 0) {
             result = dotrap(excp, MULTIPLY_TRAP, result, result);
         }
 
         return result;
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp square() {
         return multiply(this);
     }
 
     /** Divide this by divisor.
      * @param divisor divisor
      * @return quotient of this by divisor
      */
     @Override
     public Dfp divide(Dfp divisor) {
         int[] dividend; // current status of the dividend
         int[] quotient; // quotient
         int[] remainder;// remainder
         int qd;         // current quotient digit we're working with
         int nsqd;       // number of significant quotient digits we have
         int trial=0;    // trial quotient digit
         int minadj;     // minimum adjustment
         boolean trialgood; // Flag to indicate a good trail digit
         int md;         // most sig digit in result
         int excp;       // exceptions
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
         }
 
         Dfp result = newInstance(getZero());
 
         /* handle special cases */
         if (nans != FINITE || divisor.nans != FINITE) {
             if (isNaN()) {
                 return this;
             }
 
             if (divisor.isNaN()) {
                 return divisor;
             }
 
             if (nans == INFINITE && divisor.nans == FINITE) {
                 result = newInstance(this);
                 result.sign = (byte) (sign * divisor.sign);
                 return result;
             }
 
             if (divisor.nans == INFINITE && nans == FINITE) {
                 result = newInstance(getZero());
                 result.sign = (byte) (sign * divisor.sign);
                 return result;
             }
 
             if (divisor.nans == INFINITE && nans == INFINITE) {
                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
                 result = newInstance(getZero());
                 result.nans = QNAN;
                 result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
                 return result;
             }
         }
 
         /* Test for divide by zero */
         if (divisor.mant[mant.length-1] == 0) {
             field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
             result = newInstance(getZero());
             result.sign = (byte) (sign * divisor.sign);
             result.nans = INFINITE;
             result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
             return result;
         }
 
         dividend = new int[mant.length+1];  // one extra digit needed
         quotient = new int[mant.length+2];  // two extra digits needed 1 for overflow, 1 for rounding
         remainder = new int[mant.length+1]; // one extra digit needed
 
         /* Initialize our most significant digits to zero */
 
         dividend[mant.length] = 0;
         quotient[mant.length] = 0;
         quotient[mant.length+1] = 0;
         remainder[mant.length] = 0;
 
         /* copy our mantissa into the dividend, initialize the
        quotient while we are at it */
 
         for (int i = 0; i < mant.length; i++) {
             dividend[i] = mant[i];
             quotient[i] = 0;
             remainder[i] = 0;
         }
 
         /* outer loop.  Once per quotient digit */
         nsqd = 0;
         for (qd = mant.length+1; qd >= 0; qd--) {
             /* Determine outer limits of our quotient digit */
 
             // r =  most sig 2 digits of dividend
             final int divMsb = dividend[mant.length]*RADIX+dividend[mant.length-1];
             int min = divMsb       / (divisor.mant[mant.length-1]+1);
             int max = (divMsb + 1) / divisor.mant[mant.length-1];
 
             trialgood = false;
             while (!trialgood) {
                 // try the mean
                 trial = (min+max)/2;
 
                 /* Multiply by divisor and store as remainder */
                 int rh = 0;
                 for (int i = 0; i < mant.length + 1; i++) {
                     int dm = (i<mant.length)?divisor.mant[i]:0;
                     final int r = (dm * trial) + rh;
                     rh = r / RADIX;
                     remainder[i] = r - rh * RADIX;
                 }
 
                 /* subtract the remainder from the dividend */
                 rh = 1;  // carry in to aid the subtraction
                 for (int i = 0; i < mant.length + 1; i++) {
                     final int r = ((RADIX-1) - remainder[i]) + dividend[i] + rh;
                     rh = r / RADIX;
                     remainder[i] = r - rh * RADIX;
                 }
 
                 /* Lets analyze what we have here */
                 if (rh == 0) {
                     // trial is too big -- negative remainder
                     max = trial-1;
                     continue;
                 }
 
                 /* find out how far off the remainder is telling us we are */
                 minadj = (remainder[mant.length] * RADIX)+remainder[mant.length-1];
                 minadj /= divisor.mant[mant.length-1] + 1;
 
                 if (minadj >= 2) {
                     min = trial+minadj;  // update the minimum
                     continue;
                 }
 
                 /* May have a good one here, check more thoroughly.  Basically
            its a good one if it is less than the divisor */
                 trialgood = false;  // assume false
                 for (int i = mant.length - 1; i >= 0; i--) {
                     if (divisor.mant[i] > remainder[i]) {
                         trialgood = true;
                     }
                     if (divisor.mant[i] < remainder[i]) {
                         break;
                     }
                 }
 
                 if (remainder[mant.length] != 0) {
                     trialgood = false;
                 }
 
                 if (!trialgood) {
                     min = trial+1;
                 }
             }
 
             /* Great we have a digit! */
             quotient[qd] = trial;
             if (trial != 0 || nsqd != 0) {
                 nsqd++;
             }
 
             if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
                 // We have enough for this mode
                 break;
             }
 
             if (nsqd > mant.length) {
                 // We have enough digits
                 break;
             }
 
             /* move the remainder into the dividend while left shifting */
             dividend[0] = 0;
             System.arraycopy(remainder, 0, dividend, 1, mant.length);
         }
 
         /* Find the most sig digit */
         md = mant.length;  // default
         for (int i = mant.length + 1; i >= 0; i--) {
             if (quotient[i] != 0) {
                 md = i;
                 break;
             }
         }
 
         /* Copy the digits into the result */
         for (int i=0; i<mant.length; i++) {
             result.mant[mant.length-i-1] = quotient[md-i];
         }
 
         /* Fixup the exponent. */
         result.exp = exp - divisor.exp + md - mant.length;
         result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
 
         if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
             result.exp = 0;
         }
 
         if (md > (mant.length-1)) {
             excp = result.round(quotient[md-mant.length]);
         } else {
             excp = result.round(0);
         }
 
         if (excp != 0) {
             result = dotrap(excp, DIVIDE_TRAP, divisor, result);
         }
 
         return result;
     }
 
     /** Divide by a single digit less than radix.
      *  Special case, so there are speed advantages. 0 &lt;= divisor &lt; radix
      * @param divisor divisor
      * @return quotient of this by divisor
      */
     public Dfp divide(int divisor) {
 
         // Handle special cases
         if (nans != FINITE) {
             if (isNaN()) {
                 return this;
             }
 
             if (nans == INFINITE) {
                 return newInstance(this);
             }
         }
 
         // Test for divide by zero
         if (divisor == 0) {
             field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
             Dfp result = newInstance(getZero());
             result.sign = sign;
             result.nans = INFINITE;
             result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
             return result;
         }
 
         // range check divisor
         if (divisor < 0 || divisor >= RADIX) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             Dfp result = newInstance(getZero());
             result.nans = QNAN;
             result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
             return result;
         }
 
         Dfp result = newInstance(this);
 
         int rl = 0;
         for (int i = mant.length-1; i >= 0; i--) {
             final int r = rl*RADIX + result.mant[i];
             final int rh = r / divisor;
             rl = r - rh * divisor;
             result.mant[i] = rh;
         }
 
         if (result.mant[mant.length-1] == 0) {
             // normalize
             result.shiftLeft();
             final int r = rl * RADIX;        // compute the next digit and put it in
             final int rh = r / divisor;
             rl = r - rh * divisor;
             result.mant[0] = rh;
         }
 
         final int excp = result.round(rl * RADIX / divisor);  // do the rounding
         if (excp != 0) {
             result = dotrap(excp, DIVIDE_TRAP, result, result);
         }
 
         return result;
 
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp reciprocal() {
         return field.getOne().divide(this);
     }
 
     /** Compute the square root.
      * @return square root of the instance
      */
     @Override
     public Dfp sqrt() {
 
         // check for unusual cases
         if (nans == FINITE && mant[mant.length-1] == 0) {
             // if zero
             return newInstance(this);
         }
 
         if (nans != FINITE) {
             if (nans == INFINITE && sign == 1) {
                 // if positive infinity
                 return newInstance(this);
             }
 
             if (nans == QNAN) {
                 return newInstance(this);
             }
 
             if (nans == SNAN) {
                 Dfp result;
 
                 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
                 result = newInstance(this);
                 result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
                 return result;
             }
         }
 
         if (sign == -1) {
             // if negative
             Dfp result;
 
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             result = newInstance(this);
             result.nans = QNAN;
             result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
             return result;
         }
 
         Dfp x = newInstance(this);
 
         /* Lets make a reasonable guess as to the size of the square root */
         if (x.exp < -1 || x.exp > 1) {
             x.exp = this.exp / 2;
         }
 
         /* Coarsely estimate the mantissa */
         switch (x.mant[mant.length-1] / 2000) {
             case 0:
                 x.mant[mant.length-1] = x.mant[mant.length-1]/2+1;
                 break;
             case 2:
                 x.mant[mant.length-1] = 1500;
                 break;
             case 3:
                 x.mant[mant.length-1] = 2200;
                 break;
             default:
                 x.mant[mant.length-1] = 3000;
                 break;
         }
 
         /* Now that we have the first pass estimate, compute the rest
        by the formula dx = (y - x*x) / (2x); */
 
         Dfp dx;
         Dfp px  = getZero();
         Dfp ppx;
         while (x.unequal(px)) {
             dx = newInstance(x);
             dx.sign = -1;
             dx = dx.add(this.divide(x));
             dx = dx.divide(2);
             ppx = px;
             px = x;
             x = x.add(dx);
 
             if (x.equals(ppx)) {
                 // alternating between two values
                 break;
             }
 
             // if dx is zero, break.  Note testing the most sig digit
             // is a sufficient test since dx is normalized
             if (dx.mant[mant.length-1] == 0) {
                 break;
             }
         }
 
         return x;
 
     }
 
     /** Get a string representation of the instance.
      * @return string representation of the instance
      */
     @Override
     public String toString() {
         if (nans != FINITE) {
             // if non-finite exceptional cases
             if (nans == INFINITE) {
                 return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
             } else {
                 return NAN_STRING;
             }
         }
 
         if (exp > mant.length || exp < -1) {
             return dfp2sci();
         }
 
         return dfp2string();
 
     }
 
     /** Convert an instance to a string using scientific notation.
      * @return string representation of the instance in scientific notation
      */
     protected String dfp2sci() {
         char[] rawdigits = new char[mant.length * 4];
         int p;
         int e;
         int ae;
         int shf;
 
         // Get all the digits
         p = 0;
         for (int i = mant.length - 1; i >= 0; i--) {
             rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
             rawdigits[p++] = (char) (((mant[i] / 100) %10) + '0');
             rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
             rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
         }
 
         // Find the first non-zero one
         for (p = 0; p < rawdigits.length; p++) {
             if (rawdigits[p] != '0') {
                 break;
             }
         }
         shf = p;
 
         // Now do the conversion
         StringBuilder builder = new StringBuilder();
         if (sign == -1) {
             builder.append('-');
         }
 
         if (p != rawdigits.length) {
             // there are non zero digits...
             builder.append(rawdigits[p++]).append('.');
 
             while (p<rawdigits.length) {
                 builder.append(rawdigits[p++]);
             }
         } else {
             builder.append("0.0e0");
             return builder.toString();
         }
 
         builder.append('e');
 
         // Find the msd of the exponent
 
         e = exp * 4 - shf - 1;
         ae = e;
         if (e < 0) {
             ae = -e;
         }
 
         // Find the largest p such that p < e
         for (p = 1000000000; p > ae; p /= 10) { // NOPMD - empty loop is normal here
             // nothing to do
         }
 
         if (e < 0) {
             builder.append('-');
         }
 
         while (p > 0) {
             builder.append((char)(ae / p + '0'));
             ae %= p;
             p /= 10;
         }
 
         return builder.toString();
 
     }
 
     /** Convert an instance to a string using normal notation.
      * @return string representation of the instance in normal notation
      */
     protected String dfp2string() {
         final String fourZero = "0000";
         int e = exp;
         boolean pointInserted = false;
 
         StringBuilder builder = new StringBuilder(23);
 
         if (e <= 0) {
             builder.append("0.");
             pointInserted = true;
         }
 
         while (e < 0) {
             builder.append(fourZero);
             e++;
         }
 
         for (int i = mant.length - 1; i >= 0; i--) {
             builder.
                 append((char) ((mant[i] / 1000) + '0')).
                 append((char) (((mant[i] / 100) % 10) + '0')).
                 append((char) (((mant[i] / 10) % 10) + '0')).
                 append((char) (((mant[i]) % 10) + '0'));
             --e;
             if (e == 0) {
                 builder.append('.');
                 pointInserted = true;
             }
         }
 
         while (e > 0) {
             builder.append(fourZero);
             e--;
         }
 
         if (!pointInserted) {
             // Ensure we have a radix point!
             builder.append('.');
         }
 
         // Suppress leading zeros
         while (builder.charAt(0) == '0') {
             builder.deleteCharAt(0);
         }
         if (builder.charAt(0) == '.') {
             builder.insert(0, '0');
         }
 
         // Suppress trailing zeros
         while (builder.charAt(builder.length() - 1) == '0') {
             builder.deleteCharAt(builder.length() - 1);
         }
 
         // Insert sign
         if (sign < 0) {
             builder.insert(0, '-');
         }
 
         return builder.toString();
 
     }
 
     /** Raises a trap.  This does not set the corresponding flag however.
      *  @param type the trap type
      *  @param what - name of routine trap occurred in
      *  @param oper - input operator to function
      *  @param result - the result computed prior to the trap
      *  @return The suggested return value from the trap handler
      */
     public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
         Dfp def = result;
 
         switch (type) {
             case DfpField.FLAG_INVALID:
                 def = newInstance(getZero());
                 def.sign = result.sign;
                 def.nans = QNAN;
                 break;
 
             case DfpField.FLAG_DIV_ZERO:
                 if (nans == FINITE && mant[mant.length-1] != 0) {
                     // normal case, we are finite, non-zero
                     def = newInstance(getZero());
                     def.sign = (byte)(sign*oper.sign);
                     def.nans = INFINITE;
                 }
 
                 if (nans == FINITE && mant[mant.length-1] == 0) {
                     //  0/0
                     def = newInstance(getZero());
                     def.nans = QNAN;
                 }
 
                 if (nans == INFINITE || nans == QNAN) {
                     def = newInstance(getZero());
                     def.nans = QNAN;
                 }
 
                 if (nans == INFINITE || nans == SNAN) {
                     def = newInstance(getZero());
                     def.nans = QNAN;
                 }
                 break;
 
             case DfpField.FLAG_UNDERFLOW:
                 if ( (result.exp+mant.length) < MIN_EXP) {
                     def = newInstance(getZero());
                     def.sign = result.sign;
                 } else {
                     def = newInstance(result);  // gradual underflow
                 }
                 result.exp += ERR_SCALE;
                 break;
 
             case DfpField.FLAG_OVERFLOW:
                 result.exp -= ERR_SCALE;
                 def = newInstance(getZero());
                 def.sign = result.sign;
                 def.nans = INFINITE;
                 break;
 
             default: def = result; break;
         }
 
         return trap(type, what, oper, def, result);
 
     }
 
     /** Trap handler.  Subclasses may override this to provide trap
      *  functionality per IEEE 854-1987.
      *
      *  @param type  The exception type - e.g. FLAG_OVERFLOW
      *  @param what  The name of the routine we were in e.g. divide()
      *  @param oper  An operand to this function if any
      *  @param def   The default return value if trap not enabled
      *  @param result    The result that is specified to be delivered per
      *                   IEEE 854, if any
      *  @return the value that should be return by the operation triggering the trap
      */
     protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
         return def;
     }
 
     /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN.
      * @return type of the number
      */
     public int classify() {
         return nans;
     }
 
     /** Creates an instance that is the same as x except that it has the sign of y.
      * abs(x) = dfp.copysign(x, dfp.one)
      * @param x number to get the value from
      * @param y number to get the sign from
      * @return a number with the value of x and the sign of y
      */
     public static Dfp copysign(final Dfp x, final Dfp y) {
         Dfp result = x.newInstance(x);
         result.sign = y.sign;
         return result;
     }
 
     /** Returns the next number greater than this one in the direction of x.
      * If this==x then simply returns this.
      * @param x direction where to look at
      * @return closest number next to instance in the direction of x
      */
     public Dfp nextAfter(final Dfp x) {
 
         // make sure we don't mix number with different precision
         if (field.getRadixDigits() != x.field.getRadixDigits()) {
             field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
             final Dfp result = newInstance(getZero());
             result.nans = QNAN;
             return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
         }
 
         // if this is greater than x
         boolean up = false;
         if (this.lessThan(x)) {
             up = true;
         }
 
         if (compare(this, x) == 0) {
             return newInstance(x);
         }
 
         if (lessThan(getZero())) {
             up = !up;
         }
 
         final Dfp inc;
         Dfp result;
         if (up) {
             inc = newInstance(getOne());
             inc.exp = this.exp-mant.length+1;
             inc.sign = this.sign;
 
             if (this.equals(getZero())) {
                 inc.exp = MIN_EXP-mant.length;
             }
 
             result = add(inc);
         } else {
             inc = newInstance(getOne());
             inc.exp = this.exp;
             inc.sign = this.sign;
 
             if (this.equals(inc)) {
                 inc.exp = this.exp-mant.length;
             } else {
                 inc.exp = this.exp-mant.length+1;
             }
 
             if (this.equals(getZero())) {
                 inc.exp = MIN_EXP-mant.length;
             }
 
             result = this.subtract(inc);
         }
 
         if (result.classify() == INFINITE && this.classify() != INFINITE) {
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
             result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
         }
 
         if (result.equals(getZero()) && !this.equals(getZero())) {
             field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
             result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
         }
 
         return result;
 
     }
 
     /** Convert the instance into a double.
      * @return a double approximating the instance
      * @see #toSplitDouble()
      */
     public double toDouble() {
 
         if (isInfinite()) {
             if (lessThan(getZero())) {
                 return Double.NEGATIVE_INFINITY;
             } else {
                 return Double.POSITIVE_INFINITY;
             }
         }
 
         if (isNaN()) {
             return Double.NaN;
         }
 
         Dfp y = this;
         boolean negate = false;
         int cmp0 = compare(this, getZero());
         if (cmp0 == 0) {
             return sign < 0 ? -0.0 : +0.0;
         } else if (cmp0 < 0) {
             y = negate();
             negate = true;
         }
 
         /* Find the exponent, first estimate by integer log10, then adjust.
          Should be faster than doing a natural logarithm.  */
         int exponent = (int)(y.intLog10() * 3.32);
         if (exponent < 0) {
             exponent--;
         }
 
         Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
         while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
             tempDfp = tempDfp.multiply(2);
             exponent++;
         }
         exponent--;
 
         /* We have the exponent, now work on the mantissa */
 
         y = y.divide(DfpMath.pow(getTwo(), exponent));
         if (exponent > -1023) {
             y = y.subtract(getOne());
         }
 
         if (exponent < -1074) {
             return 0;
         }
 
         if (exponent > 1023) {
             return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
         }
 
 
         y = y.multiply(newInstance(4503599627370496L)).rint();
         String str = y.toString();
         str = str.substring(0, str.length()-1);
         long mantissa = Long.parseLong(str);
 
         if (mantissa == 4503599627370496L) {
             // Handle special case where we round up to next power of two
             mantissa = 0;
             exponent++;
         }
 
         /* Its going to be subnormal, so make adjustments */
         if (exponent <= -1023) {
             exponent--;
         }
 
         while (exponent < -1023) {
             exponent++;
             mantissa >>>= 1;
         }
 
         long bits = mantissa | ((exponent + 1023L) << 52);
         double x = Double.longBitsToDouble(bits);
 
         if (negate) {
             x = -x;
         }
 
         return x;
 
     }
 
     /** Convert the instance into a split double.
      * @return an array of two doubles which sum represent the instance
      * @see #toDouble()
      */
     public double[] toSplitDouble() {
         double[] split = new double[2];
         long mask = 0xffffffffc0000000L;
 
         split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
         split[1] = subtract(newInstance(split[0])).toDouble();
 
         return split;
     }
 
     /** {@inheritDoc}
      */
     @Override
     public double getReal() {
         return toDouble();
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp getAddendum() {
         return isFinite() ? subtract(getReal()) : getZero();
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp remainder(final double a) {
         return remainder(newInstance(a));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp sign() {
         if (isNaN() || isZero()) {
             return this;
         } else {
             return newInstance(sign > 0 ? +1 : -1);
         }
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp copySign(final Dfp s) {
         if ((sign >= 0 && s.sign >= 0) || (sign < 0 && s.sign < 0)) { // Sign is currently OK
             return this;
         }
         return negate(); // flip sign
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp copySign(final double s) {
         long sb = Double.doubleToLongBits(s);
         if ((sign >= 0 && sb >= 0) || (sign < 0 && sb < 0)) { // Sign is currently OK
             return this;
         }
         return negate(); // flip sign
     }
 
     /** {@inheritDoc}
      */
     @Override
     public int getExponent() {
 
         if (nans != FINITE) {
             // 2⁴³⁵⁴¹¹ < 10000³²⁷⁶⁸ < 2⁴³⁵⁴¹²
             return 435411;
         }
         if (isZero()) {
             return -435412;
         }
 
         final Dfp abs = abs();
 
         // estimate a lower bound for binary exponent
         // 13301/1001 is a continued fraction approximation of ln(10000)/ln(2)
         int p = FastMath.max(13301 * exp / 1001 - 15, -435411);
         Dfp twoP = DfpMath.pow(getTwo(), p);
         while (compare(abs, twoP) >= 0) {
             twoP = twoP.add(twoP);
             ++p;
         }
 
         return p - 1;
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp scalb(final int n) {
         return multiply(DfpMath.pow(getTwo(), n));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp ulp() {
         final Dfp result = new Dfp(field);
         result.mant[result.mant.length - 1] = 1;
         result.exp = exp - (result.mant.length - 1);
         return result;
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp hypot(final Dfp y) {
 
         if (isInfinite() || y.isInfinite()) {
             return field.newDfp(Double.POSITIVE_INFINITY);
         } else if (isNaN() || y.isNaN()) {
             return field.newDfp(Double.NaN);
         } else {
             // find scaling to avoid both overflow and underflow
             final int scalingExp = (exp + y.exp) / 2;
 
             // scale both operands
             final Dfp scaledX = new Dfp(this);
             scaledX.exp -= scalingExp;
             final Dfp scaledY = new Dfp(y);
             scaledY.exp -= scalingExp;
 
             // compute scaled hypothenuse
             final Dfp h = scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
 
             // scale result
             h.exp += scalingExp;
 
             return h;
         }
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp rootN(final int n) {
         return (sign >= 0) ?
                DfpMath.pow(this, getOne().divide(n)) :
                DfpMath.pow(negate(), getOne().divide(n)).negate();
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp pow(final double p) {
         return DfpMath.pow(this, newInstance(p));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp pow(final int n) {
         return DfpMath.pow(this, n);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp pow(final Dfp e) {
         return DfpMath.pow(this, e);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp exp() {
         return DfpMath.exp(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp expm1() {
         return DfpMath.exp(this).subtract(getOne());
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp log() {
         return DfpMath.log(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp log1p() {
         return DfpMath.log(this.add(getOne()));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp log10() {
         return DfpMath.log(this).divide(DfpMath.log(newInstance(10)));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp cos() {
         return DfpMath.cos(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp sin() {
         return DfpMath.sin(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp tan() {
         return DfpMath.tan(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp acos() {
         return DfpMath.acos(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp asin() {
         return DfpMath.asin(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp atan() {
         return DfpMath.atan(this);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp atan2(final Dfp x)
         throws MathIllegalArgumentException {
 
         // compute r = sqrt(x^2+y^2)
         final Dfp r = x.square().add(multiply(this)).sqrt();
         if (r.isZero()) {
             // special cases handling
             if (x.sign >= 0) {
                 return this; // ±0.0
             } else {
                 return newInstance((sign <= 0) ? -FastMath.PI : FastMath.PI); // ±π
             }
         }
 
         if (x.sign >= 0) {
 
             // compute atan2(y, x) = 2 atan(y / (r + x))
             return getTwo().multiply(divide(r.add(x)).atan());
 
         } else {
 
             // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
             final Dfp tmp = getTwo().multiply(divide(r.subtract(x)).atan());
             final Dfp pmPi = newInstance((tmp.sign <= 0) ? -FastMath.PI : FastMath.PI);
             return pmPi.subtract(tmp);
 
         }
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp cosh() {
         return DfpMath.exp(this).add(DfpMath.exp(negate())).multiply(0.5);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp sinh() {
         return DfpMath.exp(this).subtract(DfpMath.exp(negate())).multiply(0.5);
     }
 
     /** {@inheritDoc}
      */
     @Override
     public FieldSinhCosh<Dfp> sinhCosh() {
         final Dfp p = DfpMath.exp(this);
         final Dfp m = DfpMath.exp(negate());
         return new FieldSinhCosh<>(p.subtract(m).multiply(0.5), p.add(m).multiply(0.5));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp tanh() {
         final Dfp ePlus  = DfpMath.exp(this);
         final Dfp eMinus = DfpMath.exp(negate());
         return ePlus.subtract(eMinus).divide(ePlus.add(eMinus));
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp acosh() {
         return square().subtract(getOne()).sqrt().add(this).log();
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp asinh() {
         return square().add(getOne()).sqrt().add(this).log();
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp atanh() {
         return getOne().add(this).divide(getOne().subtract(this)).log().divide(2);
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp toDegrees() {
         return multiply(field.getRadToDeg());
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp toRadians() {
         return multiply(field.getDegToRad());
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final Dfp[] a, final Dfp[] b)
         throws MathIllegalArgumentException {
 
         MathUtils.checkDimension(a.length, b.length);
 
         // compute in extended accuracy
         final DfpField extendedField = a[0].field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         Dfp r = extendedField.getZero();
         for (int i = 0; i < a.length; ++i) {
             final Dfp aiExt = a[i].newInstance(extendedField, null);
             final Dfp biExt = b[i].newInstance(extendedField, null);
             r = r.add(aiExt.multiply(biExt));
         }
 
         // back to normal accuracy
         return r.newInstance(a[0].field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final double[] a, final Dfp[] b)
         throws MathIllegalArgumentException {
 
         MathUtils.checkDimension(a.length, b.length);
 
         // compute in extended accuracy
         final DfpField extendedField = b[0].field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         Dfp r = extendedField.getZero();
         for (int i = 0; i < a.length; ++i) {
             final Dfp biExt = b[i].newInstance(extendedField, null);
             r = r.add(biExt.multiply(a[i]));
         }
 
         // back to normal accuracy
         return r.newInstance(b[0].field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2) {
 
         // switch to extended accuracy
         final DfpField extendedField = a1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         final Dfp a1Ext = a1.newInstance(extendedField, null);
         final Dfp b1Ext = b1.newInstance(extendedField, null);
         final Dfp a2Ext = a2.newInstance(extendedField, null);
         final Dfp b2Ext = b2.newInstance(extendedField, null);
 
         // compute linear combination in extended accuracy
         final Dfp resultExt = a1Ext.multiply(b1Ext).
                           add(a2Ext.multiply(b2Ext));
 
         // back to normal accuracy
         return resultExt.newInstance(a1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2) {
 
         // switch to extended accuracy
         final DfpField extendedField = b1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         final Dfp b1Ext = b1.newInstance(extendedField, null);
         final Dfp b2Ext = b2.newInstance(extendedField, null);
 
         // compute linear combination in extended accuracy
         final Dfp resultExt = b1Ext.multiply(a1).
                           add(b2Ext.multiply(a2));
 
         // back to normal accuracy
         return resultExt.newInstance(b1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final Dfp a1, final Dfp b1,
                                  final Dfp a2, final Dfp b2,
                                  final Dfp a3, final Dfp b3) {
 
         // switch to extended accuracy
         final DfpField extendedField = a1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         final Dfp a1Ext = a1.newInstance(extendedField, null);
         final Dfp b1Ext = b1.newInstance(extendedField, null);
         final Dfp a2Ext = a2.newInstance(extendedField, null);
         final Dfp b2Ext = b2.newInstance(extendedField, null);
         final Dfp a3Ext = a3.newInstance(extendedField, null);
         final Dfp b3Ext = b3.newInstance(extendedField, null);
 
         // compute linear combination in extended accuracy
         final Dfp resultExt = a1Ext.multiply(b1Ext).
                           add(a2Ext.multiply(b2Ext)).
                           add(a3Ext.multiply(b3Ext));
 
         // back to normal accuracy
         return resultExt.newInstance(a1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final double a1, final Dfp b1,
                                  final double a2, final Dfp b2,
                                  final double a3, final Dfp b3) {
 
         // switch to extended accuracy
         final DfpField extendedField = b1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         final Dfp b1Ext = b1.newInstance(extendedField, null);
         final Dfp b2Ext = b2.newInstance(extendedField, null);
         final Dfp b3Ext = b3.newInstance(extendedField, null);
 
         // compute linear combination in extended accuracy
         final Dfp resultExt = b1Ext.multiply(a1).
                           add(b2Ext.multiply(a2)).
                           add(b3Ext.multiply(a3));
 
         // back to normal accuracy
         return resultExt.newInstance(b1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2,
                                  final Dfp a3, final Dfp b3, final Dfp a4, final Dfp b4) {
 
         // switch to extended accuracy
         final DfpField extendedField = a1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         final Dfp a1Ext = a1.newInstance(extendedField, null);
         final Dfp b1Ext = b1.newInstance(extendedField, null);
         final Dfp a2Ext = a2.newInstance(extendedField, null);
         final Dfp b2Ext = b2.newInstance(extendedField, null);
         final Dfp a3Ext = a3.newInstance(extendedField, null);
         final Dfp b3Ext = b3.newInstance(extendedField, null);
         final Dfp a4Ext = a4.newInstance(extendedField, null);
         final Dfp b4Ext = b4.newInstance(extendedField, null);
 
         // compute linear combination in extended accuracy
         final Dfp resultExt = a1Ext.multiply(b1Ext).
                           add(a2Ext.multiply(b2Ext)).
                           add(a3Ext.multiply(b3Ext)).
                           add(a4Ext.multiply(b4Ext));
 
         // back to normal accuracy
         return resultExt.newInstance(a1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc}
      */
     @Override
     public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2,
                                  final double a3, final Dfp b3, final double a4, final Dfp b4) {
 
         // switch to extended accuracy
         final DfpField extendedField = b1.field.getExtendedField(LINEAR_COMBINATION_DIGITS_FACTOR, false);
         final Dfp b1Ext = b1.newInstance(extendedField, null);
         final Dfp b2Ext = b2.newInstance(extendedField, null);
         final Dfp b3Ext = b3.newInstance(extendedField, null);
         final Dfp b4Ext = b4.newInstance(extendedField, null);
 
         // compute linear combination in extended accuracy
         final Dfp resultExt = b1Ext.multiply(a1).
                           add(b2Ext.multiply(a2)).
                           add(b3Ext.multiply(a3)).
                           add(b4Ext.multiply(a4));
 
         // back to normal accuracy
         return resultExt.newInstance(b1.field, DfpField.RoundingMode.ROUND_HALF_EVEN);
 
     }
 
     /** {@inheritDoc} */
     @Override
     public Dfp getPi() {
         return field.getPi();
     }
 
 }