home | career | drupal | java | mac | mysql | perl | scala | uml | unix  

Java example source code file (Dfp.java)

This example Java source code file (Dfp.java) is included in the alvinalexander.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

Learn more about this Java project at its project page.

Java - Java tags/keywords

add_trap, dfp, dfpfield, divide_trap, finite, infinite, less_than_trap, multiply_trap, next_after_trap, qnan, radix, snan, string, trunc_trap, util

The Dfp.java Java example source code

/*
 * 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.RealFieldElement;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.util.FastMath;

/**
 *  Decimal floating point library for Java
 *
 *  <p>Another floating point class.  This one is built using radix 10000
 *  which is 10<sup>4, so its almost decimal.

* * <p>The design goals here are: * <ol> * <li>Decimal math, or close to it * <li>Settable precision (but no mix between numbers using different settings) * <li>Portability. Code should be kept as portable as possible. * <li>Performance * <li>Accuracy - Results should always be +/- 1 ULP for basic * algebraic operation</li> * <li>Comply with IEEE 854-1987 as much as possible. * (See IEEE 854-1987 notes below)</li> * </ol>

* * <p>Trade offs: * <ol> * <li>Memory foot print. I'm using more memory than necessary to * represent numbers to get better performance.</li> * <li>Digits are bigger, so rounding is a greater loss. So, if you * really need 12 decimal digits, better use 4 base 10000 digits * there can be one partially filled.</li> * </ol>

* * <p>Numbers are represented in the following form: * <pre> * n = sign × mant × (radix)<sup>exp;

* </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>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 * @since 2.2 */ public class Dfp implements RealFieldElement<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 * @since 3.2 */ 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 + (isZero() ? 0 : (sign << 8)) + (nans << 16) + exp + Arrays.hashCode(mant); } /** Check if instance is not equal to x. * @param x number to check instance against * @return true if instance is not equal to x and neither are NaN, false otherwise */ public boolean unequal(final Dfp x) { if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) { return false; } return greaterThan(x) || lessThan(x); } /** Compare two instances. * @param a first instance in comparison * @param b second instance in comparison * @return -1 if a<b, 1 if a>b and 0 if a==b * Note this method does not properly handle NaNs or numbers with different precision. */ private static int compare(final Dfp a, final Dfp b) { // Ignore the sign of zero if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 && a.nans == FINITE && b.nans == FINITE) { return 0; } if (a.sign != b.sign) { if (a.sign == -1) { return -1; } else { return 1; } } // deal with the infinities if (a.nans == INFINITE && b.nans == FINITE) { return a.sign; } if (a.nans == FINITE && b.nans == INFINITE) { return -b.sign; } if (a.nans == INFINITE && b.nans == INFINITE) { return 0; } // Handle special case when a or b is zero, by ignoring the exponents if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) { if (a.exp < b.exp) { return -a.sign; } if (a.exp > b.exp) { return a.sign; } } // compare the mantissas for (int i = a.mant.length - 1; i >= 0; i--) { if (a.mant[i] > b.mant[i]) { return a.sign; } if (a.mant[i] < b.mant[i]) { return -a.sign; } } return 0; } /** Round to nearest integer using the round-half-even method. * That is round to nearest integer unless both are equidistant. * In which case round to the even one. * @return rounded value * @since 3.2 */ 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 * @since 3.2 */ 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 * @since 3.2 */ 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 * @since 3.2 */ 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 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 */ 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 * @since 3.2 */ public int intLog10() { if (mant[mant.length-1] > 1000) { return exp * 4 - 1; } if (mant[mant.length-1] > 100) { return exp * 4 - 2; } if (mant[mant.length-1] > 10) { return exp * 4 - 3; } return exp * 4 - 4; } /** Return the specified power of 10. * @param e desired power * @return 10<sup>e */ 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 -= 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 += 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 /= 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 * @since 3.2 */ 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 %= 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 += ERR_SCALE; break; case DfpField.FLAG_OVERFLOW: result.exp -= ERR_SCALE; def = newInstance(getZero()); def.sign = result.sign; def.nans = INFINITE; break; default: def = result; break; } return trap(type, what, oper, def, result); } /** Trap handler. Subclasses may override this to provide trap * functionality per IEEE 854-1987. * * @param type The exception type - e.g. FLAG_OVERFLOW * @param what The name of the routine we were in e.g. divide() * @param oper An operand to this function if any * @param def The default return value if trap not enabled * @param result The result that is specified to be delivered per * IEEE 854, if any * @return the value that should be return by the operation triggering the trap */ protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) { return def; } /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN. * @return type of the number */ public int classify() { return nans; } /** Creates an instance that is the same as x except that it has the sign of y. * abs(x) = dfp.copysign(x, dfp.one) * @param x number to get the value from * @param y number to get the sign from * @return a number with the value of x and the sign of y */ public static Dfp copysign(final Dfp x, final Dfp y) { Dfp result = x.newInstance(x); result.sign = y.sign; return result; } /** Returns the next number greater than this one in the direction of x. * If this==x then simply returns this. * @param x direction where to look at * @return closest number next to instance in the direction of x */ public Dfp nextAfter(final Dfp x) { // make sure we don't mix number with different precision if (field.getRadixDigits() != x.field.getRadixDigits()) { field.setIEEEFlagsBits(DfpField.FLAG_INVALID); final Dfp result = newInstance(getZero()); result.nans = QNAN; return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result); } // if this is greater than x boolean up = false; if (this.lessThan(x)) { up = true; } if (compare(this, x) == 0) { return newInstance(x); } if (lessThan(getZero())) { up = !up; } final Dfp inc; Dfp result; if (up) { inc = newInstance(getOne()); inc.exp = this.exp-mant.length+1; inc.sign = this.sign; if (this.equals(getZero())) { inc.exp = MIN_EXP-mant.length; } result = add(inc); } else { inc = newInstance(getOne()); inc.exp = this.exp; inc.sign = this.sign; if (this.equals(inc)) { inc.exp = this.exp-mant.length; } else { inc.exp = this.exp-mant.length+1; } if (this.equals(getZero())) { inc.exp = MIN_EXP-mant.length; } result = this.subtract(inc); } if (result.classify() == INFINITE && this.classify() != INFINITE) { field.setIEEEFlagsBits(DfpField.FLAG_INEXACT); result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result); } if (result.equals(getZero()) && this.equals(getZero()) == 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.intLog10() * 3.32); if (exponent < 0) { exponent--; } Dfp tempDfp = DfpMath.pow(getTwo(), exponent); while (tempDfp.lessThan(y) || tempDfp.equals(y)) { tempDfp = tempDfp.multiply(2); exponent++; } exponent--; /* We have the exponent, now work on the mantissa */ y = y.divide(DfpMath.pow(getTwo(), exponent)); if (exponent > -1023) { y = y.subtract(getOne()); } if (exponent < -1074) { return 0; } if (exponent > 1023) { return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; } y = y.multiply(newInstance(4503599627370496l)).rint(); String str = y.toString(); str = str.substring(0, str.length()-1); long mantissa = Long.parseLong(str); if (mantissa == 4503599627370496L) { // Handle special case where we round up to next power of two mantissa = 0; exponent++; } /* Its going to be subnormal, so make adjustments */ if (exponent <= -1023) { exponent--; } while (exponent < -1023) { exponent++; mantissa >>>= 1; } long bits = mantissa | ((exponent + 1023L) << 52); double x = Double.longBitsToDouble(bits); if (negate) { x = -x; } return x; } /** Convert the instance into a split double. * @return an array of two doubles which sum represent the instance * @see #toDouble() */ public double[] toSplitDouble() { double split[] = new double[2]; long mask = 0xffffffffc0000000L; split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask); split[1] = subtract(newInstance(split[0])).toDouble(); return split; } /** {@inheritDoc} * @since 3.2 */ public double getReal() { return toDouble(); } /** {@inheritDoc} * @since 3.2 */ public Dfp add(final double a) { return add(newInstance(a)); } /** {@inheritDoc} * @since 3.2 */ public Dfp subtract(final double a) { return subtract(newInstance(a)); } /** {@inheritDoc} * @since 3.2 */ public Dfp multiply(final double a) { return multiply(newInstance(a)); } /** {@inheritDoc} * @since 3.2 */ public Dfp divide(final double a) { return divide(newInstance(a)); } /** {@inheritDoc} * @since 3.2 */ public Dfp remainder(final double a) { return remainder(newInstance(a)); } /** {@inheritDoc} * @since 3.2 */ public long round() { return FastMath.round(toDouble()); } /** {@inheritDoc} * @since 3.2 */ public Dfp signum() { if (isNaN() || isZero()) { return this; } else { return newInstance(sign > 0 ? +1 : -1); } } /** {@inheritDoc} * @since 3.2 */ public Dfp copySign(final Dfp s) { if ((sign >= 0 && s.sign >= 0) || (sign < 0 && s.sign < 0)) { // Sign is currently OK return this; } return negate(); // flip sign } /** {@inheritDoc} * @since 3.2 */ public Dfp copySign(final double s) { long sb = Double.doubleToLongBits(s); if ((sign >= 0 && sb >= 0) || (sign < 0 && sb < 0)) { // Sign is currently OK return this; } return negate(); // flip sign } /** {@inheritDoc} * @since 3.2 */ public Dfp scalb(final int n) { return multiply(DfpMath.pow(getTwo(), n)); } /** {@inheritDoc} * @since 3.2 */ public Dfp hypot(final Dfp y) { return multiply(this).add(y.multiply(y)).sqrt(); } /** {@inheritDoc} * @since 3.2 */ public Dfp cbrt() { return rootN(3); } /** {@inheritDoc} * @since 3.2 */ public Dfp rootN(final int n) { return (sign >= 0) ? DfpMath.pow(this, getOne().divide(n)) : DfpMath.pow(negate(), getOne().divide(n)).negate(); } /** {@inheritDoc} * @since 3.2 */ public Dfp pow(final double p) { return DfpMath.pow(this, newInstance(p)); } /** {@inheritDoc} * @since 3.2 */ public Dfp pow(final int n) { return DfpMath.pow(this, n); } /** {@inheritDoc} * @since 3.2 */ public Dfp pow(final Dfp e) { return DfpMath.pow(this, e); } /** {@inheritDoc} * @since 3.2 */ public Dfp exp() { return DfpMath.exp(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp expm1() { return DfpMath.exp(this).subtract(getOne()); } /** {@inheritDoc} * @since 3.2 */ public Dfp log() { return DfpMath.log(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp log1p() { return DfpMath.log(this.add(getOne())); } // TODO: deactivate this implementation (and return type) in 4.0 /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this). * @return integer base 10 logarithm * @deprecated as of 3.2, replaced by {@link #intLog10()}, in 4.0 the return type * will be changed to Dfp */ @Deprecated public int log10() { return intLog10(); } // TODO: activate this implementation (and return type) in 4.0 // /** {@inheritDoc} // * @since 3.2 // */ // public Dfp log10() { // return DfpMath.log(this).divide(DfpMath.log(newInstance(10))); // } /** {@inheritDoc} * @since 3.2 */ public Dfp cos() { return DfpMath.cos(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp sin() { return DfpMath.sin(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp tan() { return DfpMath.tan(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp acos() { return DfpMath.acos(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp asin() { return DfpMath.asin(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp atan() { return DfpMath.atan(this); } /** {@inheritDoc} * @since 3.2 */ public Dfp atan2(final Dfp x) throws DimensionMismatchException { // compute r = sqrt(x^2+y^2) final Dfp r = x.multiply(x).add(multiply(this)).sqrt(); if (x.sign >= 0) { // compute atan2(y, x) = 2 atan(y / (r + x)) return getTwo().multiply(divide(r.add(x)).atan()); } else { // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x)) final Dfp tmp = getTwo().multiply(divide(r.subtract(x)).atan()); final Dfp pmPi = newInstance((tmp.sign <= 0) ? -FastMath.PI : FastMath.PI); return pmPi.subtract(tmp); } } /** {@inheritDoc} * @since 3.2 */ public Dfp cosh() { return DfpMath.exp(this).add(DfpMath.exp(negate())).divide(2); } /** {@inheritDoc} * @since 3.2 */ public Dfp sinh() { return DfpMath.exp(this).subtract(DfpMath.exp(negate())).divide(2); } /** {@inheritDoc} * @since 3.2 */ public Dfp tanh() { final Dfp ePlus = DfpMath.exp(this); final Dfp eMinus = DfpMath.exp(negate()); return ePlus.subtract(eMinus).divide(ePlus.add(eMinus)); } /** {@inheritDoc} * @since 3.2 */ public Dfp acosh() { return multiply(this).subtract(getOne()).sqrt().add(this).log(); } /** {@inheritDoc} * @since 3.2 */ public Dfp asinh() { return multiply(this).add(getOne()).sqrt().add(this).log(); } /** {@inheritDoc} * @since 3.2 */ public Dfp atanh() { return getOne().add(this).divide(getOne().subtract(this)).log().divide(2); } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final Dfp[] a, final Dfp[] b) throws DimensionMismatchException { if (a.length != b.length) { throw new DimensionMismatchException(a.length, b.length); } Dfp r = getZero(); for (int i = 0; i < a.length; ++i) { r = r.add(a[i].multiply(b[i])); } return r; } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final double[] a, final Dfp[] b) throws DimensionMismatchException { if (a.length != b.length) { throw new DimensionMismatchException(a.length, b.length); } Dfp r = getZero(); for (int i = 0; i < a.length; ++i) { r = r.add(b[i].multiply(a[i])); } return r; } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2) { return a1.multiply(b1).add(a2.multiply(b2)); } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2) { return b1.multiply(a1).add(b2.multiply(a2)); } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2, final Dfp a3, final Dfp b3) { return a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)); } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2, final double a3, final Dfp b3) { return b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)); } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2, final Dfp a3, final Dfp b3, final Dfp a4, final Dfp b4) { return a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4)); } /** {@inheritDoc} * @since 3.2 */ public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2, final double a3, final Dfp b3, final double a4, final Dfp b4) { return b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4)); } }

Other Java examples (source code examples)

Here is a short list of links related to this Java Dfp.java source code file:



my book on functional programming

 

new blog posts

 

Copyright 1998-2019 Alvin Alexander, alvinalexander.com
All Rights Reserved.

A percentage of advertising revenue from
pages under the /java/jwarehouse URI on this website is
paid back to open source projects.