From bc7a1145b59255a7f7566a27a63e352517a8f9df Mon Sep 17 00:00:00 2001 From: tdurieux Date: Tue, 7 Mar 2017 13:17:08 +0100 Subject: [PATCH] buggy files form Math #54 --- .../54/org/apache/commons/math/dfp/Dfp.java | 2400 +++++++++++++++++ 1 file changed, 2400 insertions(+) create mode 100644 projects/Math/54/org/apache/commons/math/dfp/Dfp.java diff --git a/projects/Math/54/org/apache/commons/math/dfp/Dfp.java b/projects/Math/54/org/apache/commons/math/dfp/Dfp.java new file mode 100644 index 0000000..73f1ea3 --- /dev/null +++ b/projects/Math/54/org/apache/commons/math/dfp/Dfp.java @@ -0,0 +1,2400 @@ +/* + * 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 + * + *

Another floating point class. This one is built using radix 10000 + * which is 104, so its almost decimal.

+ * + *

The design goals here are: + *

    + *
  1. Decimal math, or close to it
    2. + *
  2. Settable precision (but no mix between numbers using different settings)
    3. + *
  3. Portability. Code should be keep as portable as possible.
    4. + *
  4. Performance
    5. + *
  5. Accuracy - Results should always be +/- 1 ULP for basic + * algebraic operation
    6. + *
  6. Comply with IEEE 854-1987 as much as possible. + * (See IEEE 854-1987 notes below)
    7. + *

+ * + *

Trade offs: + *

    + *
  1. Memory foot print. I'm using more memory than necessary to + * represent numbers to get better performance.
    2. + *
  2. 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.
    3. + *

+ * + *

Numbers are represented in the following form: + * 

+ *  n  =  sign × mant × (radix)exp;

+ *
+ * where sign is ±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

+ * + *

IEEE 854-1987 Notes and differences

+ * + *

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.

+ * + *

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.

+ * + *

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.

+ * + *

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.

+ * + *

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.

+ * @see DfpField + * @version $Revision$ $Date$ + * @since 2.2 + */ +public class Dfp implements FieldElement { + + /** 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: 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()]; + 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) { + // make sure 0 has the right sign + 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= '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. + *

+ * The field is linked to the number of digits and acts as a factory + * for {@link Dfp} instances. + *

+ * @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 ab 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 × 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 106 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 10000e + */ + 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 10e + */ + 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= 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<=x<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= 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-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 <= divisor < 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 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; + } + +}