Java tutorial
/* * 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.math3.dfp; import java.util.Arrays; import org.apache.commons.math3.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 kept as portable as possible.</li> * <li>Performance</li> * <li>Accuracy - Results should always be +/- 1 ULP for basic * algebraic operation</li> * <li>Comply with IEEE 854-1987 as much as possible. * (See IEEE 854-1987 notes below)</li> * </ol></p> * * <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 × mant × (radix)<sup>exp</sup>;</p> * </pre> * 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</p> * * <p>IEEE 854-1987 Notes and differences</p> * * <p>IEEE 854 requires the radix to be either 2 or 10. The radix here is * 10000, so that requirement is not met, but it is possible that a * subclassed can be made to make it behave as a radix 10 * number. It is my opinion that if it looks and behaves as a radix * 10 number then it is one and that requirement would be met.</p> * * <p>The radix of 10000 was chosen because it should be faster to operate * on 4 decimal digits at once instead of one at a time. Radix 10 behavior * can be realized by adding an additional rounding step to ensure that * the number of decimal digits represented is constant.</p> * * <p>The IEEE standard specifically leaves out internal data encoding, * so it is reasonable to conclude that such a subclass of this radix * 10000 system is merely an encoding of a radix 10 system.</p> * * <p>IEEE 854 also specifies the existence of "sub-normal" numbers. This * class does not contain any such entities. The most significant radix * 10000 digit is always non-zero. Instead, we support "gradual underflow" * by raising the underflow flag for numbers less with exponent less than * expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits. * Thus the smallest number we can represent would be: * 1E(-(MIN_EXP-digits-1)*4), eg, for digits=5, MIN_EXP=-32767, that would * be 1e-131092.</p> * * <p>IEEE 854 defines that the implied radix point lies just to the right * of the most significant digit and to the left of the remaining digits. * This implementation puts the implied radix point to the left of all * digits including the most significant one. The most significant digit * here is the one just to the right of the radix point. This is a fine * detail and is really only a matter of definition. Any side effects of * this can be rendered invisible by a subclass.</p> * @see DfpField * @version $Id: Dfp.java 1416643 2012-12-03 19:37:14Z tn $ * @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: 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 if ((bits & 0x8000000000000000L) != 0) { sign = -1; } return; } exponent++; // Normalize the subnormal number while ((mantissa & 0x0010000000000000L) == 0) { exponent--; mantissa <<= 1; } mantissa &= 0x000fffffffffffffL; } if (exponent == 1024) { // infinity or NAN if (x != x) { sign = (byte) 1; nans = QNAN; } else if (x < 0) { sign = (byte) -1; nans = INFINITE; } else { sign = (byte) 1; nans = INFINITE; } return; } Dfp xdfp = new Dfp(field, mantissa); xdfp = xdfp.divide(new Dfp(field, 4503599627370496l)).add(field.getOne()); // Divide by 2^52, then add one xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent)); if ((bits & 0x8000000000000000L) != 0) { xdfp = xdfp.negate(); } System.arraycopy(xdfp.mant, 0, mant, 0, mant.length); sign = xdfp.sign; exp = xdfp.exp; nans = xdfp.nans; } /** Copy constructor. * @param d instance to copy */ public Dfp(final Dfp d) { mant = d.mant.clone(); sign = d.sign; exp = d.exp; nans = d.nans; field = d.field; } /** Create an instance from a String representation. * @param field field to which this instance belongs * @param s string representation of the instance */ protected Dfp(final DfpField field, final String s) { // initialize as if 0 mant = new int[field.getRadixDigits()]; sign = 1; 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; while (p > q) { if (striped[p] != '0') { break; } p--; } // Make sure the decimal is on a mod 10000 boundary int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize; q -= i; decimalPos += i; // Make the mantissa length right by adding zeros at the end if necessary while ((p - q) < (mant.length * rsize)) { for (i = 0; i < rsize; i++) { striped[++p] = '0'; } } // Ok, now we know how many trailing zeros there are, // and where the least significant digit is for (i = mant.length - 1; i >= 0; i--) { mant[i] = (striped[q] - '0') * 1000 + (striped[q + 1] - '0') * 100 + (striped[q + 2] - '0') * 10 + (striped[q + 3] - '0'); q += 4; } exp = (decimalPos + sciexp) / rsize; if (q < striped.length) { // Is there possible another digit? round((striped[q] - '0') * 1000); } } /** Creates an instance with a non-finite value. * @param field field to which this instance belongs * @param sign sign of the Dfp to create * @param nans code of the value, must be one of {@link #INFINITE}, * {@link #SNAN}, {@link #QNAN} */ protected Dfp(final DfpField field, final byte sign, final byte nans) { this.field = field; this.mant = new int[field.getRadixDigits()]; this.sign = sign; this.exp = 0; this.nans = nans; } /** Create an instance with a value of 0. * Use this internally in preference to constructors to facilitate subclasses * @return a new instance with a value of 0 */ public Dfp newInstance() { return new Dfp(getField()); } /** Create an instance from a byte value. * @param x value to convert to an instance * @return a new instance with value x */ public Dfp newInstance(final byte x) { return new Dfp(getField(), x); } /** Create an instance from an int value. * @param x value to convert to an instance * @return a new instance with value x */ public Dfp newInstance(final int x) { return new Dfp(getField(), x); } /** Create an instance from a long value. * @param x value to convert to an instance * @return a new instance with value x */ public Dfp newInstance(final long x) { return new Dfp(getField(), x); } /** 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.math3.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.math3.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 less than or equal to 0. * @return true if instance is not NaN and less than or equal to 0, false otherwise */ public boolean negativeOrNull() { if (isNaN()) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero())); return false; } return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite()); } /** Check if instance is strictly less than 0. * @return true if instance is not NaN and less than or equal to 0, false otherwise */ public boolean strictlyNegative() { if (isNaN()) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero())); return false; } return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite()); } /** Check if instance is greater than or equal to 0. * @return true if instance is not NaN and greater than or equal to 0, false otherwise */ public boolean positiveOrNull() { if (isNaN()) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero())); return false; } return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite()); } /** Check if instance is strictly greater than 0. * @return true if instance is not NaN and greater than or equal to 0, false otherwise */ public boolean strictlyPositive() { if (isNaN()) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero())); return false; } return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite()); } /** Get the absolute value of instance. * @return absolute value of instance */ public Dfp abs() { Dfp result = newInstance(this); result.sign = 1; return result; } /** 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 zero. * @return true if instance is equal to zero */ public boolean isZero() { if (isNaN()) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero())); return false; } return (mant[mant.length - 1] == 0) && !isInfinite(); } /** Check if instance is equal to x. * @param other object to check instance against * @return true if instance is equal to x and neither are NaN, false otherwise */ @Override public boolean equals(final Object other) { if (other instanceof Dfp) { final Dfp x = (Dfp) other; if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) { return false; } return compare(this, x) == 0; } return false; } /** * Gets a hashCode for the instance. * @return a hash code value for this object */ @Override public int hashCode() { return 17 + (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 × 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 x. * @param x multiplicand * @return product of this and x */ public Dfp multiply(final int x) { if (x >= 0 && x < RADIX) { return multiplyFast(x); } else { return multiply(newInstance(x)); } } /** Multiply this by a single digit 0<=x<radix. * There are speed advantages in this special case. * @param x multiplicand * @return product of this and x */ private Dfp multiplyFast(final int x) { Dfp result = newInstance(this); /* handle special cases */ if (nans != FINITE) { if (isNaN()) { return this; } if (nans == INFINITE && x != 0) { result = newInstance(this); return result; } if (nans == INFINITE && x == 0) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); result = newInstance(getZero()); result.nans = QNAN; result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result); return result; } } /* range check x */ if (x < 0 || x >= RADIX) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); result = newInstance(getZero()); result.nans = QNAN; result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result); return result; } int rh = 0; for (int i = 0; i < mant.length; i++) { final int r = mant[i] * x + rh; rh = r / RADIX; result.mant[i] = r - rh * RADIX; } int lostdigit = 0; if (rh != 0) { lostdigit = result.mant[0]; result.shiftRight(); result.mant[mant.length - 1] = rh; } if (result.mant[mant.length - 1] == 0) { // if result is zero, set exp to zero result.exp = 0; } final int excp = result.round(lostdigit); if (excp != 0) { result = dotrap(excp, MULTIPLY_TRAP, result, result); } return result; } /** 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 <= 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; } /** {@inheritDoc} */ public Dfp reciprocal() { return field.getOne().divide(this); } /** 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; int cmp0 = compare(this, getZero()); if (cmp0 == 0) { return sign < 0 ? -0.0 : +0.0; } else if (cmp0 < 0) { y = negate(); negate = true; } /* Find the exponent, first estimate by integer log10, then adjust. Should be faster than doing a natural logarithm. */ int exponent = (int) (y.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; } }