Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added the dfp library 

The Dfp class represent the high precision numbers, it implements our existing Field interface and hence each instance is associated with a DfpField that provides the constants at the required precision as well as factory methods. This allowed to remove the compile-time constraint in the library. Users can even use at the same time a field for 20 decimals digits precision and another field for 100 digits precision. Dfp instances with different precision CANNOT be mixed in the same computation (doing so creates a NaN). A few utility methods have been added, like constructors from integral types, isInfinite and isNaN methods, equal has been renames to equals and its signature changed to match the general Object method (a hashcode method has been added too).
JIRA: MATH-412

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@992697 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2010-09-04 22:59:21 +00:00
parent 18f77ee26c
commit 85d35b0ff6
10 changed files with 7323 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2400
src/main/java/org/apache/commons/math/dfp/Dfp.java Normal file
View File

File diff suppressed because it is too large Load Diff

359
src/main/java/org/apache/commons/math/dfp/DfpDec.java Normal file
View File

@ -0,0 +1,359 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.dfp;
/** 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.
* @version $Revision$ $Date$
* @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} */
public Dfp newInstance() {
return new DfpDec(getField());
}
/** {@inheritDoc} */
public Dfp newInstance(final byte x) {
return new DfpDec(getField(), x);
}
/** {@inheritDoc} */
public Dfp newInstance(final int x) {
return new DfpDec(getField(), x);
}
/** {@inheritDoc} */
public Dfp newInstance(final long x) {
return new DfpDec(getField(), x);
}
/** {@inheritDoc} */
public Dfp newInstance(final double x) {
return new DfpDec(getField(), x);
}
/** {@inheritDoc} */
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} */
public Dfp newInstance(final String s) {
return new DfpDec(getField(), s);
}
/** {@inheritDoc} */
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} */
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} */
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(log10() - 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(log10());
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;
}
}

742
src/main/java/org/apache/commons/math/dfp/DfpField.java Normal file
View File

@ -0,0 +1,742 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.dfp;
import org.apache.commons.math.Field;
/** Field for Decimal floating point instances.
* @version $Revision$ $Date$
* @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 &radic;2. */
private static String sqr2String;
/** High precision string representation of &radic;2 / 2. */
private static String sqr2ReciprocalString;
/** High precision string representation of &radic;3. */
private static String sqr3String;
/** High precision string representation of &radic;3 / 3. */
private static String sqr3ReciprocalString;
/** High precision string representation of &pi;. */
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 &radic;2. */
private final Dfp sqr2;
/** A two elements {@link Dfp} array with value &radic;2 split in two pieces. */
private final Dfp[] sqr2Split;
/** A {@link Dfp} with value &radic;2 / 2. */
private final Dfp sqr2Reciprocal;
/** A {@link Dfp} with value &radic;3. */
private final Dfp sqr3;
/** A {@link Dfp} with value &radic;3 / 3. */
private final Dfp sqr3Reciprocal;
/** A {@link Dfp} with value &pi;. */
private final Dfp pi;
/** A two elements {@link Dfp} array with value &pi; 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.
* </p>
* @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.
* </p>
* @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)
*/
public DfpField(final int decimalDigits, final boolean computeConstants) {
this.radixDigits = (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)}
* </p>
* @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;
}
/** Get the constant 2.
* @return a {@link Dfp} with value 2
*/
public Dfp getTwo() {
return two;
}
/** Get the constant &radic;2.
* @return a {@link Dfp} with value &radic;2
*/
public Dfp getSqr2() {
return sqr2;
}
/** Get the constant &radic;2 split in two pieces.
* @return a {@link Dfp} with value &radic;2 split in two pieces
*/
public Dfp[] getSqr2Split() {
return sqr2Split.clone();
}
/** Get the constant &radic;2 / 2.
* @return a {@link Dfp} with value &radic;2 / 2
*/
public Dfp getSqr2Reciprocal() {
return sqr2Reciprocal;
}
/** Get the constant &radic;3.
* @return a {@link Dfp} with value &radic;3
*/
public Dfp getSqr3() {
return sqr3;
}
/** Get the constant &radic;3 / 3.
* @return a {@link Dfp} with value &radic;3 / 3
*/
public Dfp getSqr3Reciprocal() {
return sqr3Reciprocal;
}
/** Get the constant &pi;.
* @return a {@link Dfp} with value &pi;
*/
public Dfp getPi() {
return pi;
}
/** Get the constant &pi; split in two pieces.
* @return a {@link Dfp} with value &pi; 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 &pi; by atan(1/&radic;(3)) = &pi;/6.
* @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 &pi;
*/
private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {
Dfp x = three;
x = x.sqrt();
x = one.divide(x);
Dfp denom = one;
Dfp py = new Dfp(x);
Dfp y = new Dfp(x);
for (int i = 1; i < 10000; i++) {
x = x.divide(three);
denom = denom.add(two);
if ((i&1) != 0) {
y = y.subtract(x.divide(denom));
} else {
y = y.add(x.divide(denom));
}
if (y.equals(py)) {
break;
}
py = new Dfp(y);
}
return y.multiply(new Dfp(one.getField(), 6));
}
/** 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).
*
* Let f(x) = ln(x),
*
* We know that f'(x) = 1/x, thus from Taylor's theorem we have:
*
* ----- n+1 n
* f(x) = \ (-1) (x - 1)
* / ---------------- for 1 <= n <= infinity
* ----- n
*
* or
* 2 3 4
* (x-1) (x-1) (x-1)
* ln(x) = (x-1) - ----- + ------ - ------ + ...
* 2 3 4
*
* alternatively,
*
* 2 3 4
* x x x
* ln(x+1) = x - - + - - - + ...
* 2 3 4
*
* This series can be used to compute ln(x), but it converges too slowly.
*
* If we substitute -x for x above, we get
*
* 2 3 4
* x x x
* ln(1-x) = -x - - - - - - + ...
* 2 3 4
*
* 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
*
* 3 5 7
* 2x 2x 2x
* ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
* 3 5 7
*
* By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
*
* 3 5 7
* x+1 / x x x \
* ln ----- = 2 * | x + ---- + ---- + ---- + ... |
* x-1 \ 3 5 7 /
*
* 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 = den + 2;
Dfp t = num.divide(den);
y = y.add(t);
if (y.equals(py)) {
break;
}
py = new Dfp(y);
}
return y.multiply(two);
}
}

969
src/main/java/org/apache/commons/math/dfp/DfpMath.java Normal file
View File

@ -0,0 +1,969 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.dfp;
/** Mathematical routines for use with {@link Dfp}.
* The constants are defined in {@link DfpField}
* @version $Revision$ $Date$
* @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 &times; 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 = 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 = trial * 2;
} while (a>trial);
r = prevr;
trial = prevtrial;
a = 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> &times; 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, (byte) 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.equals(a) == false)) {
// negative, zero or NaN
a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, (byte) 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 StringBuffer().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 StringBuffer().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),
*
* We know that f'(x) = 1/x, thus from Taylor's theorum we have:
*
* ----- n+1 n
* f(x) = \ (-1) (x - 1)
* / ---------------- for 1 <= n <= infinity
* ----- n
*
* or
* 2 3 4
* (x-1) (x-1) (x-1)
* ln(x) = (x-1) - ----- + ------ - ------ + ...
* 2 3 4
*
* alternatively,
*
* 2 3 4
* x x x
* ln(x+1) = x - - + - - - + ...
* 2 3 4
*
* This series can be used to compute ln(x), but it converges too slowly.
*
* If we substitute -x for x above, we get
*
* 2 3 4
* x x x
* ln(1-x) = -x - - - - - - + ...
* 2 3 4
*
* 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
*
* 3 5 7
* 2x 2x 2x
* ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
* 3 5 7
*
* By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
*
* 3 5 7
* x+1 / x x x \
* ln ----- = 2 * | x + ---- + ---- + ---- + ... |
* x-1 \ 3 5 7 /
*
* 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 = 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.equals(x)) {
// Test for NaNs
x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
}
return x;
}
if (!x.equals(x) || !y.equals(y)) {
// Test for NaNs
x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, (byte) 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, (byte)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, (byte)Dfp.INFINITE));
}
} else {
// Y is not odd integer
if (y.greaterThan(zero)) {
return x.newInstance(zero);
} else {
return x.newInstance(x.newInstance((byte)1, (byte)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, (byte) 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, (byte)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, (byte)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, (byte) 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) {
// if y is odd integer
if (y.rint().equals(y) && !y.remainder(two).equals(zero)) {
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))) {
Dfp c[] = new Dfp[2];
c[0] = x;
c[1] = zero;
//y = sinInternal(c);
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
*
* Uses the typical taylor series
*
* but may reduce arguments using the following identity
* tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
*
* since tan(PI/8) = sqrt(2)-1,
*
* 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);
}
}

88
src/main/java/org/apache/commons/math/dfp/package.html Normal file
View File

@ -0,0 +1,88 @@
<html>
<!--
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.
-->
<!-- $Revision$ $Date$ -->
<body>
Decimal floating point library for Java
<p>Another floating point class. This one is built using radix 10000
which is 10<sup>4</sup>, so its almost decimal.</p>
<p>The design goals here are:
<ol>
<li>Decimal math, or close to it</li>
<li>Settable precision (but no mix between numbers using different settings)</li>
<li>Portability. Code should be keep as portable as possible.</li>
<li>Performance</li>
<li>Accuracy - Results should always be +/- 1 ULP for basic
algebraic operation</li>
<li>Comply with IEEE 854-1987 as much as possible.
(See IEEE 854-1987 notes below)</li>
</ol></p>
<p>Trade offs:
<ol>
<li>Memory foot print. I'm using more memory than necessary to
represent numbers to get better performance.</li>
<li>Digits are bigger, so rounding is a greater loss. So, if you
really need 12 decimal digits, better use 4 base 10000 digits
there can be one partially filled.</li>
</ol></p>
<p>Numbers are represented in the following form:
<pre>
n = sign &times; mant &times; (radix)<sup>exp</sup>;</p>
</pre>
where sign is &plusmn;1, mantissa represents a fractional number between
zero and one. mant[0] is the least significant digit.
exp is in the range of -32767 to 32768</p>
<p>IEEE 854-1987 Notes and differences</p>
<p>IEEE 854 requires the radix to be either 2 or 10. The radix here is
10000, so that requirement is not met, but it is possible that a
subclassed can be made to make it behave as a radix 10
number. It is my opinion that if it looks and behaves as a radix
10 number then it is one and that requirement would be met.</p>
<p>The radix of 10000 was chosen because it should be faster to operate
on 4 decimal digits at once instead of one at a time. Radix 10 behavior
can be realized by add an additional rounding step to ensure that
the number of decimal digits represented is constant.</p>
<p>The IEEE standard specifically leaves out internal data encoding,
so it is reasonable to conclude that such a subclass of this radix
10000 system is merely an encoding of a radix 10 system.</p>
<p>IEEE 854 also specifies the existence of "sub-normal" numbers. This
class does not contain any such entities. The most significant radix
10000 digit is always non-zero. Instead, we support "gradual underflow"
by raising the underflow flag for numbers less with exponent less than
expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
Thus the smallest number we can represent would be:
1E(-(MIN_EXP-digits-1)*4), eg, for digits=5, MIN_EXP=-32767, that would
be 1e-131092.</p>
<p>IEEE 854 defines that the implied radix point lies just to the right
of the most significant digit and to the left of the remaining digits.
This implementation puts the implied radix point to the left of all
digits including the most significant one. The most significant digit
here is the one just to the right of the radix point. This is a fine
detail and is really only a matter of definition. Any side effects of
this can be rendered invisible by a subclass.</p>
</body>
</html>

7
src/site/xdoc/changes.xml
View File

@ -71,6 +71,13 @@ The <action> type attribute can be add,update,fix,remove.
</action>
</release>
<release version="2.2" date="TBD" description="TBD">
<action dev="luc" type="fix" issue="MATH-412" due-to="Bill Rossi">
Added the dfp library providing arbitrary precision floating point computation in the spirit of
IEEE 854-1987 (not exactly as it uses base 1000 instead of base 10). In addition to finite numbers,
infinities and NaNs are available (but there are no subnormals). All IEEE 854-1987 rounding modes and
signaling flags are supported. The available operations are +, -, *, / and the available functions
are sqrt, sin, cos, tan, asin, acos, atan, exp, log.
</action>
<action dev="luc" type="fix" issue="MATH-375" due-to="Bill Rossi">
Added faster and more accurate version of traditional mathematical functions in a FastMath
class intended to be a drop-in replacement for java.util.Math at source-level. Some functions

90
src/test/java/org/apache/commons/math/dfp/Decimal10.java Normal file
View File

@ -0,0 +1,90 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.dfp;
public class Decimal10 extends DfpDec {
Decimal10(final DfpField factory) {
super(factory);
}
Decimal10(final DfpField factory, final byte x) {
super(factory, x);
}
Decimal10(final DfpField factory, final int x) {
super(factory, x);
}
Decimal10(final DfpField factory, final long x) {
super(factory, x);
}
Decimal10(final DfpField factory, final double x) {
super(factory, x);
}
public Decimal10(final Dfp d) {
super(d);
}
public Decimal10(final DfpField factory, final String s) {
super(factory, s);
}
protected Decimal10(final DfpField factory, final byte sign, final byte nans) {
super(factory, sign, nans);
}
public Dfp newInstance() {
return new Decimal10(getField());
}
public Dfp newInstance(final byte x) {
return new Decimal10(getField(), x);
}
public Dfp newInstance(final int x) {
return new Decimal10(getField(), x);
}
public Dfp newInstance(final long x) {
return new Decimal10(getField(), x);
}
public Dfp newInstance(final double x) {
return new Decimal10(getField(), x);
}
public Dfp newInstance(final Dfp d) {
return new Decimal10(d);
}
public Dfp newInstance(final String s) {
return new Decimal10(getField(), s);
}
public Dfp newInstance(final byte sign, final byte nans) {
return new Decimal10(getField(), sign, nans);
}
protected int getDecimalDigits() {
return 10;
}
}

574
src/test/java/org/apache/commons/math/dfp/DfpDecTest.java Normal file
View File

@ -0,0 +1,574 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.dfp;
import org.junit.After;
import org.junit.Assert;
import org.junit.Before;
import org.junit.Test;
public class DfpDecTest {
private DfpField field;
private Dfp pinf;
private Dfp ninf;
private Dfp nan;
private Dfp snan;
private Dfp qnan;
@Before
public void setUp() {
// Some basic setup. Define some constants and clear the status flags
field = new DfpField(20);
pinf = new DfpDec(field, 1).divide(new DfpDec(field, 0));
ninf = new DfpDec(field, -1).divide(new DfpDec(field, 0));
nan = new DfpDec(field, 0).divide(new DfpDec(field, 0));
snan = field.newDfp((byte)1, Dfp.SNAN);
qnan = field.newDfp((byte)1, Dfp.QNAN);
ninf.getField().clearIEEEFlags();
}
@After
public void tearDown() {
field = null;
pinf = null;
ninf = null;
nan = null;
snan = null;
qnan = null;
}
// Generic test function. Takes params x and y and tests them for
// equality. Then checks the status flags against the flags argument.
// If the test fail, it prints the desc string
private void test(Dfp x, Dfp y, int flags, String desc) {
boolean b = x.equals(y);
if (!x.equals(y) && !x.unequal(y)) // NaNs involved
b = (x.toString().equals(y.toString()));
if (x.equals(new DfpDec(field, 0))) // distinguish +/- zero
b = (b && (x.toString().equals(y.toString())));
b = (b && x.getField().getIEEEFlags() == flags);
if (!b)
Assert.assertTrue("assersion failed "+desc+" x = "+x.toString()+" flags = "+x.getField().getIEEEFlags(), b);
x.getField().clearIEEEFlags();
}
@Test
public void testRound()
{
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_EVEN);
test(new DfpDec(field, "12345678901234567890"),
new DfpDec(field, "12345678901234568000"),
DfpField.FLAG_INEXACT, "Round #1");
test(new DfpDec(field, "0.12345678901234567890"),
new DfpDec(field, "0.12345678901234568"),
DfpField.FLAG_INEXACT, "Round #2");
test(new DfpDec(field, "0.12345678901234567500"),
new DfpDec(field, "0.12345678901234568"),
DfpField.FLAG_INEXACT, "Round #3");
test(new DfpDec(field, "0.12345678901234568500"),
new DfpDec(field, "0.12345678901234568"),
DfpField.FLAG_INEXACT, "Round #4");
test(new DfpDec(field, "0.12345678901234568501"),
new DfpDec(field, "0.12345678901234569"),
DfpField.FLAG_INEXACT, "Round #5");
test(new DfpDec(field, "0.12345678901234568499"),
new DfpDec(field, "0.12345678901234568"),
DfpField.FLAG_INEXACT, "Round #6");
test(new DfpDec(field, "1.2345678901234567890"),
new DfpDec(field, "1.2345678901234568"),
DfpField.FLAG_INEXACT, "Round #7");
test(new DfpDec(field, "1.2345678901234567500"),
new DfpDec(field, "1.2345678901234568"),
DfpField.FLAG_INEXACT, "Round #8");
test(new DfpDec(field, "1.2345678901234568500"),
new DfpDec(field, "1.2345678901234568"),
DfpField.FLAG_INEXACT, "Round #9");
test(new DfpDec(field, "1.2345678901234568000").add(new DfpDec(field, ".0000000000000000501")),
new DfpDec(field, "1.2345678901234569"),
DfpField.FLAG_INEXACT, "Round #10");
test(new DfpDec(field, "1.2345678901234568499"),
new DfpDec(field, "1.2345678901234568"),
DfpField.FLAG_INEXACT, "Round #11");
test(new DfpDec(field, "12.345678901234567890"),
new DfpDec(field, "12.345678901234568"),
DfpField.FLAG_INEXACT, "Round #12");
test(new DfpDec(field, "12.345678901234567500"),
new DfpDec(field, "12.345678901234568"),
DfpField.FLAG_INEXACT, "Round #13");
test(new DfpDec(field, "12.345678901234568500"),
new DfpDec(field, "12.345678901234568"),
DfpField.FLAG_INEXACT, "Round #14");
test(new DfpDec(field, "12.345678901234568").add(new DfpDec(field, ".000000000000000501")),
new DfpDec(field, "12.345678901234569"),
DfpField.FLAG_INEXACT, "Round #15");
test(new DfpDec(field, "12.345678901234568499"),
new DfpDec(field, "12.345678901234568"),
DfpField.FLAG_INEXACT, "Round #16");
test(new DfpDec(field, "123.45678901234567890"),
new DfpDec(field, "123.45678901234568"),
DfpField.FLAG_INEXACT, "Round #17");
test(new DfpDec(field, "123.45678901234567500"),
new DfpDec(field, "123.45678901234568"),
DfpField.FLAG_INEXACT, "Round #18");
test(new DfpDec(field, "123.45678901234568500"),
new DfpDec(field, "123.45678901234568"),
DfpField.FLAG_INEXACT, "Round #19");
test(new DfpDec(field, "123.456789012345685").add(new DfpDec(field, ".00000000000000501")),
new DfpDec(field, "123.45678901234569"),
DfpField.FLAG_INEXACT, "Round #20");
test(new DfpDec(field, "123.45678901234568499"),
new DfpDec(field, "123.45678901234568"),
DfpField.FLAG_INEXACT, "Round #21");
field.setRoundingMode(DfpField.RoundingMode.ROUND_DOWN);
// Round down
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.9")),
new DfpDec(field, "12345678901234567"),
DfpField.FLAG_INEXACT, "Round #22");
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.99999999")),
new DfpDec(field, "12345678901234567"),
DfpField.FLAG_INEXACT, "Round #23");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.99999999")),
new DfpDec(field, "-12345678901234567"),
DfpField.FLAG_INEXACT, "Round #24");
field.setRoundingMode(DfpField.RoundingMode.ROUND_UP);
// Round up
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.1")),
new DfpDec(field, "12345678901234568"),
DfpField.FLAG_INEXACT, "Round #25");
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.0001")),
new DfpDec(field, "12345678901234568"),
DfpField.FLAG_INEXACT, "Round #26");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.1")),
new DfpDec(field, "-12345678901234568"),
DfpField.FLAG_INEXACT, "Round #27");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.0001")),
new DfpDec(field, "-12345678901234568"),
DfpField.FLAG_INEXACT, "Round #28");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "0")),
new DfpDec(field, "-12345678901234567"),
0, "Round #28.5");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_UP);
// Round half up
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.499999999999")),
new DfpDec(field, "12345678901234567"),
DfpField.FLAG_INEXACT, "Round #29");
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.50000001")),
new DfpDec(field, "12345678901234568"),
DfpField.FLAG_INEXACT, "Round #30");
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.5")),
new DfpDec(field, "12345678901234568"),
DfpField.FLAG_INEXACT, "Round #30.5");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.499999999999")),
new DfpDec(field, "-12345678901234567"),
DfpField.FLAG_INEXACT, "Round #31");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.50000001")),
new DfpDec(field, "-12345678901234568"),
DfpField.FLAG_INEXACT, "Round #32");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_DOWN);
// Round half down
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.5001")),
new DfpDec(field, "12345678901234568"),
DfpField.FLAG_INEXACT, "Round #33");
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.5000")),
new DfpDec(field, "12345678901234567"),
DfpField.FLAG_INEXACT, "Round #34");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.5001")),
new DfpDec(field, "-12345678901234568"),
DfpField.FLAG_INEXACT, "Round #35");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.6")),
new DfpDec(field, "-12345678901234568"),
DfpField.FLAG_INEXACT, "Round #35.5");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.5000")),
new DfpDec(field, "-12345678901234567"),
DfpField.FLAG_INEXACT, "Round #36");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_ODD);
// Round half odd
test(new DfpDec(field, "12345678901234568").add(new DfpDec(field, "0.5000")),
new DfpDec(field, "12345678901234569"),
DfpField.FLAG_INEXACT, "Round #37");
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.5000")),
new DfpDec(field, "12345678901234567"),
DfpField.FLAG_INEXACT, "Round #38");
test(new DfpDec(field, "-12345678901234568").add(new DfpDec(field, "-0.5000")),
new DfpDec(field, "-12345678901234569"),
DfpField.FLAG_INEXACT, "Round #39");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.5000")),
new DfpDec(field, "-12345678901234567"),
DfpField.FLAG_INEXACT, "Round #40");
field.setRoundingMode(DfpField.RoundingMode.ROUND_CEIL);
// Round ceil
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.0001")),
new DfpDec(field, "12345678901234568"),
DfpField.FLAG_INEXACT, "Round #41");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.9999")),
new DfpDec(field, "-12345678901234567"),
DfpField.FLAG_INEXACT, "Round #42");
field.setRoundingMode(DfpField.RoundingMode.ROUND_FLOOR);
// Round floor
test(new DfpDec(field, "12345678901234567").add(new DfpDec(field, "0.9999")),
new DfpDec(field, "12345678901234567"),
DfpField.FLAG_INEXACT, "Round #43");
test(new DfpDec(field, "-12345678901234567").add(new DfpDec(field, "-0.0001")),
new DfpDec(field, "-12345678901234568"),
DfpField.FLAG_INEXACT, "Round #44");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_EVEN); // reset
}
@Test
public void testRoundDecimal10()
{
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_EVEN);
test(new Decimal10(field, "1234567891234567890"),
new Decimal10(field, "1234567891000000000"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #1");
test(new Decimal10(field, "0.1234567891634567890"),
new Decimal10(field, "0.1234567892"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #2");
test(new Decimal10(field, "0.1234567891500000000"),
new Decimal10(field, "0.1234567892"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #3");
test(new Decimal10(field, "0.1234567890500"),
new Decimal10(field, "0.1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #4");
test(new Decimal10(field, "0.1234567890501"),
new Decimal10(field, "0.1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #5");
test(new Decimal10(field, "0.1234567890499"),
new Decimal10(field, "0.1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #6");
test(new Decimal10(field, "1.234567890890"),
new Decimal10(field, "1.234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #7");
test(new Decimal10(field, "1.234567891500"),
new Decimal10(field, "1.234567892"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #8");
test(new Decimal10(field, "1.234567890500"),
new Decimal10(field, "1.234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #9");
test(new Decimal10(field, "1.234567890000").add(new Decimal10(field, ".000000000501")),
new Decimal10(field, "1.234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #10");
test(new Decimal10(field, "1.234567890499"),
new Decimal10(field, "1.234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #11");
test(new Decimal10(field, "12.34567890890"),
new Decimal10(field, "12.34567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #12");
test(new Decimal10(field, "12.34567891500"),
new Decimal10(field, "12.34567892"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #13");
test(new Decimal10(field, "12.34567890500"),
new Decimal10(field, "12.34567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #14");
test(new Decimal10(field, "12.34567890").add(new Decimal10(field, ".00000000501")),
new Decimal10(field, "12.34567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #15");
test(new Decimal10(field, "12.34567890499"),
new Decimal10(field, "12.34567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #16");
test(new Decimal10(field, "123.4567890890"),
new Decimal10(field, "123.4567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #17");
test(new Decimal10(field, "123.4567891500"),
new Decimal10(field, "123.4567892"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #18");
test(new Decimal10(field, "123.4567890500"),
new Decimal10(field, "123.4567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #19");
test(new Decimal10(field, "123.4567890").add(new Decimal10(field, ".0000000501")),
new Decimal10(field, "123.4567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #20");
test(new Decimal10(field, "123.4567890499"),
new Decimal10(field, "123.4567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #21");
field.setRoundingMode(DfpField.RoundingMode.ROUND_DOWN);
// RoundDecimal10 down
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.9")),
new Decimal10(field, "1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #22");
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.99999999")),
new Decimal10(field, "1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #23");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.99999999")),
new Decimal10(field, "-1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #24");
field.setRoundingMode(DfpField.RoundingMode.ROUND_UP);
// RoundDecimal10 up
test(new Decimal10(field, 1234567890).add(new Decimal10(field, "0.1")),
new Decimal10(field, 1234567891l),
DfpField.FLAG_INEXACT, "RoundDecimal10 #25");
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.0001")),
new Decimal10(field, "1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #26");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.1")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #27");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.0001")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #28");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "0")),
new Decimal10(field, "-1234567890"),
0, "RoundDecimal10 #28.5");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_UP);
// RoundDecimal10 half up
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.4999999999")),
new Decimal10(field, "1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #29");
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.50000001")),
new Decimal10(field, "1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #30");
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.5")),
new Decimal10(field, "1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #30.5");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.4999999999")),
new Decimal10(field, "-1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #31");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.50000001")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #32");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_DOWN);
// RoundDecimal10 half down
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.5001")),
new Decimal10(field, "1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #33");
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.5000")),
new Decimal10(field, "1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #34");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.5001")),
new Decimal10(field, "-1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #35");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.6")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #35.5");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.5000")),
new Decimal10(field, "-1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #36");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_ODD);
// RoundDecimal10 half odd
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.5000")),
new Decimal10(field, "1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #37");
test(new Decimal10(field, "1234567891").add(new Decimal10(field, "0.5000")),
new Decimal10(field, "1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #38");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.5000")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #39");
test(new Decimal10(field, "-1234567891").add(new Decimal10(field, "-0.5000")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #40");
field.setRoundingMode(DfpField.RoundingMode.ROUND_CEIL);
// RoundDecimal10 ceil
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.0001")),
new Decimal10(field, "1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #41");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.9999")),
new Decimal10(field, "-1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #42");
field.setRoundingMode(DfpField.RoundingMode.ROUND_FLOOR);
// RoundDecimal10 floor
test(new Decimal10(field, "1234567890").add(new Decimal10(field, "0.9999")),
new Decimal10(field, "1234567890"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #43");
test(new Decimal10(field, "-1234567890").add(new Decimal10(field, "-0.0001")),
new Decimal10(field, "-1234567891"),
DfpField.FLAG_INEXACT, "RoundDecimal10 #44");
field.setRoundingMode(DfpField.RoundingMode.ROUND_HALF_EVEN); // reset
}
@Test
public void testNextAfter()
{
test(new DfpDec(field, 1).nextAfter(pinf),
new DfpDec(field, "1.0000000000000001"),
0, "NextAfter #1");
test(new DfpDec(field, "1.0000000000000001").nextAfter(ninf),
new DfpDec(field, 1),
0, "NextAfter #1.5");
test(new DfpDec(field, 1).nextAfter(ninf),
new DfpDec(field, "0.99999999999999999"),
0, "NextAfter #2");
test(new DfpDec(field, "0.99999999999999999").nextAfter(new DfpDec(field, 2)),
new DfpDec(field, 1),
0, "NextAfter #3");
test(new DfpDec(field, -1).nextAfter(ninf),
new DfpDec(field, "-1.0000000000000001"),
0, "NextAfter #4");
test(new DfpDec(field, -1).nextAfter(pinf),
new DfpDec(field, "-0.99999999999999999"),
0, "NextAfter #5");
test(new DfpDec(field, "-0.99999999999999999").nextAfter(new DfpDec(field, -2)),
new DfpDec(field, (byte) -1),
0, "NextAfter #6");
test(new DfpDec(field, (byte) 2).nextAfter(new DfpDec(field, 2)),
new DfpDec(field, 2l),
0, "NextAfter #7");
test(new DfpDec(field, 0).nextAfter(new DfpDec(field, 0)),
new DfpDec(field, 0),
0, "NextAfter #8");
test(new DfpDec(field, -2).nextAfter(new DfpDec(field, -2)),
new DfpDec(field, -2),
0, "NextAfter #9");
test(new DfpDec(field, 0).nextAfter(new DfpDec(field, 1)),
new DfpDec(field, "1e-131092"),
DfpField.FLAG_UNDERFLOW, "NextAfter #10");
test(new DfpDec(field, 0).nextAfter(new DfpDec(field, -1)),
new DfpDec(field, "-1e-131092"),
DfpField.FLAG_UNDERFLOW, "NextAfter #11");
test(new DfpDec(field, "-1e-131092").nextAfter(pinf),
new DfpDec(field, "-0"),
DfpField.FLAG_UNDERFLOW|DfpField.FLAG_INEXACT, "Next After #12");
test(new DfpDec(field, "1e-131092").nextAfter(ninf),
new DfpDec(field, "0"),
DfpField.FLAG_UNDERFLOW|DfpField.FLAG_INEXACT, "Next After #13");
test(new DfpDec(field, "9.9999999999999999e131078").nextAfter(pinf),
pinf,
DfpField.FLAG_OVERFLOW|DfpField.FLAG_INEXACT, "Next After #14");
}
}

587
src/test/java/org/apache/commons/math/dfp/DfpMathTest.java Normal file
View File

@ -0,0 +1,587 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.dfp;
import org.junit.After;
import org.junit.Assert;
import org.junit.Before;
public class DfpMathTest {
private DfpField factory;
private Dfp pinf;
private Dfp ninf;
private Dfp nan;
private Dfp snan;
private Dfp qnan;
@Before
public void setUp() {
// Some basic setup. Define some constants and clear the status flags
factory = new DfpField(20);
pinf = factory.newDfp("1").divide(factory.newDfp("0"));
ninf = factory.newDfp("-1").divide(factory.newDfp("0"));
nan = factory.newDfp("0").divide(factory.newDfp("0"));
snan = factory.newDfp((byte)1, Dfp.SNAN);
qnan = factory.newDfp((byte)1, Dfp.QNAN);
ninf.getField().clearIEEEFlags();
// force loading of dfpmath
Dfp pi = factory.getPi();
pi.getField().clearIEEEFlags();
}
@After
public void tearDown() {
pinf = null;
ninf = null;
nan = null;
snan = null;
qnan = null;
}
// Generic test function. Takes params x and y and tests them for
// equality. Then checks the status flags against the flags argument.
// If the test fail, it prints the desc string
private void test(Dfp x, Dfp y, int flags, String desc)
{
boolean b = x.equals(y);
if (!x.equals(y) && !x.unequal(y)) // NaNs involved
b = (x.toString().equals(y.toString()));
if (x.equals(factory.newDfp("0"))) // distinguish +/- zero
b = (b && (x.toString().equals(y.toString())));
b = (b && x.getField().getIEEEFlags() == flags);
if (!b)
Assert.assertTrue("assersion failed "+desc+" x = "+x.toString()+" flags = "+x.getField().getIEEEFlags(), b);
x.getField().clearIEEEFlags();
}
public void testPow()
{
// Test special cases exponent of zero
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("0")),
factory.newDfp("1"),
0, "pow #1");
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("-0")),
factory.newDfp("1"),
0, "pow #2");
test(DfpMath.pow(factory.newDfp("2"), factory.newDfp("0")),
factory.newDfp("1"),
0, "pow #3");
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("-0")),
factory.newDfp("1"),
0, "pow #4");
test(DfpMath.pow(pinf, factory.newDfp("-0")),
factory.newDfp("1"),
0, "pow #5");
test(DfpMath.pow(pinf, factory.newDfp("0")),
factory.newDfp("1"),
0, "pow #6");
test(DfpMath.pow(ninf, factory.newDfp("-0")),
factory.newDfp("1"),
0, "pow #7");
test(DfpMath.pow(ninf, factory.newDfp("0")),
factory.newDfp("1"),
0, "pow #8");
test(DfpMath.pow(qnan, factory.newDfp("0")),
factory.newDfp("1"),
0, "pow #8");
// exponent of one
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("1")),
factory.newDfp("0"),
0, "pow #9");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("1")),
factory.newDfp("-0"),
0, "pow #10");
test(DfpMath.pow(factory.newDfp("2"), factory.newDfp("1")),
factory.newDfp("2"),
0, "pow #11");
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("1")),
factory.newDfp("-2"),
0, "pow #12");
test(DfpMath.pow(pinf, factory.newDfp("1")),
pinf,
0, "pow #13");
test(DfpMath.pow(ninf, factory.newDfp("1")),
ninf,
0, "pow #14");
test(DfpMath.pow(qnan, factory.newDfp("1")),
qnan,
DfpField.FLAG_INVALID, "pow #14.1");
// exponent of NaN
test(DfpMath.pow(factory.newDfp("0"), qnan),
qnan,
DfpField.FLAG_INVALID, "pow #15");
test(DfpMath.pow(factory.newDfp("-0"), qnan),
qnan,
DfpField.FLAG_INVALID, "pow #16");
test(DfpMath.pow(factory.newDfp("2"), qnan),
qnan,
DfpField.FLAG_INVALID, "pow #17");
test(DfpMath.pow(factory.newDfp("-2"), qnan),
qnan,
DfpField.FLAG_INVALID, "pow #18");
test(DfpMath.pow(pinf, qnan),
qnan,
DfpField.FLAG_INVALID, "pow #19");
test(DfpMath.pow(ninf, qnan),
qnan,
DfpField.FLAG_INVALID, "pow #20");
test(DfpMath.pow(qnan, qnan),
qnan,
DfpField.FLAG_INVALID, "pow #21");
// radix of NaN
test(DfpMath.pow(qnan, factory.newDfp("1")),
qnan,
DfpField.FLAG_INVALID, "pow #22");
test(DfpMath.pow(qnan, factory.newDfp("-1")),
qnan,
DfpField.FLAG_INVALID, "pow #23");
test(DfpMath.pow(qnan, pinf),
qnan,
DfpField.FLAG_INVALID, "pow #24");
test(DfpMath.pow(qnan, ninf),
qnan,
DfpField.FLAG_INVALID, "pow #25");
test(DfpMath.pow(qnan, qnan),
qnan,
DfpField.FLAG_INVALID, "pow #26");
// (x > 1) ^ pinf = pinf, (x < -1) ^ pinf = pinf
test(DfpMath.pow(factory.newDfp("2"), pinf),
pinf,
0, "pow #27");
test(DfpMath.pow(factory.newDfp("-2"), pinf),
pinf,
0, "pow #28");
test(DfpMath.pow(pinf, pinf),
pinf,
0, "pow #29");
test(DfpMath.pow(ninf, pinf),
pinf,
0, "pow #30");
// (x > 1) ^ ninf = +0, (x < -1) ^ ninf = +0
test(DfpMath.pow(factory.newDfp("2"), ninf),
factory.getZero(),
0, "pow #31");
test(DfpMath.pow(factory.newDfp("-2"), ninf),
factory.getZero(),
0, "pow #32");
test(DfpMath.pow(pinf, ninf),
factory.getZero(),
0, "pow #33");
test(DfpMath.pow(ninf, ninf),
factory.getZero(),
0, "pow #34");
// (-1 < x < 1) ^ pinf = 0
test(DfpMath.pow(factory.newDfp("0.5"), pinf),
factory.getZero(),
0, "pow #35");
test(DfpMath.pow(factory.newDfp("-0.5"), pinf),
factory.getZero(),
0, "pow #36");
// (-1 < x < 1) ^ ninf = pinf
test(DfpMath.pow(factory.newDfp("0.5"), ninf),
pinf,
0, "pow #37");
test(DfpMath.pow(factory.newDfp("-0.5"), ninf),
pinf,
0, "pow #38");
// +/- 1 ^ +/-inf = NaN
test(DfpMath.pow(factory.getOne(), pinf),
qnan,
DfpField.FLAG_INVALID, "pow #39");
test(DfpMath.pow(factory.getOne(), ninf),
qnan,
DfpField.FLAG_INVALID, "pow #40");
test(DfpMath.pow(factory.newDfp("-1"), pinf),
qnan,
DfpField.FLAG_INVALID, "pow #41");
test(DfpMath.pow(factory.getOne().negate(), ninf),
qnan,
DfpField.FLAG_INVALID, "pow #42");
// +0 ^ +anything except 0, NAN = +0
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("1")),
factory.newDfp("0"),
0, "pow #43");
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("1e30")),
factory.newDfp("0"),
0, "pow #44");
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("1e-30")),
factory.newDfp("0"),
0, "pow #45");
test(DfpMath.pow(factory.newDfp("0"), pinf),
factory.newDfp("0"),
0, "pow #46");
// -0 ^ +anything except 0, NAN, odd integer = +0
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("2")),
factory.newDfp("0"),
0, "pow #47");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("1e30")),
factory.newDfp("0"),
0, "pow #48");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("1e-30")),
factory.newDfp("0"),
DfpField.FLAG_INEXACT, "pow #49");
test(DfpMath.pow(factory.newDfp("-0"), pinf),
factory.newDfp("0"),
0, "pow #50");
// +0 ^ -anything except 0, NAN = +INF
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("-1")),
pinf,
0, "pow #51");
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("-1e30")),
pinf,
0, "pow #52");
test(DfpMath.pow(factory.newDfp("0"), factory.newDfp("-1e-30")),
pinf,
0, "pow #53");
test(DfpMath.pow(factory.newDfp("0"), ninf),
pinf,
0, "pow #54");
// -0 ^ -anything except 0, NAN, odd integer = +INF
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("-2")),
pinf,
0, "pow #55");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("-1e30")),
pinf,
0, "pow #56");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("-1e-30")),
pinf,
DfpField.FLAG_INEXACT, "pow #57");
test(DfpMath.pow(factory.newDfp("-0"), ninf),
pinf,
0, "pow #58");
// -0 ^ -odd integer = -INF
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("-1")),
ninf,
DfpField.FLAG_INEXACT, "pow #59");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("-12345")),
ninf,
DfpField.FLAG_INEXACT, "pow #60");
// -0 ^ +odd integer = -0
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("3")),
factory.newDfp("-0"),
DfpField.FLAG_INEXACT, "pow #61");
test(DfpMath.pow(factory.newDfp("-0"), factory.newDfp("12345")),
factory.newDfp("-0"),
DfpField.FLAG_INEXACT, "pow #62");
// pinf ^ +anything = pinf
test(DfpMath.pow(pinf, factory.newDfp("3")),
pinf,
0, "pow #63");
test(DfpMath.pow(pinf, factory.newDfp("1e30")),
pinf,
0, "pow #64");
test(DfpMath.pow(pinf, factory.newDfp("1e-30")),
pinf,
0, "pow #65");
test(DfpMath.pow(pinf, pinf),
pinf,
0, "pow #66");
// pinf ^ -anything = +0
test(DfpMath.pow(pinf, factory.newDfp("-3")),
factory.getZero(),
0, "pow #67");
test(DfpMath.pow(pinf, factory.newDfp("-1e30")),
factory.getZero(),
0, "pow #68");
test(DfpMath.pow(pinf, factory.newDfp("-1e-30")),
factory.getZero(),
0, "pow #69");
test(DfpMath.pow(pinf, ninf),
factory.getZero(),
0, "pow #70");
// ninf ^ anything = -0 ^ -anything
// ninf ^ -anything except 0, NAN, odd integer = +0
test(DfpMath.pow(ninf, factory.newDfp("-2")),
factory.newDfp("0"),
0, "pow #71");
test(DfpMath.pow(ninf, factory.newDfp("-1e30")),
factory.newDfp("0"),
0, "pow #72");
test(DfpMath.pow(ninf, factory.newDfp("-1e-30")),
factory.newDfp("0"),
DfpField.FLAG_INEXACT, "pow #73");
test(DfpMath.pow(ninf, ninf),
factory.newDfp("0"),
0, "pow #74");
// ninf ^ +anything except 0, NAN, odd integer = +INF
test(DfpMath.pow(ninf, factory.newDfp("2")),
pinf,
0, "pow #75");
test(DfpMath.pow(ninf, factory.newDfp("1e30")),
pinf,
0, "pow #76");
test(DfpMath.pow(ninf, factory.newDfp("1e-30")),
pinf,
DfpField.FLAG_INEXACT, "pow #77");
test(DfpMath.pow(ninf, pinf),
pinf,
0, "pow #78");
// ninf ^ +odd integer = -INF
test(DfpMath.pow(ninf, factory.newDfp("3")),
ninf,
DfpField.FLAG_INEXACT, "pow #79");
test(DfpMath.pow(ninf, factory.newDfp("12345")),
ninf,
DfpField.FLAG_INEXACT, "pow #80");
// ninf ^ -odd integer = -0
test(DfpMath.pow(ninf, factory.newDfp("-3")),
factory.newDfp("-0"),
DfpField.FLAG_INEXACT, "pow #81");
test(DfpMath.pow(ninf, factory.newDfp("-12345")),
factory.newDfp("-0"),
DfpField.FLAG_INEXACT, "pow #82");
// -anything ^ integer
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("3")),
factory.newDfp("-8"),
DfpField.FLAG_INEXACT, "pow #83");
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("16")),
factory.newDfp("65536"),
0, "pow #84");
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("-3")),
factory.newDfp("-0.125"),
DfpField.FLAG_INEXACT, "pow #85");
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("-4")),
factory.newDfp("0.0625"),
0, "pow #86");
// -anything ^ noninteger = NaN
test(DfpMath.pow(factory.newDfp("-2"), factory.newDfp("-4.1")),
qnan,
DfpField.FLAG_INVALID|DfpField.FLAG_INEXACT, "pow #87");
// Some fractional cases.
test(DfpMath.pow(factory.newDfp("2"),factory.newDfp("1.5")),
factory.newDfp("2.8284271247461901"),
DfpField.FLAG_INEXACT, "pow #88");
}
public void testSin()
{
test(DfpMath.sin(pinf),
nan,
DfpField.FLAG_INVALID|DfpField.FLAG_INEXACT, "sin #1");
test(DfpMath.sin(nan),
nan,
DfpField.FLAG_INVALID|DfpField.FLAG_INEXACT, "sin #2");
test(DfpMath.sin(factory.getZero()),
factory.getZero(),
DfpField.FLAG_INEXACT, "sin #3");
test(DfpMath.sin(factory.getPi()),
factory.getZero(),
DfpField.FLAG_INEXACT, "sin #4");
test(DfpMath.sin(factory.getPi().negate()),
factory.newDfp("-0"),
DfpField.FLAG_INEXACT, "sin #5");
test(DfpMath.sin(factory.getPi().multiply(2)),
factory.getZero(),
DfpField.FLAG_INEXACT, "sin #6");
test(DfpMath.sin(factory.getPi().divide(2)),
factory.getOne(),
DfpField.FLAG_INEXACT, "sin #7");
test(DfpMath.sin(factory.getPi().divide(2).negate()),
factory.getOne().negate(),
DfpField.FLAG_INEXACT, "sin #8");
test(DfpMath.sin(DfpMath.atan(factory.getOne())), // pi/4
factory.newDfp("0.5").sqrt(),
DfpField.FLAG_INEXACT, "sin #9");
test(DfpMath.sin(DfpMath.atan(factory.getOne())).negate(), // -pi/4
factory.newDfp("0.5").sqrt().negate(),
DfpField.FLAG_INEXACT, "sin #10");
test(DfpMath.sin(DfpMath.atan(factory.getOne())).negate(), // -pi/4
factory.newDfp("0.5").sqrt().negate(),
DfpField.FLAG_INEXACT, "sin #11");
test(DfpMath.sin(factory.newDfp("0.1")),
factory.newDfp("0.0998334166468281523"),
DfpField.FLAG_INEXACT, "sin #12");
test(DfpMath.sin(factory.newDfp("0.2")),
factory.newDfp("0.19866933079506121546"),
DfpField.FLAG_INEXACT, "sin #13");
test(DfpMath.sin(factory.newDfp("0.3")),
factory.newDfp("0.2955202066613395751"),
DfpField.FLAG_INEXACT, "sin #14");
test(DfpMath.sin(factory.newDfp("0.4")),
factory.newDfp("0.38941834230865049166"),
DfpField.FLAG_INEXACT, "sin #15");
test(DfpMath.sin(factory.newDfp("0.5")),
factory.newDfp("0.47942553860420300026"), // off by one ULP
DfpField.FLAG_INEXACT, "sin #16");
test(DfpMath.sin(factory.newDfp("0.6")),
factory.newDfp("0.56464247339503535721"), // off by one ULP
DfpField.FLAG_INEXACT, "sin #17");
test(DfpMath.sin(factory.newDfp("0.7")),
factory.newDfp("0.64421768723769105367"),
DfpField.FLAG_INEXACT, "sin #18");
test(DfpMath.sin(factory.newDfp("0.8")),
factory.newDfp("0.71735609089952276163"),
DfpField.FLAG_INEXACT, "sin #19");
test(DfpMath.sin(factory.newDfp("0.9")), // off by one ULP
factory.newDfp("0.78332690962748338847"),
DfpField.FLAG_INEXACT, "sin #20");
test(DfpMath.sin(factory.newDfp("1.0")),
factory.newDfp("0.84147098480789650666"),
DfpField.FLAG_INEXACT, "sin #21");
test(DfpMath.sin(factory.newDfp("1.1")),
factory.newDfp("0.89120736006143533995"),
DfpField.FLAG_INEXACT, "sin #22");
test(DfpMath.sin(factory.newDfp("1.2")),
factory.newDfp("0.93203908596722634968"),
DfpField.FLAG_INEXACT, "sin #23");
test(DfpMath.sin(factory.newDfp("1.3")),
factory.newDfp("0.9635581854171929647"),
DfpField.FLAG_INEXACT, "sin #24");
test(DfpMath.sin(factory.newDfp("1.4")),
factory.newDfp("0.98544972998846018066"),
DfpField.FLAG_INEXACT, "sin #25");
test(DfpMath.sin(factory.newDfp("1.5")),
factory.newDfp("0.99749498660405443096"),
DfpField.FLAG_INEXACT, "sin #26");
test(DfpMath.sin(factory.newDfp("1.6")),
factory.newDfp("0.99957360304150516323"),
DfpField.FLAG_INEXACT, "sin #27");
}
}

1507
src/test/java/org/apache/commons/math/dfp/DfpTest.java Normal file
View File

File diff suppressed because it is too large Load Diff