diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/dfp')
7 files changed, 5547 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java b/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java new file mode 100644 index 0000000..23d4861 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/BracketingNthOrderBrentSolverDFP.java @@ -0,0 +1,161 @@ +/* + * 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 org.apache.commons.math3.analysis.RealFieldUnivariateFunction; +import org.apache.commons.math3.analysis.solvers.AllowedSolution; +import org.apache.commons.math3.analysis.solvers.FieldBracketingNthOrderBrentSolver; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.util.MathUtils; + +/** + * This class implements a modification of the <a + * href="http://mathworld.wolfram.com/BrentsMethod.html">Brent algorithm</a>. + * + * <p>The changes with respect to the original Brent algorithm are: + * + * <ul> + * <li>the returned value is chosen in the current interval according to user specified {@link + * AllowedSolution}, + * <li>the maximal order for the invert polynomial root search is user-specified instead of being + * invert quadratic only + * </ul> + * + * The given interval must bracket the root. + * + * @deprecated as of 3.6 replaced with {@link FieldBracketingNthOrderBrentSolver} + */ +@Deprecated +public class BracketingNthOrderBrentSolverDFP extends FieldBracketingNthOrderBrentSolver<Dfp> { + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Function value accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolverDFP( + final Dfp relativeAccuracy, + final Dfp absoluteAccuracy, + final Dfp functionValueAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy, maximalOrder); + } + + /** + * Get the absolute accuracy. + * + * @return absolute accuracy + */ + @Override + public Dfp getAbsoluteAccuracy() { + return super.getAbsoluteAccuracy(); + } + + /** + * Get the relative accuracy. + * + * @return relative accuracy + */ + @Override + public Dfp getRelativeAccuracy() { + return super.getRelativeAccuracy(); + } + + /** + * Get the function accuracy. + * + * @return function accuracy + */ + @Override + public Dfp getFunctionValueAccuracy() { + return super.getFunctionValueAccuracy(); + } + + /** + * Solve for a zero in the given interval. A solver may require that the interval brackets a + * single zero root. Solvers that do require bracketing should be able to handle the case where + * one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param allowedSolution The kind of solutions that the root-finding algorithm may accept as + * solutions. + * @return a value where the function is zero. + * @exception NullArgumentException if f is null. + * @exception NoBracketingException if root cannot be bracketed + */ + public Dfp solve( + final int maxEval, + final UnivariateDfpFunction f, + final Dfp min, + final Dfp max, + final AllowedSolution allowedSolution) + throws NullArgumentException, NoBracketingException { + return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution); + } + + /** + * Solve for a zero in the given interval, start at {@code startValue}. A solver may require + * that the interval brackets a single zero root. Solvers that do require bracketing should be + * able to handle the case where one of the endpoints is itself a root. + * + * @param maxEval Maximum number of evaluations. + * @param f Function to solve. + * @param min Lower bound for the interval. + * @param max Upper bound for the interval. + * @param startValue Start value to use. + * @param allowedSolution The kind of solutions that the root-finding algorithm may accept as + * solutions. + * @return a value where the function is zero. + * @exception NullArgumentException if f is null. + * @exception NoBracketingException if root cannot be bracketed + */ + public Dfp solve( + final int maxEval, + final UnivariateDfpFunction f, + final Dfp min, + final Dfp max, + final Dfp startValue, + final AllowedSolution allowedSolution) + throws NullArgumentException, NoBracketingException { + + // checks + MathUtils.checkNotNull(f); + + // wrap the function + RealFieldUnivariateFunction<Dfp> fieldF = + new RealFieldUnivariateFunction<Dfp>() { + + /** {@inheritDoc} */ + public Dfp value(final Dfp x) { + return f.value(x); + } + }; + + // delegate to general field solver + return solve(maxEval, fieldF, min, max, startValue, allowedSolution); + } +} diff --git a/src/main/java/org/apache/commons/math3/dfp/Dfp.java b/src/main/java/org/apache/commons/math3/dfp/Dfp.java new file mode 100644 index 0000000..4f59a19 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/Dfp.java @@ -0,0 +1,3082 @@ +/* + * 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 org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.util.FastMath; + +import java.util.Arrays; + +/** + * 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>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>Comply with IEEE 854-1987 as much as possible. (See IEEE 854-1987 notes below) + * </ol> + * + * <p>Trade offs: + * + * <ol> + * <li>Memory foot print. I'm using more memory than necessary to represent numbers to get better + * performance. + * <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. + * </ol> + * + * <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>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>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>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>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>IEEE 854 defines that the implied radix point lies just to the right of the most significant + * digit and to the left of the remaining digits. This implementation puts the implied radix point + * to the left of all digits including the most significant one. The most significant digit here is + * the one just to the right of the radix point. This is a fine detail and is really only a matter + * of definition. Any side effects of this can be rendered invisible by a subclass. + * + * @see DfpField + * @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. + * + * @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</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 + * @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</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 -= 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)); + } +} diff --git a/src/main/java/org/apache/commons/math3/dfp/DfpDec.java b/src/main/java/org/apache/commons/math3/dfp/DfpDec.java new file mode 100644 index 0000000..5571c2d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/DfpDec.java @@ -0,0 +1,388 @@ +/* + * 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; + +/** + * Subclass of {@link Dfp} which hides the radix-10000 artifacts of the superclass. This should give + * outward appearances of being a decimal number with DIGITS*4-3 decimal digits. This class can be + * subclassed to appear to be an arbitrary number of decimal digits less than DIGITS*4-3. + * + * @since 2.2 + */ +public class DfpDec extends Dfp { + + /** + * Makes an instance with a value of zero. + * + * @param factory factory linked to this instance + */ + protected DfpDec(final DfpField factory) { + super(factory); + } + + /** + * Create an instance from a byte value. + * + * @param factory factory linked to this instance + * @param x value to convert to an instance + */ + protected DfpDec(final DfpField factory, byte x) { + super(factory, x); + } + + /** + * Create an instance from an int value. + * + * @param factory factory linked to this instance + * @param x value to convert to an instance + */ + protected DfpDec(final DfpField factory, int x) { + super(factory, x); + } + + /** + * Create an instance from a long value. + * + * @param factory factory linked to this instance + * @param x value to convert to an instance + */ + protected DfpDec(final DfpField factory, long x) { + super(factory, x); + } + + /** + * Create an instance from a double value. + * + * @param factory factory linked to this instance + * @param x value to convert to an instance + */ + protected DfpDec(final DfpField factory, double x) { + super(factory, x); + round(0); + } + + /** + * Copy constructor. + * + * @param d instance to copy + */ + public DfpDec(final Dfp d) { + super(d); + round(0); + } + + /** + * Create an instance from a String representation. + * + * @param factory factory linked to this instance + * @param s string representation of the instance + */ + protected DfpDec(final DfpField factory, final String s) { + super(factory, s); + round(0); + } + + /** + * Creates an instance with a non-finite value. + * + * @param factory factory linked to this instance + * @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 DfpDec(final DfpField factory, final byte sign, final byte nans) { + super(factory, sign, nans); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance() { + return new DfpDec(getField()); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final byte x) { + return new DfpDec(getField(), x); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final int x) { + return new DfpDec(getField(), x); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final long x) { + return new DfpDec(getField(), x); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final double x) { + return new DfpDec(getField(), x); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final Dfp d) { + + // make sure we don't mix number with different precision + if (getField().getRadixDigits() != d.getField().getRadixDigits()) { + getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + final Dfp result = newInstance(getZero()); + result.nans = QNAN; + return dotrap(DfpField.FLAG_INVALID, "newInstance", d, result); + } + + return new DfpDec(d); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final String s) { + return new DfpDec(getField(), s); + } + + /** {@inheritDoc} */ + @Override + public Dfp newInstance(final byte sign, final byte nans) { + return new DfpDec(getField(), sign, nans); + } + + /** + * Get the number of decimal digits this class is going to represent. Default implementation + * returns {@link #getRadixDigits()}*4-3. Subclasses can override this to return something less. + * + * @return number of decimal digits this class is going to represent + */ + protected int getDecimalDigits() { + return getRadixDigits() * 4 - 3; + } + + /** {@inheritDoc} */ + @Override + protected int round(int in) { + + int msb = mant[mant.length - 1]; + if (msb == 0) { + // special case -- this == zero + return 0; + } + + int cmaxdigits = mant.length * 4; + int lsbthreshold = 1000; + while (lsbthreshold > msb) { + lsbthreshold /= 10; + cmaxdigits--; + } + + final int digits = getDecimalDigits(); + final int lsbshift = cmaxdigits - digits; + final int lsd = lsbshift / 4; + + lsbthreshold = 1; + for (int i = 0; i < lsbshift % 4; i++) { + lsbthreshold *= 10; + } + + final int lsb = mant[lsd]; + + if (lsbthreshold <= 1 && digits == 4 * mant.length - 3) { + return super.round(in); + } + + int discarded = in; // not looking at this after this point + final int n; + if (lsbthreshold == 1) { + // look to the next digit for rounding + n = (mant[lsd - 1] / 1000) % 10; + mant[lsd - 1] %= 1000; + discarded |= mant[lsd - 1]; + } else { + n = (lsb * 10 / lsbthreshold) % 10; + discarded |= lsb % (lsbthreshold / 10); + } + + for (int i = 0; i < lsd; i++) { + discarded |= mant[i]; // need to know if there are any discarded bits + mant[i] = 0; + } + + mant[lsd] = lsb / lsbthreshold * lsbthreshold; + + final boolean inc; + switch (getField().getRoundingMode()) { + case ROUND_DOWN: + inc = false; + break; + + case ROUND_UP: + inc = (n != 0) || (discarded != 0); // round up if n!=0 + break; + + case ROUND_HALF_UP: + inc = n >= 5; // round half up + break; + + case ROUND_HALF_DOWN: + inc = n > 5; // round half down + break; + + case ROUND_HALF_EVEN: + inc = + (n > 5) + || (n == 5 && discarded != 0) + || (n == 5 + && discarded == 0 + && ((lsb / lsbthreshold) & 1) == 1); // round half-even + break; + + case ROUND_HALF_ODD: + inc = + (n > 5) + || (n == 5 && discarded != 0) + || (n == 5 + && discarded == 0 + && ((lsb / lsbthreshold) & 1) == 0); // round half-odd + break; + + case ROUND_CEIL: + inc = (sign == 1) && (n != 0 || discarded != 0); // round ceil + break; + + case ROUND_FLOOR: + default: + inc = (sign == -1) && (n != 0 || discarded != 0); // round floor + break; + } + + if (inc) { + // increment if necessary + int rh = lsbthreshold; + for (int i = lsd; i < mant.length; i++) { + final int r = mant[i] + rh; + rh = r / RADIX; + mant[i] = r % 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 + getField().setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW); + return DfpField.FLAG_UNDERFLOW; + } + + if (exp > MAX_EXP) { + // Overflow + getField().setIEEEFlagsBits(DfpField.FLAG_OVERFLOW); + return DfpField.FLAG_OVERFLOW; + } + + if (n != 0 || discarded != 0) { + // Inexact + getField().setIEEEFlagsBits(DfpField.FLAG_INEXACT); + return DfpField.FLAG_INEXACT; + } + return 0; + } + + /** {@inheritDoc} */ + @Override + public Dfp nextAfter(Dfp x) { + + final String trapName = "nextAfter"; + + // make sure we don't mix number with different precision + if (getField().getRadixDigits() != x.getField().getRadixDigits()) { + getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + final Dfp result = newInstance(getZero()); + result.nans = QNAN; + return dotrap(DfpField.FLAG_INVALID, trapName, x, result); + } + + boolean up = false; + Dfp result; + Dfp inc; + + // if this is greater than x + if (this.lessThan(x)) { + up = true; + } + + if (equals(x)) { + return newInstance(x); + } + + if (lessThan(getZero())) { + up = !up; + } + + if (up) { + inc = power10(intLog10() - getDecimalDigits() + 1); + inc = copysign(inc, this); + + if (this.equals(getZero())) { + inc = power10K(MIN_EXP - mant.length - 1); + } + + if (inc.equals(getZero())) { + result = copysign(newInstance(getZero()), this); + } else { + result = add(inc); + } + } else { + inc = power10(intLog10()); + inc = copysign(inc, this); + + if (this.equals(inc)) { + inc = inc.divide(power10(getDecimalDigits())); + } else { + inc = inc.divide(power10(getDecimalDigits() - 1)); + } + + if (this.equals(getZero())) { + inc = power10K(MIN_EXP - mant.length - 1); + } + + if (inc.equals(getZero())) { + result = copysign(newInstance(getZero()), this); + } else { + result = subtract(inc); + } + } + + if (result.classify() == INFINITE && this.classify() != INFINITE) { + getField().setIEEEFlagsBits(DfpField.FLAG_INEXACT); + result = dotrap(DfpField.FLAG_INEXACT, trapName, x, result); + } + + if (result.equals(getZero()) && this.equals(getZero()) == false) { + getField().setIEEEFlagsBits(DfpField.FLAG_INEXACT); + result = dotrap(DfpField.FLAG_INEXACT, trapName, x, result); + } + + return result; + } +} diff --git a/src/main/java/org/apache/commons/math3/dfp/DfpField.java b/src/main/java/org/apache/commons/math3/dfp/DfpField.java new file mode 100644 index 0000000..bd5a7eb --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/DfpField.java @@ -0,0 +1,826 @@ +/* + * 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 org.apache.commons.math3.Field; +import org.apache.commons.math3.FieldElement; + +/** + * Field for Decimal floating point instances. + * + * @since 2.2 + */ +public class DfpField implements Field<Dfp> { + + /** Enumerate for rounding modes. */ + public enum RoundingMode { + + /** Rounds toward zero (truncation). */ + ROUND_DOWN, + + /** Rounds away from zero if discarded digit is non-zero. */ + ROUND_UP, + + /** + * Rounds towards nearest unless both are equidistant in which case it rounds away from + * zero. + */ + ROUND_HALF_UP, + + /** + * Rounds towards nearest unless both are equidistant in which case it rounds toward zero. + */ + ROUND_HALF_DOWN, + + /** + * Rounds towards nearest unless both are equidistant in which case it rounds toward the + * even neighbor. This is the default as specified by IEEE 854-1987 + */ + ROUND_HALF_EVEN, + + /** + * Rounds towards nearest unless both are equidistant in which case it rounds toward the odd + * neighbor. + */ + ROUND_HALF_ODD, + + /** Rounds towards positive infinity. */ + ROUND_CEIL, + + /** Rounds towards negative infinity. */ + ROUND_FLOOR; + } + + /** IEEE 854-1987 flag for invalid operation. */ + public static final int FLAG_INVALID = 1; + + /** IEEE 854-1987 flag for division by zero. */ + public static final int FLAG_DIV_ZERO = 2; + + /** IEEE 854-1987 flag for overflow. */ + public static final int FLAG_OVERFLOW = 4; + + /** IEEE 854-1987 flag for underflow. */ + public static final int FLAG_UNDERFLOW = 8; + + /** IEEE 854-1987 flag for inexact result. */ + public static final int FLAG_INEXACT = 16; + + /** High precision string representation of √2. */ + private static String sqr2String; + + // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class") + + /** High precision string representation of √2 / 2. */ + private static String sqr2ReciprocalString; + + /** High precision string representation of √3. */ + private static String sqr3String; + + /** High precision string representation of √3 / 3. */ + private static String sqr3ReciprocalString; + + /** High precision string representation of π. */ + private static String piString; + + /** High precision string representation of e. */ + private static String eString; + + /** High precision string representation of ln(2). */ + private static String ln2String; + + /** High precision string representation of ln(5). */ + private static String ln5String; + + /** High precision string representation of ln(10). */ + private static String ln10String; + + /** + * The number of radix digits. Note these depend on the radix which is 10000 digits, so each one + * is equivalent to 4 decimal digits. + */ + private final int radixDigits; + + /** A {@link Dfp} with value 0. */ + private final Dfp zero; + + /** A {@link Dfp} with value 1. */ + private final Dfp one; + + /** A {@link Dfp} with value 2. */ + private final Dfp two; + + /** A {@link Dfp} with value √2. */ + private final Dfp sqr2; + + /** A two elements {@link Dfp} array with value √2 split in two pieces. */ + private final Dfp[] sqr2Split; + + /** A {@link Dfp} with value √2 / 2. */ + private final Dfp sqr2Reciprocal; + + /** A {@link Dfp} with value √3. */ + private final Dfp sqr3; + + /** A {@link Dfp} with value √3 / 3. */ + private final Dfp sqr3Reciprocal; + + /** A {@link Dfp} with value π. */ + private final Dfp pi; + + /** A two elements {@link Dfp} array with value π split in two pieces. */ + private final Dfp[] piSplit; + + /** A {@link Dfp} with value e. */ + private final Dfp e; + + /** A two elements {@link Dfp} array with value e split in two pieces. */ + private final Dfp[] eSplit; + + /** A {@link Dfp} with value ln(2). */ + private final Dfp ln2; + + /** A two elements {@link Dfp} array with value ln(2) split in two pieces. */ + private final Dfp[] ln2Split; + + /** A {@link Dfp} with value ln(5). */ + private final Dfp ln5; + + /** A two elements {@link Dfp} array with value ln(5) split in two pieces. */ + private final Dfp[] ln5Split; + + /** A {@link Dfp} with value ln(10). */ + private final Dfp ln10; + + /** Current rounding mode. */ + private RoundingMode rMode; + + /** IEEE 854-1987 signals. */ + private int ieeeFlags; + + /** + * Create a factory for the specified number of radix digits. + * + * <p>Note that since the {@link Dfp} class uses 10000 as its radix, each radix digit is + * equivalent to 4 decimal digits. This implies that asking for 13, 14, 15 or 16 decimal digits + * will really lead to a 4 radix 10000 digits in all cases. + * + * @param decimalDigits minimal number of decimal digits. + */ + public DfpField(final int decimalDigits) { + this(decimalDigits, true); + } + + /** + * Create a factory for the specified number of radix digits. + * + * <p>Note that since the {@link Dfp} class uses 10000 as its radix, each radix digit is + * equivalent to 4 decimal digits. This implies that asking for 13, 14, 15 or 16 decimal digits + * will really lead to a 4 radix 10000 digits in all cases. + * + * @param decimalDigits minimal number of decimal digits + * @param computeConstants if true, the transcendental constants for the given precision must be + * computed (setting this flag to false is RESERVED for the internal recursive call) + */ + private DfpField(final int decimalDigits, final boolean computeConstants) { + + this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4; + this.rMode = RoundingMode.ROUND_HALF_EVEN; + this.ieeeFlags = 0; + this.zero = new Dfp(this, 0); + this.one = new Dfp(this, 1); + this.two = new Dfp(this, 2); + + if (computeConstants) { + // set up transcendental constants + synchronized (DfpField.class) { + + // as a heuristic to circumvent Table-Maker's Dilemma, we set the string + // representation of the constants to be at least 3 times larger than the + // number of decimal digits, also as an attempt to really compute these + // constants only once, we set a minimum number of digits + computeStringConstants((decimalDigits < 67) ? 200 : (3 * decimalDigits)); + + // set up the constants at current field accuracy + sqr2 = new Dfp(this, sqr2String); + sqr2Split = split(sqr2String); + sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString); + sqr3 = new Dfp(this, sqr3String); + sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString); + pi = new Dfp(this, piString); + piSplit = split(piString); + e = new Dfp(this, eString); + eSplit = split(eString); + ln2 = new Dfp(this, ln2String); + ln2Split = split(ln2String); + ln5 = new Dfp(this, ln5String); + ln5Split = split(ln5String); + ln10 = new Dfp(this, ln10String); + } + } else { + // dummy settings for unused constants + sqr2 = null; + sqr2Split = null; + sqr2Reciprocal = null; + sqr3 = null; + sqr3Reciprocal = null; + pi = null; + piSplit = null; + e = null; + eSplit = null; + ln2 = null; + ln2Split = null; + ln5 = null; + ln5Split = null; + ln10 = null; + } + } + + /** + * Get the number of radix digits of the {@link Dfp} instances built by this factory. + * + * @return number of radix digits + */ + public int getRadixDigits() { + return radixDigits; + } + + /** + * Set the rounding mode. If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}. + * + * @param mode desired rounding mode Note that the rounding mode is common to all {@link Dfp} + * instances belonging to the current {@link DfpField} in the system and will affect all + * future calculations. + */ + public void setRoundingMode(final RoundingMode mode) { + rMode = mode; + } + + /** + * Get the current rounding mode. + * + * @return current rounding mode + */ + public RoundingMode getRoundingMode() { + return rMode; + } + + /** + * Get the IEEE 854 status flags. + * + * @return IEEE 854 status flags + * @see #clearIEEEFlags() + * @see #setIEEEFlags(int) + * @see #setIEEEFlagsBits(int) + * @see #FLAG_INVALID + * @see #FLAG_DIV_ZERO + * @see #FLAG_OVERFLOW + * @see #FLAG_UNDERFLOW + * @see #FLAG_INEXACT + */ + public int getIEEEFlags() { + return ieeeFlags; + } + + /** + * Clears the IEEE 854 status flags. + * + * @see #getIEEEFlags() + * @see #setIEEEFlags(int) + * @see #setIEEEFlagsBits(int) + * @see #FLAG_INVALID + * @see #FLAG_DIV_ZERO + * @see #FLAG_OVERFLOW + * @see #FLAG_UNDERFLOW + * @see #FLAG_INEXACT + */ + public void clearIEEEFlags() { + ieeeFlags = 0; + } + + /** + * Sets the IEEE 854 status flags. + * + * @param flags desired value for the flags + * @see #getIEEEFlags() + * @see #clearIEEEFlags() + * @see #setIEEEFlagsBits(int) + * @see #FLAG_INVALID + * @see #FLAG_DIV_ZERO + * @see #FLAG_OVERFLOW + * @see #FLAG_UNDERFLOW + * @see #FLAG_INEXACT + */ + public void setIEEEFlags(final int flags) { + ieeeFlags = + flags + & (FLAG_INVALID + | FLAG_DIV_ZERO + | FLAG_OVERFLOW + | FLAG_UNDERFLOW + | FLAG_INEXACT); + } + + /** + * Sets some bits in the IEEE 854 status flags, without changing the already set bits. + * + * <p>Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)} + * + * @param bits bits to set + * @see #getIEEEFlags() + * @see #clearIEEEFlags() + * @see #setIEEEFlags(int) + * @see #FLAG_INVALID + * @see #FLAG_DIV_ZERO + * @see #FLAG_OVERFLOW + * @see #FLAG_UNDERFLOW + * @see #FLAG_INEXACT + */ + public void setIEEEFlagsBits(final int bits) { + ieeeFlags |= + bits + & (FLAG_INVALID + | FLAG_DIV_ZERO + | FLAG_OVERFLOW + | FLAG_UNDERFLOW + | FLAG_INEXACT); + } + + /** + * Makes a {@link Dfp} with a value of 0. + * + * @return a new {@link Dfp} with a value of 0 + */ + public Dfp newDfp() { + return new Dfp(this); + } + + /** + * Create an instance from a byte value. + * + * @param x value to convert to an instance + * @return a new {@link Dfp} with the same value as x + */ + public Dfp newDfp(final byte x) { + return new Dfp(this, x); + } + + /** + * Create an instance from an int value. + * + * @param x value to convert to an instance + * @return a new {@link Dfp} with the same value as x + */ + public Dfp newDfp(final int x) { + return new Dfp(this, x); + } + + /** + * Create an instance from a long value. + * + * @param x value to convert to an instance + * @return a new {@link Dfp} with the same value as x + */ + public Dfp newDfp(final long x) { + return new Dfp(this, x); + } + + /** + * Create an instance from a double value. + * + * @param x value to convert to an instance + * @return a new {@link Dfp} with the same value as x + */ + public Dfp newDfp(final double x) { + return new Dfp(this, x); + } + + /** + * Copy constructor. + * + * @param d instance to copy + * @return a new {@link Dfp} with the same value as d + */ + public Dfp newDfp(Dfp d) { + return new Dfp(d); + } + + /** + * Create a {@link Dfp} given a String representation. + * + * @param s string representation of the instance + * @return a new {@link Dfp} parsed from specified string + */ + public Dfp newDfp(final String s) { + return new Dfp(this, s); + } + + /** + * Creates a {@link Dfp} with a non-finite value. + * + * @param sign sign of the Dfp to create + * @param nans code of the value, must be one of {@link Dfp#INFINITE}, {@link Dfp#SNAN}, {@link + * Dfp#QNAN} + * @return a new {@link Dfp} with a non-finite value + */ + public Dfp newDfp(final byte sign, final byte nans) { + return new Dfp(this, sign, nans); + } + + /** + * Get the constant 0. + * + * @return a {@link Dfp} with value 0 + */ + public Dfp getZero() { + return zero; + } + + /** + * Get the constant 1. + * + * @return a {@link Dfp} with value 1 + */ + public Dfp getOne() { + return one; + } + + /** {@inheritDoc} */ + public Class<? extends FieldElement<Dfp>> getRuntimeClass() { + return Dfp.class; + } + + /** + * Get the constant 2. + * + * @return a {@link Dfp} with value 2 + */ + public Dfp getTwo() { + return two; + } + + /** + * Get the constant √2. + * + * @return a {@link Dfp} with value √2 + */ + public Dfp getSqr2() { + return sqr2; + } + + /** + * Get the constant √2 split in two pieces. + * + * @return a {@link Dfp} with value √2 split in two pieces + */ + public Dfp[] getSqr2Split() { + return sqr2Split.clone(); + } + + /** + * Get the constant √2 / 2. + * + * @return a {@link Dfp} with value √2 / 2 + */ + public Dfp getSqr2Reciprocal() { + return sqr2Reciprocal; + } + + /** + * Get the constant √3. + * + * @return a {@link Dfp} with value √3 + */ + public Dfp getSqr3() { + return sqr3; + } + + /** + * Get the constant √3 / 3. + * + * @return a {@link Dfp} with value √3 / 3 + */ + public Dfp getSqr3Reciprocal() { + return sqr3Reciprocal; + } + + /** + * Get the constant π. + * + * @return a {@link Dfp} with value π + */ + public Dfp getPi() { + return pi; + } + + /** + * Get the constant π split in two pieces. + * + * @return a {@link Dfp} with value π split in two pieces + */ + public Dfp[] getPiSplit() { + return piSplit.clone(); + } + + /** + * Get the constant e. + * + * @return a {@link Dfp} with value e + */ + public Dfp getE() { + return e; + } + + /** + * Get the constant e split in two pieces. + * + * @return a {@link Dfp} with value e split in two pieces + */ + public Dfp[] getESplit() { + return eSplit.clone(); + } + + /** + * Get the constant ln(2). + * + * @return a {@link Dfp} with value ln(2) + */ + public Dfp getLn2() { + return ln2; + } + + /** + * Get the constant ln(2) split in two pieces. + * + * @return a {@link Dfp} with value ln(2) split in two pieces + */ + public Dfp[] getLn2Split() { + return ln2Split.clone(); + } + + /** + * Get the constant ln(5). + * + * @return a {@link Dfp} with value ln(5) + */ + public Dfp getLn5() { + return ln5; + } + + /** + * Get the constant ln(5) split in two pieces. + * + * @return a {@link Dfp} with value ln(5) split in two pieces + */ + public Dfp[] getLn5Split() { + return ln5Split.clone(); + } + + /** + * Get the constant ln(10). + * + * @return a {@link Dfp} with value ln(10) + */ + public Dfp getLn10() { + return ln10; + } + + /** + * Breaks a string representation up into two {@link Dfp}'s. The split is such that the sum of + * them is equivalent to the input string, but has higher precision than using a single Dfp. + * + * @param a string representation of the number to split + * @return an array of two {@link Dfp Dfp} instances which sum equals a + */ + private Dfp[] split(final String a) { + Dfp result[] = new Dfp[2]; + boolean leading = true; + int sp = 0; + int sig = 0; + + char[] buf = new char[a.length()]; + + for (int i = 0; i < buf.length; i++) { + buf[i] = a.charAt(i); + + if (buf[i] >= '1' && buf[i] <= '9') { + leading = false; + } + + if (buf[i] == '.') { + sig += (400 - sig) % 4; + leading = false; + } + + if (sig == (radixDigits / 2) * 4) { + sp = i; + break; + } + + if (buf[i] >= '0' && buf[i] <= '9' && !leading) { + sig++; + } + } + + result[0] = new Dfp(this, new String(buf, 0, sp)); + + for (int i = 0; i < buf.length; i++) { + buf[i] = a.charAt(i); + if (buf[i] >= '0' && buf[i] <= '9' && i < sp) { + buf[i] = '0'; + } + } + + result[1] = new Dfp(this, new String(buf)); + + return result; + } + + /** + * Recompute the high precision string constants. + * + * @param highPrecisionDecimalDigits precision at which the string constants mus be computed + */ + private static void computeStringConstants(final int highPrecisionDecimalDigits) { + if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) { + + // recompute the string representation of the transcendental constants + final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false); + final Dfp highPrecisionOne = new Dfp(highPrecisionField, 1); + final Dfp highPrecisionTwo = new Dfp(highPrecisionField, 2); + final Dfp highPrecisionThree = new Dfp(highPrecisionField, 3); + + final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt(); + sqr2String = highPrecisionSqr2.toString(); + sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString(); + + final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt(); + sqr3String = highPrecisionSqr3.toString(); + sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString(); + + piString = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString(); + eString = computeExp(highPrecisionOne, highPrecisionOne).toString(); + ln2String = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString(); + ln5String = + computeLn(new Dfp(highPrecisionField, 5), highPrecisionOne, highPrecisionTwo) + .toString(); + ln10String = + computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo) + .toString(); + } + } + + /** + * Compute π using Jonathan and Peter Borwein quartic formula. + * + * @param one constant with value 1 at desired precision + * @param two constant with value 2 at desired precision + * @param three constant with value 3 at desired precision + * @return π + */ + private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) { + + Dfp sqrt2 = two.sqrt(); + Dfp yk = sqrt2.subtract(one); + Dfp four = two.add(two); + Dfp two2kp3 = two; + Dfp ak = two.multiply(three.subtract(two.multiply(sqrt2))); + + // The formula converges quartically. This means the number of correct + // digits is multiplied by 4 at each iteration! Five iterations are + // sufficient for about 160 digits, eight iterations give about + // 10000 digits (this has been checked) and 20 iterations more than + // 160 billions of digits (this has NOT been checked). + // So the limit here is considered sufficient for most purposes ... + for (int i = 1; i < 20; i++) { + final Dfp ykM1 = yk; + + final Dfp y2 = yk.multiply(yk); + final Dfp oneMinusY4 = one.subtract(y2.multiply(y2)); + final Dfp s = oneMinusY4.sqrt().sqrt(); + yk = one.subtract(s).divide(one.add(s)); + + two2kp3 = two2kp3.multiply(four); + + final Dfp p = one.add(yk); + final Dfp p2 = p.multiply(p); + ak = + ak.multiply(p2.multiply(p2)) + .subtract( + two2kp3.multiply(yk) + .multiply(one.add(yk).add(yk.multiply(yk)))); + + if (yk.equals(ykM1)) { + break; + } + } + + return one.divide(ak); + } + + /** + * Compute exp(a). + * + * @param a number for which we want the exponential + * @param one constant with value 1 at desired precision + * @return exp(a) + */ + public static Dfp computeExp(final Dfp a, final Dfp one) { + + Dfp y = new Dfp(one); + Dfp py = new Dfp(one); + Dfp f = new Dfp(one); + Dfp fi = new Dfp(one); + Dfp x = new Dfp(one); + + for (int i = 0; i < 10000; i++) { + x = x.multiply(a); + y = y.add(x.divide(f)); + fi = fi.add(one); + f = f.multiply(fi); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + return y; + } + + /** + * Compute ln(a). + * + * <p>Let f(x) = ln(x), + * + * <p>We know that f'(x) = 1/x, thus from Taylor's theorem we have: + * + * <p>----- n+1 n f(x) = \ (-1) (x - 1) / ---------------- for 1 <= n <= infinity ----- n + * + * <p>or 2 3 4 (x-1) (x-1) (x-1) ln(x) = (x-1) - ----- + ------ - ------ + ... 2 3 4 + * + * <p>alternatively, + * + * <p>2 3 4 x x x ln(x+1) = x - - + - - - + ... 2 3 4 + * + * <p>This series can be used to compute ln(x), but it converges too slowly. + * + * <p>If we substitute -x for x above, we get + * + * <p>2 3 4 x x x ln(1-x) = -x - - - - - - + ... 2 3 4 + * + * <p>Note that all terms are now negative. Because the even powered ones absorbed the sign. + * Now, subtract the series above from the previous one to get ln(x+1) - ln(1-x). Note the even + * terms cancel out leaving only the odd ones + * + * <p>3 5 7 2x 2x 2x ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... 3 5 7 + * + * <p>By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: + * + * <p>3 5 7 x+1 / x x x \ ln ----- = 2 * | x + ---- + ---- + ---- + ... | x-1 \ 3 5 7 / + * + * <p>But now we want to find ln(a), so we need to find the value of x such that a = + * (x+1)/(x-1). This is easily solved to find that x = (a-1)/(a+1). + * + * @param a number for which we want the exponential + * @param one constant with value 1 at desired precision + * @param two constant with value 2 at desired precision + * @return ln(a) + */ + public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) { + + int den = 1; + Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one)); + + Dfp y = new Dfp(x); + Dfp num = new Dfp(x); + Dfp py = new Dfp(y); + for (int i = 0; i < 10000; i++) { + num = num.multiply(x); + num = num.multiply(x); + den += 2; + Dfp t = num.divide(den); + y = y.add(t); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + return y.multiply(two); + } +} diff --git a/src/main/java/org/apache/commons/math3/dfp/DfpMath.java b/src/main/java/org/apache/commons/math3/dfp/DfpMath.java new file mode 100644 index 0000000..ba50d05 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/DfpMath.java @@ -0,0 +1,970 @@ +/* + * 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; + +/** + * Mathematical routines for use with {@link Dfp}. The constants are defined in {@link DfpField} + * + * @since 2.2 + */ +public class DfpMath { + + /** Name for traps triggered by pow. */ + private static final String POW_TRAP = "pow"; + + /** Private Constructor. */ + private DfpMath() {} + + /** + * Breaks a string representation up into two dfp's. + * + * <p>The two dfp are such that the sum of them is equivalent to the input string, but has + * higher precision than using a single dfp. This is useful for improving accuracy of + * exponentiation and critical multiplies. + * + * @param field field to which the Dfp must belong + * @param a string representation to split + * @return an array of two {@link Dfp} which sum is a + */ + protected static Dfp[] split(final DfpField field, final String a) { + Dfp result[] = new Dfp[2]; + char[] buf; + boolean leading = true; + int sp = 0; + int sig = 0; + + buf = new char[a.length()]; + + for (int i = 0; i < buf.length; i++) { + buf[i] = a.charAt(i); + + if (buf[i] >= '1' && buf[i] <= '9') { + leading = false; + } + + if (buf[i] == '.') { + sig += (400 - sig) % 4; + leading = false; + } + + if (sig == (field.getRadixDigits() / 2) * 4) { + sp = i; + break; + } + + if (buf[i] >= '0' && buf[i] <= '9' && !leading) { + sig++; + } + } + + result[0] = field.newDfp(new String(buf, 0, sp)); + + for (int i = 0; i < buf.length; i++) { + buf[i] = a.charAt(i); + if (buf[i] >= '0' && buf[i] <= '9' && i < sp) { + buf[i] = '0'; + } + } + + result[1] = field.newDfp(new String(buf)); + + return result; + } + + /** + * Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link + * Dfp}. + * + * @param a number to split + * @return two elements array containing the split number + */ + protected static Dfp[] split(final Dfp a) { + final Dfp[] result = new Dfp[2]; + final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2)); + result[0] = a.add(shift).subtract(shift); + result[1] = a.subtract(result[0]); + return result; + } + + /** + * Multiply two numbers that are split in to two pieces that are meant to be added together. Use + * binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1 Store the first term in + * result0, the rest in result1 + * + * @param a first factor of the multiplication, in split form + * @param b second factor of the multiplication, in split form + * @return a × b, in split form + */ + protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) { + final Dfp[] result = new Dfp[2]; + + result[1] = a[0].getZero(); + result[0] = a[0].multiply(b[0]); + + /* If result[0] is infinite or zero, don't compute result[1]. + * Attempting to do so may produce NaNs. + */ + + if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) { + return result; + } + + result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1])); + + return result; + } + + /** + * Divide two numbers that are split in to two pieces that are meant to be added together. + * Inverse of split multiply above: (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) ) + * + * @param a dividend, in split form + * @param b divisor, in split form + * @return a / b, in split form + */ + protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) { + final Dfp[] result; + + result = new Dfp[2]; + + result[0] = a[0].divide(b[0]); + result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1])); + result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1]))); + + return result; + } + + /** + * Raise a split base to the a power. + * + * @param base number to raise + * @param a power + * @return base<sup>a</sup> + */ + protected static Dfp splitPow(final Dfp[] base, int a) { + boolean invert = false; + + Dfp[] r = new Dfp[2]; + + Dfp[] result = new Dfp[2]; + result[0] = base[0].getOne(); + result[1] = base[0].getZero(); + + if (a == 0) { + // Special case a = 0 + return result[0].add(result[1]); + } + + if (a < 0) { + // If a is less than zero + invert = true; + a = -a; + } + + // Exponentiate by successive squaring + do { + r[0] = new Dfp(base[0]); + r[1] = new Dfp(base[1]); + int trial = 1; + + int prevtrial; + while (true) { + prevtrial = trial; + trial *= 2; + if (trial > a) { + break; + } + r = splitMult(r, r); + } + + trial = prevtrial; + + a -= trial; + result = splitMult(result, r); + + } while (a >= 1); + + result[0] = result[0].add(result[1]); + + if (invert) { + result[0] = base[0].getOne().divide(result[0]); + } + + return result[0]; + } + + /** + * Raises base to the power a by successive squaring. + * + * @param base number to raise + * @param a power + * @return base<sup>a</sup> + */ + public static Dfp pow(Dfp base, int a) { + boolean invert = false; + + Dfp result = base.getOne(); + + if (a == 0) { + // Special case + return result; + } + + if (a < 0) { + invert = true; + a = -a; + } + + // Exponentiate by successive squaring + do { + Dfp r = new Dfp(base); + Dfp prevr; + int trial = 1; + int prevtrial; + + do { + prevr = new Dfp(r); + prevtrial = trial; + r = r.multiply(r); + trial *= 2; + } while (a > trial); + + r = prevr; + trial = prevtrial; + + a -= trial; + result = result.multiply(r); + + } while (a >= 1); + + if (invert) { + result = base.getOne().divide(result); + } + + return base.newInstance(result); + } + + /** + * Computes e to the given power. a is broken into two parts, such that a = n+m where n is an + * integer. We use pow() to compute e<sup>n</sup> and a Taylor series to compute e<sup>m</sup>. + * We return e*<sup>n</sup> × e<sup>m</sup> + * + * @param a power at which e should be raised + * @return e<sup>a</sup> + */ + public static Dfp exp(final Dfp a) { + + final Dfp inta = a.rint(); + final Dfp fraca = a.subtract(inta); + + final int ia = inta.intValue(); + if (ia > 2147483646) { + // return +Infinity + return a.newInstance((byte) 1, Dfp.INFINITE); + } + + if (ia < -2147483646) { + // return 0; + return a.newInstance(); + } + + final Dfp einta = splitPow(a.getField().getESplit(), ia); + final Dfp efraca = expInternal(fraca); + + return einta.multiply(efraca); + } + + /** + * Computes e to the given power. Where -1 < a < 1. Use the classic Taylor series. 1 + x**2/2! + + * x**3/3! + x**4/4! ... + * + * @param a power at which e should be raised + * @return e<sup>a</sup> + */ + protected static Dfp expInternal(final Dfp a) { + Dfp y = a.getOne(); + Dfp x = a.getOne(); + Dfp fact = a.getOne(); + Dfp py = new Dfp(y); + + for (int i = 1; i < 90; i++) { + x = x.multiply(a); + fact = fact.divide(i); + y = y.add(x.multiply(fact)); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + return y; + } + + /** + * Returns the natural logarithm of a. a is first split into three parts such that a = + * (10000^h)(2^j)k. ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k) k is in the range + * 2/3 < k <4/3 and is passed on to a series expansion. + * + * @param a number from which logarithm is requested + * @return log(a) + */ + public static Dfp log(Dfp a) { + int lr; + Dfp x; + int ix; + int p2 = 0; + + // Check the arguments somewhat here + if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) { + // negative, zero or NaN + a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte) 1, Dfp.QNAN)); + } + + if (a.classify() == Dfp.INFINITE) { + return a; + } + + x = new Dfp(a); + lr = x.log10K(); + + x = x.divide(pow(a.newInstance(10000), lr)); /* This puts x in the range 0-10000 */ + ix = x.floor().intValue(); + + while (ix > 2) { + ix >>= 1; + p2++; + } + + Dfp[] spx = split(x); + Dfp[] spy = new Dfp[2]; + spy[0] = pow(a.getTwo(), p2); // use spy[0] temporarily as a divisor + spx[0] = spx[0].divide(spy[0]); + spx[1] = spx[1].divide(spy[0]); + + spy[0] = a.newInstance("1.33333"); // Use spy[0] for comparison + while (spx[0].add(spx[1]).greaterThan(spy[0])) { + spx[0] = spx[0].divide(2); + spx[1] = spx[1].divide(2); + p2++; + } + + // X is now in the range of 2/3 < x < 4/3 + Dfp[] spz = logInternal(spx); + + spx[0] = a.newInstance(new StringBuilder().append(p2 + 4 * lr).toString()); + spx[1] = a.getZero(); + spy = splitMult(a.getField().getLn2Split(), spx); + + spz[0] = spz[0].add(spy[0]); + spz[1] = spz[1].add(spy[1]); + + spx[0] = a.newInstance(new StringBuilder().append(4 * lr).toString()); + spx[1] = a.getZero(); + spy = splitMult(a.getField().getLn5Split(), spx); + + spz[0] = spz[0].add(spy[0]); + spz[1] = spz[1].add(spy[1]); + + return a.newInstance(spz[0].add(spz[1])); + } + + /** + * Computes the natural log of a number between 0 and 2. Let f(x) = ln(x), + * + * <p>We know that f'(x) = 1/x, thus from Taylor's theorum we have: + * + * <p>----- n+1 n f(x) = \ (-1) (x - 1) / ---------------- for 1 <= n <= infinity ----- n + * + * <p>or 2 3 4 (x-1) (x-1) (x-1) ln(x) = (x-1) - ----- + ------ - ------ + ... 2 3 4 + * + * <p>alternatively, + * + * <p>2 3 4 x x x ln(x+1) = x - - + - - - + ... 2 3 4 + * + * <p>This series can be used to compute ln(x), but it converges too slowly. + * + * <p>If we substitute -x for x above, we get + * + * <p>2 3 4 x x x ln(1-x) = -x - - - - - - + ... 2 3 4 + * + * <p>Note that all terms are now negative. Because the even powered ones absorbed the sign. + * Now, subtract the series above from the previous one to get ln(x+1) - ln(1-x). Note the even + * terms cancel out leaving only the odd ones + * + * <p>3 5 7 2x 2x 2x ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ... 3 5 7 + * + * <p>By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have: + * + * <p>3 5 7 x+1 / x x x \ ln ----- = 2 * | x + ---- + ---- + ---- + ... | x-1 \ 3 5 7 / + * + * <p>But now we want to find ln(a), so we need to find the value of x such that a = + * (x+1)/(x-1). This is easily solved to find that x = (a-1)/(a+1). + * + * @param a number from which logarithm is requested, in split form + * @return log(a) + */ + protected static Dfp[] logInternal(final Dfp a[]) { + + /* Now we want to compute x = (a-1)/(a+1) but this is prone to + * loss of precision. So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4) + */ + Dfp t = a[0].divide(4).add(a[1].divide(4)); + Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25"))); + + Dfp y = new Dfp(x); + Dfp num = new Dfp(x); + Dfp py = new Dfp(y); + int den = 1; + for (int i = 0; i < 10000; i++) { + num = num.multiply(x); + num = num.multiply(x); + den += 2; + t = num.divide(den); + y = y.add(t); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + y = y.multiply(a[0].getTwo()); + + return split(y); + } + + /** + * Computes x to the y power. + * + * <p>Uses the following method: + * + * <p> + * + * <ol> + * <li>Set u = rint(y), v = y-u + * <li>Compute a = v * ln(x) + * <li>Compute b = rint( a/ln(2) ) + * <li>Compute c = a - b*ln(2) + * <li>x<sup>y</sup> = x<sup>u</sup> * 2<sup>b</sup> * e<sup>c</sup> + * </ol> + * + * if |y| > 1e8, then we compute by exp(y*ln(x)) + * + * <p><b>Special Cases</b> + * + * <p> + * + * <ul> + * <li>if y is 0.0 or -0.0 then result is 1.0 + * <li>if y is 1.0 then result is x + * <li>if y is NaN then result is NaN + * <li>if x is NaN and y is not zero then result is NaN + * <li>if |x| > 1.0 and y is +Infinity then result is +Infinity + * <li>if |x| < 1.0 and y is -Infinity then result is +Infinity + * <li>if |x| > 1.0 and y is -Infinity then result is +0 + * <li>if |x| < 1.0 and y is +Infinity then result is +0 + * <li>if |x| = 1.0 and y is +/-Infinity then result is NaN + * <li>if x = +0 and y > 0 then result is +0 + * <li>if x = +Inf and y < 0 then result is +0 + * <li>if x = +0 and y < 0 then result is +Inf + * <li>if x = +Inf and y > 0 then result is +Inf + * <li>if x = -0 and y > 0, finite, not odd integer then result is +0 + * <li>if x = -0 and y < 0, finite, and odd integer then result is -Inf + * <li>if x = -Inf and y > 0, finite, and odd integer then result is -Inf + * <li>if x = -0 and y < 0, not finite odd integer then result is +Inf + * <li>if x = -Inf and y > 0, not finite odd integer then result is +Inf + * <li>if x < 0 and y > 0, finite, and odd integer then result is -(|x|<sup>y</sup>) + * <li>if x < 0 and y > 0, finite, and not integer then result is NaN + * </ul> + * + * @param x base to be raised + * @param y power to which base should be raised + * @return x<sup>y</sup> + */ + public static Dfp pow(Dfp x, final Dfp y) { + + // make sure we don't mix number with different precision + if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) { + x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + final Dfp result = x.newInstance(x.getZero()); + result.nans = Dfp.QNAN; + return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result); + } + + final Dfp zero = x.getZero(); + final Dfp one = x.getOne(); + final Dfp two = x.getTwo(); + boolean invert = false; + int ui; + + /* Check for special cases */ + if (y.equals(zero)) { + return x.newInstance(one); + } + + if (y.equals(one)) { + if (x.isNaN()) { + // Test for NaNs + x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x); + } + return x; + } + + if (x.isNaN() || y.isNaN()) { + // Test for NaNs + x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte) 1, Dfp.QNAN)); + } + + // X == 0 + if (x.equals(zero)) { + if (Dfp.copysign(one, x).greaterThan(zero)) { + // X == +0 + if (y.greaterThan(zero)) { + return x.newInstance(zero); + } else { + return x.newInstance(x.newInstance((byte) 1, Dfp.INFINITE)); + } + } else { + // X == -0 + if (y.classify() == Dfp.FINITE + && y.rint().equals(y) + && !y.remainder(two).equals(zero)) { + // If y is odd integer + if (y.greaterThan(zero)) { + return x.newInstance(zero.negate()); + } else { + return x.newInstance(x.newInstance((byte) -1, Dfp.INFINITE)); + } + } else { + // Y is not odd integer + if (y.greaterThan(zero)) { + return x.newInstance(zero); + } else { + return x.newInstance(x.newInstance((byte) 1, Dfp.INFINITE)); + } + } + } + } + + if (x.lessThan(zero)) { + // Make x positive, but keep track of it + x = x.negate(); + invert = true; + } + + if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) { + if (y.greaterThan(zero)) { + return y; + } else { + return x.newInstance(zero); + } + } + + if (x.lessThan(one) && y.classify() == Dfp.INFINITE) { + if (y.greaterThan(zero)) { + return x.newInstance(zero); + } else { + return x.newInstance(Dfp.copysign(y, one)); + } + } + + if (x.equals(one) && y.classify() == Dfp.INFINITE) { + x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte) 1, Dfp.QNAN)); + } + + if (x.classify() == Dfp.INFINITE) { + // x = +/- inf + if (invert) { + // negative infinity + if (y.classify() == Dfp.FINITE + && y.rint().equals(y) + && !y.remainder(two).equals(zero)) { + // If y is odd integer + if (y.greaterThan(zero)) { + return x.newInstance(x.newInstance((byte) -1, Dfp.INFINITE)); + } else { + return x.newInstance(zero.negate()); + } + } else { + // Y is not odd integer + if (y.greaterThan(zero)) { + return x.newInstance(x.newInstance((byte) 1, Dfp.INFINITE)); + } else { + return x.newInstance(zero); + } + } + } else { + // positive infinity + if (y.greaterThan(zero)) { + return x; + } else { + return x.newInstance(zero); + } + } + } + + if (invert && !y.rint().equals(y)) { + x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID); + return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte) 1, Dfp.QNAN)); + } + + // End special cases + + Dfp r; + if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) { + final Dfp u = y.rint(); + ui = u.intValue(); + + final Dfp v = y.subtract(u); + + if (v.unequal(zero)) { + final Dfp a = v.multiply(log(x)); + final Dfp b = a.divide(x.getField().getLn2()).rint(); + + final Dfp c = a.subtract(b.multiply(x.getField().getLn2())); + r = splitPow(split(x), ui); + r = r.multiply(pow(two, b.intValue())); + r = r.multiply(exp(c)); + } else { + r = splitPow(split(x), ui); + } + } else { + // very large exponent. |y| > 1e8 + r = exp(log(x).multiply(y)); + } + + if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) { + // if y is odd integer + r = r.negate(); + } + + return x.newInstance(r); + } + + /** + * Computes sin(a) Used when 0 < a < pi/4. Uses the classic Taylor series. x - x**3/3! + x**5/5! + * ... + * + * @param a number from which sine is desired, in split form + * @return sin(a) + */ + protected static Dfp sinInternal(Dfp a[]) { + + Dfp c = a[0].add(a[1]); + Dfp y = c; + c = c.multiply(c); + Dfp x = y; + Dfp fact = a[0].getOne(); + Dfp py = new Dfp(y); + + for (int i = 3; i < 90; i += 2) { + x = x.multiply(c); + x = x.negate(); + + fact = fact.divide((i - 1) * i); // 1 over fact + y = y.add(x.multiply(fact)); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + return y; + } + + /** + * Computes cos(a) Used when 0 < a < pi/4. Uses the classic Taylor series for cosine. 1 - + * x**2/2! + x**4/4! ... + * + * @param a number from which cosine is desired, in split form + * @return cos(a) + */ + protected static Dfp cosInternal(Dfp a[]) { + final Dfp one = a[0].getOne(); + + Dfp x = one; + Dfp y = one; + Dfp c = a[0].add(a[1]); + c = c.multiply(c); + + Dfp fact = one; + Dfp py = new Dfp(y); + + for (int i = 2; i < 90; i += 2) { + x = x.multiply(c); + x = x.negate(); + + fact = fact.divide((i - 1) * i); // 1 over fact + + y = y.add(x.multiply(fact)); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + return y; + } + + /** + * computes the sine of the argument. + * + * @param a number from which sine is desired + * @return sin(a) + */ + public static Dfp sin(final Dfp a) { + final Dfp pi = a.getField().getPi(); + final Dfp zero = a.getField().getZero(); + boolean neg = false; + + /* First reduce the argument to the range of +/- PI */ + Dfp x = a.remainder(pi.multiply(2)); + + /* if x < 0 then apply identity sin(-x) = -sin(x) */ + /* This puts x in the range 0 < x < PI */ + if (x.lessThan(zero)) { + x = x.negate(); + neg = true; + } + + /* Since sine(x) = sine(pi - x) we can reduce the range to + * 0 < x < pi/2 + */ + + if (x.greaterThan(pi.divide(2))) { + x = pi.subtract(x); + } + + Dfp y; + if (x.lessThan(pi.divide(4))) { + y = sinInternal(split(x)); + } else { + final Dfp c[] = new Dfp[2]; + final Dfp[] piSplit = a.getField().getPiSplit(); + c[0] = piSplit[0].divide(2).subtract(x); + c[1] = piSplit[1].divide(2); + y = cosInternal(c); + } + + if (neg) { + y = y.negate(); + } + + return a.newInstance(y); + } + + /** + * computes the cosine of the argument. + * + * @param a number from which cosine is desired + * @return cos(a) + */ + public static Dfp cos(Dfp a) { + final Dfp pi = a.getField().getPi(); + final Dfp zero = a.getField().getZero(); + boolean neg = false; + + /* First reduce the argument to the range of +/- PI */ + Dfp x = a.remainder(pi.multiply(2)); + + /* if x < 0 then apply identity cos(-x) = cos(x) */ + /* This puts x in the range 0 < x < PI */ + if (x.lessThan(zero)) { + x = x.negate(); + } + + /* Since cos(x) = -cos(pi - x) we can reduce the range to + * 0 < x < pi/2 + */ + + if (x.greaterThan(pi.divide(2))) { + x = pi.subtract(x); + neg = true; + } + + Dfp y; + if (x.lessThan(pi.divide(4))) { + Dfp c[] = new Dfp[2]; + c[0] = x; + c[1] = zero; + + y = cosInternal(c); + } else { + final Dfp c[] = new Dfp[2]; + final Dfp[] piSplit = a.getField().getPiSplit(); + c[0] = piSplit[0].divide(2).subtract(x); + c[1] = piSplit[1].divide(2); + y = sinInternal(c); + } + + if (neg) { + y = y.negate(); + } + + return a.newInstance(y); + } + + /** + * computes the tangent of the argument. + * + * @param a number from which tangent is desired + * @return tan(a) + */ + public static Dfp tan(final Dfp a) { + return sin(a).divide(cos(a)); + } + + /** + * computes the arc-tangent of the argument. + * + * @param a number from which arc-tangent is desired + * @return atan(a) + */ + protected static Dfp atanInternal(final Dfp a) { + + Dfp y = new Dfp(a); + Dfp x = new Dfp(y); + Dfp py = new Dfp(y); + + for (int i = 3; i < 90; i += 2) { + x = x.multiply(a); + x = x.multiply(a); + x = x.negate(); + y = y.add(x.divide(i)); + if (y.equals(py)) { + break; + } + py = new Dfp(y); + } + + return y; + } + + /** + * computes the arc tangent of the argument + * + * <p>Uses the typical taylor series + * + * <p>but may reduce arguments using the following identity tan(x+y) = (tan(x) + tan(y)) / (1 - + * tan(x)*tan(y)) + * + * <p>since tan(PI/8) = sqrt(2)-1, + * + * <p>atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0 + * + * @param a number from which arc-tangent is desired + * @return atan(a) + */ + public static Dfp atan(final Dfp a) { + final Dfp zero = a.getField().getZero(); + final Dfp one = a.getField().getOne(); + final Dfp[] sqr2Split = a.getField().getSqr2Split(); + final Dfp[] piSplit = a.getField().getPiSplit(); + boolean recp = false; + boolean neg = false; + boolean sub = false; + + final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]); + + Dfp x = new Dfp(a); + if (x.lessThan(zero)) { + neg = true; + x = x.negate(); + } + + if (x.greaterThan(one)) { + recp = true; + x = one.divide(x); + } + + if (x.greaterThan(ty)) { + Dfp sty[] = new Dfp[2]; + sub = true; + + sty[0] = sqr2Split[0].subtract(one); + sty[1] = sqr2Split[1]; + + Dfp[] xs = split(x); + + Dfp[] ds = splitMult(xs, sty); + ds[0] = ds[0].add(one); + + xs[0] = xs[0].subtract(sty[0]); + xs[1] = xs[1].subtract(sty[1]); + + xs = splitDiv(xs, ds); + x = xs[0].add(xs[1]); + + // x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty))); + } + + Dfp y = atanInternal(x); + + if (sub) { + y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8)); + } + + if (recp) { + y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2)); + } + + if (neg) { + y = y.negate(); + } + + return a.newInstance(y); + } + + /** + * computes the arc-sine of the argument. + * + * @param a number from which arc-sine is desired + * @return asin(a) + */ + public static Dfp asin(final Dfp a) { + return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt())); + } + + /** + * computes the arc-cosine of the argument. + * + * @param a number from which arc-cosine is desired + * @return acos(a) + */ + public static Dfp acos(Dfp a) { + Dfp result; + boolean negative = false; + + if (a.lessThan(a.getZero())) { + negative = true; + } + + a = Dfp.copysign(a, a.getOne()); // absolute value + + result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a)); + + if (negative) { + result = a.getField().getPi().subtract(result); + } + + return a.newInstance(result); + } +} diff --git a/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java b/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java new file mode 100644 index 0000000..4507b41 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java @@ -0,0 +1,40 @@ +/* + * 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; + +/** + * An interface representing a univariate {@link Dfp} function. + * + * @deprecated as of 3.6, replaced with {@link + * org.apache.commons.math3.analysis.RealFieldUnivariateFunction} + */ +@Deprecated +public interface UnivariateDfpFunction { + + /** + * Compute the value of the function. + * + * @param x Point at which the function value should be computed. + * @return the value. + * @throws IllegalArgumentException when the activated method itself can ascertain that + * preconditions, specified in the API expressed at the level of the activated method, have + * been violated. In the vast majority of cases where Commons-Math throws + * IllegalArgumentException, it is the result of argument checking of actual parameters + * immediately passed to a method. + */ + Dfp value(Dfp x); +} diff --git a/src/main/java/org/apache/commons/math3/dfp/package-info.java b/src/main/java/org/apache/commons/math3/dfp/package-info.java new file mode 100644 index 0000000..f4a7ca8 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/dfp/package-info.java @@ -0,0 +1,80 @@ +/* + * 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. + */ +/** + * 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>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 keep as portable as possible. + * <li>Performance + * <li>Accuracy - Results should always be +/- 1 ULP for basic algebraic operation + * <li>Comply with IEEE 854-1987 as much as possible. (See IEEE 854-1987 notes below) + * </ol> + * + * <p>Trade offs: + * + * <ol> + * <li>Memory foot print. I'm using more memory than necessary to represent numbers to get better + * performance. + * <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. + * </ol> + * + * <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>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>The radix of 10000 was chosen because it should be faster to operate on 4 decimal digits at + * once instead of one at a time. Radix 10 behavior can be realized by add an additional rounding + * step to ensure that the number of decimal digits represented is constant. + * + * <p>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>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>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. + */ +package org.apache.commons.math3.dfp; |