src/main/java/org/apache/commons/math/dfp/Dfp.java

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

package org.apache.commons.math.dfp;

import java.util.Arrays;

import org.apache.commons.math.FieldElement;

/**
 *  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:
 *  <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 keep 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>
 *
 *  <p>Trade offs:
 *  <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>
 *
 *  <p>Numbers are represented  in the following form:
 *  <pre>
 *  n  =  sign &times; mant &times; (radix)<sup>exp</sup>;</p>
 *  </pre>
 *  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 add 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
 * @version $Revision: 1003889 $ $Date: 2010-10-02 23:11:55 +0200 (sam. 02 oct. 2010) $
 * @since 2.2
 */
public class Dfp implements FieldElement<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";

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

    /** Sign bit: & 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()];
        sign = 1;
        exp = 0;
        nans = FINITE;
        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) {
                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;
        exp = 0;
        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 (s.equals(POS_INFINITY_STRING)) {
            sign = (byte) 1;
            nans = INFINITE;
            return;
        }

        if (s.equals(NEG_INFINITY_STRING)) {
            sign = (byte) -1;
            nans = INFINITE;
            return;
        }

        if (s.equals(NAN_STRING)) {
            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;

        int trailingZeros = 0;
        while (p > q) {
            if (striped[p] != '0') {
                break;
            }
            trailingZeros++;
            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);
    }

    /** Create an instance from a double value.
     * @param x value to convert to an instance
     * @return a new instance with value x
     */
    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);
    }

    /** Get the {@link org.apache.commons.math.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.apache.commons.math.Field Field} (really a {@link DfpField}) to which the instance belongs
     */
    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 infinite.
     * @return true if instance is infinite
     */
    public boolean isInfinite() {
        return nans == INFINITE;
    }

    /** Check if instance is not a number.
     * @return true if instance is not a number
     */
    public boolean isNaN() {
        return (nans == QNAN) || (nans == SNAN);
    }

    /** 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 + (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
     */
    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
     */
    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
     */
    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
     */
    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 zero */
        if (exp < 0) {
            field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
            Dfp result = newInstance(getZero());
            result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
            return 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 log10()  {
        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);
        }

        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 = 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
     */
    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
     */
    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
     */
    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 = false;
        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
     */
    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 = 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 0&lt;=x&lt;radix.
     * There are speed advantages in this special case
     * @param x multiplicand
     * @return product of this and x
     */
    public Dfp multiply(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;
    }

    /** Divide this by divisor.
     * @param divisor divisor
     * @return quotient of this by divisor
     */
    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=0;       // 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 = 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 == false) {
                    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;
            for (int i = 0; i < mant.length; i++) {
                dividend[i + 1] = remainder[i];
            }
        }

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

    }

    /** Compute the square root.
     * @return square root of the instance
     */
    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;
        }

        Dfp dx = newInstance(x);

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

        Dfp px  = getZero();
        Dfp ppx = getZero();
        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];
        char outputbuffer[] = new char[mant.length * 4 + 20];
        int p;
        int q;
        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
        q = 0;
        if (sign == -1) {
            outputbuffer[q++] = '-';
        }

        if (p != rawdigits.length) {
            // there are non zero digits...
            outputbuffer[q++] = rawdigits[p++];
            outputbuffer[q++] = '.';

            while (p<rawdigits.length) {
                outputbuffer[q++] = rawdigits[p++];
            }
        } else {
            outputbuffer[q++] = '0';
            outputbuffer[q++] = '.';
            outputbuffer[q++] = '0';
            outputbuffer[q++] = 'e';
            outputbuffer[q++] = '0';
            return new String(outputbuffer, 0, 5);
        }

        outputbuffer[q++] = '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) {
            // nothing to do
        }

        if (e < 0) {
            outputbuffer[q++] = '-';
        }

        while (p > 0) {
            outputbuffer[q++] = (char)(ae / p + '0');
            ae = ae % p;
            p = p / 10;
        }

        return new String(outputbuffer, 0, q);

    }

    /** Convert an instance to a string using normal notation.
     * @return string representation of the instance in normal notation
     */
    protected String dfp2string() {
        char buffer[] = new char[mant.length*4 + 20];
        int p = 1;
        int q;
        int e = exp;
        boolean pointInserted = false;

        buffer[0] = ' ';

        if (e <= 0) {
            buffer[p++] = '0';
            buffer[p++] = '.';
            pointInserted = true;
        }

        while (e < 0) {
            buffer[p++] = '0';
            buffer[p++] = '0';
            buffer[p++] = '0';
            buffer[p++] = '0';
            e++;
        }

        for (int i = mant.length - 1; i >= 0; i--) {
            buffer[p++] = (char) ((mant[i] / 1000) + '0');
            buffer[p++] = (char) (((mant[i] / 100) % 10) + '0');
            buffer[p++] = (char) (((mant[i] / 10) % 10) + '0');
            buffer[p++] = (char) (((mant[i]) % 10) + '0');
            if (--e == 0) {
                buffer[p++] = '.';
                pointInserted = true;
            }
        }

        while (e > 0) {
            buffer[p++] = '0';
            buffer[p++] = '0';
            buffer[p++] = '0';
            buffer[p++] = '0';
            e--;
        }

        if (!pointInserted) {
            // Ensure we have a radix point!
            buffer[p++] = '.';
        }

        // Suppress leading zeros
        q = 1;
        while (buffer[q] == '0') {
            q++;
        }
        if (buffer[q] == '.') {
            q--;
        }

        // Suppress trailing zeros
        while (buffer[p-1] == '0') {
            p--;
        }

        // Insert sign
        if (sign < 0) {
            buffer[--q] = '-';
        }

        return new String(buffer, q, p - q);

    }

    /** 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 = result.exp + ERR_SCALE;
                break;

            case DfpField.FLAG_OVERFLOW:
                result.exp = 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()) == false) {
            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;
        if (lessThan(getZero())) {
            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.log10() * 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;
    }

}