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

Merge branch 'pow-reimplementation' into h10-builds 

Conflicts:
	src/main/java/org/apache/commons/math4/util/FastMath.java
	src/test/java/org/apache/commons/math4/util/FastMathTest.java
This commit is contained in:
Luc Maisonobe 2015-05-15 18:26:23 +02:00
parent 809ce8a304 175d7a1242
commit a6804eaa74
13 changed files with 452 additions and 259 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
doap_math.rdf
View File

@ -41,66 +41,92 @@
<created>2015-04-17</created>
<revision>3.5</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2015-01-11</created>
<revision>3.4.1</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2014-12-26</created>
<revision>3.4</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2014-05-14</created>
<revision>3.3</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2013-04-06</created>
<revision>3.2</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2013-01-13</created>
<revision>3.1.1</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2012-12-23</created>
<revision>3.1</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2012-03-08</created>
<revision>3.0</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2011-03-02</created>
<revision>2.2</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2010-04-02</created>
<revision>2.1</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2009-08-03</created>
<revision>2.0</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2008-02-24</created>
<revision>1.2</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2005-12-19</created>
<revision>1.1</revision>
</Version>
</release>
<release>
<Version>
<name>commons-math</name>
<created>2004-12-06</created>

3
pom.xml
View File

@ -296,6 +296,9 @@
<contributor>
<name>S&#233;bastien Riou</name>
</contributor>
<contributor>
<name>Benedikt Ritter</name>
</contributor>
<contributor>
<name>Bill Rossi</name>
</contributor>

13
src/changes/changes.xml
View File

@ -54,6 +54,19 @@ If the output is not quite correct, check for invisible trailing spaces!
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
<action dev="luc" type="fix" issue="MATH-1222" due-to="Benedikt Ritter">
Use Double.isNaN rather than x != x in FastMath.
</action>
<action dev="tn" type="fix"> <!-- backported to 3.6 -->
Fix potential branching errors in "FastMath#pow(double, double)" when
passing special values, i.e. infinity, due to erroneous JIT optimization.
</action>
<action dev="luc" type="fix" issue="MATH-1118" > <!-- backported to 3.6 -->
Fixed equals/hashcode contract failure for Dfp.
</action>
<action dev="luc" type="fix" issue="MATH-1223"> <!-- backported to 3.6 -->
Fixed wrong splitting of huge number in extended accuracy algorithms.
</action>
<action dev="luc" type="fix" issue="MATH-1143">
Added helper methods to FunctionUtils for univariate and multivariate differentiable functions conversion.
</action>

2
src/main/java/org/apache/commons/math4/dfp/Dfp.java
View File

@ -924,7 +924,7 @@ public class Dfp implements RealFieldElement<Dfp> {
*/
@Override
public int hashCode() {
return 17 + (sign << 8) + (nans << 16) + exp + Arrays.hashCode(mant);
return 17 + (isZero() ? 0 : (sign << 8)) + (nans << 16) + exp + Arrays.hashCode(mant);
}
/** Check if instance is not equal to x.

461
src/main/java/org/apache/commons/math4/util/FastMath.java
View File

@ -315,10 +315,17 @@ public class FastMath {
/** Mask used to clear the non-sign part of a long. */
private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
/** Mask used to extract exponent from double bits. */
private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
/** Mask used to extract mantissa from double bits. */
private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
/** Mask used to add implicit high order bit for normalized double. */
private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
/** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
private static final double TWO_POWER_52 = 4503599627370496.0;
/** 2^53 - double numbers this large must be even. */
private static final double TWO_POWER_53 = 2 * TWO_POWER_52;
/** Constant: {@value}. */
private static final double F_1_3 = 1d / 3d;
@ -392,7 +399,7 @@ public class FastMath {
* @return hyperbolic cosine of x
*/
public static double cosh(double x) {
if (x != x) {
if (Double.isNaN(x)) {
return x;
}
@ -462,7 +469,7 @@ public class FastMath {
*/
public static double sinh(double x) {
boolean negate = false;
if (x != x) {
if (Double.isNaN(x)) {
return x;
}
@ -588,7 +595,7 @@ public class FastMath {
public static double tanh(double x) {
boolean negate = false;
if (x != x) {
if (Double.isNaN(x)) {
return x;
}
@ -991,7 +998,7 @@ public class FastMath {
* @return exp(x) - 1
*/
private static double expm1(double x, double hiPrecOut[]) {
if (x != x || x == 0.0) { // NaN or zero
if (Double.isNaN(x) || x == 0.0) { // NaN or zero
return x;
}
@ -1155,7 +1162,7 @@ public class FastMath {
long bits = Double.doubleToRawLongBits(x);
/* Handle special cases of negative input, and NaN */
if (((bits & 0x8000000000000000L) != 0 || x != x) && x != 0.0) {
if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
if (hiPrec != null) {
hiPrec[0] = Double.NaN;
}
@ -1458,141 +1465,141 @@ public class FastMath {
* @return double
*/
public static double pow(final double x, final double y) {
final double lns[] = new double[2];
if (y == 0.0) {
if (y == 0) {
// y = -0 or y = +0
return 1.0;
} else if (x != x) { // X is NaN
return x;
} else if (y != y) { // y is NaN
return y;
} else if (x == 0) {
long bits = Double.doubleToRawLongBits(x);
if ((bits & 0x8000000000000000L) != 0) {
// -zero
long yi = (long) y;
if (y < 0 && y == yi && (yi & 1) == 1) {
return Double.NEGATIVE_INFINITY;
}
if (y > 0 && y == yi && (yi & 1) == 1) {
return -0.0;
}
}
if (y < 0) {
return Double.POSITIVE_INFINITY;
} else {
return 0.0;
}
} else if (x == Double.POSITIVE_INFINITY) {
if (y < 0.0) {
return 0.0;
} else {
return Double.POSITIVE_INFINITY;
}
} else if (y == Double.POSITIVE_INFINITY) {
if (x * x == 1.0) {
return Double.NaN;
}
if (x * x > 1.0) {
return Double.POSITIVE_INFINITY;
} else {
return 0.0;
}
} else if (x == Double.NEGATIVE_INFINITY) {
if (y < 0) {
long yi = (long) y;
if (y == yi && (yi & 1) == 1) {
return -0.0;
}
return 0.0;
}
if (y > 0) {
long yi = (long) y;
if (y == yi && (yi & 1) == 1) {
return Double.NEGATIVE_INFINITY;
}
return Double.POSITIVE_INFINITY;
}
} else if (y == Double.NEGATIVE_INFINITY) {
if (x * x == 1.0) {
return Double.NaN;
}
if (x * x < 1.0) {
return Double.POSITIVE_INFINITY;
} else {
return 0.0;
}
} else if (x < 0) { /* Handle special case x<0 */
// y is an even integer in this case
if (y >= TWO_POWER_53 || y <= -TWO_POWER_53) {
return pow(-x, y);
}
if (y == (long) y) {
// If y is an integer
return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
} else {
return Double.NaN;
}
}
/* Split y into ya and yb such that y = ya+yb */
double ya;
double yb;
if (y < 8e298 && y > -8e298) {
double tmp1 = y * HEX_40000000;
ya = y + tmp1 - tmp1;
yb = y - ya;
} else {
double tmp1 = y * 9.31322574615478515625E-10;
double tmp2 = tmp1 * 9.31322574615478515625E-10;
ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
yb = y - ya;
final long yBits = Double.doubleToRawLongBits(y);
final int yRawExp = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
final long xBits = Double.doubleToRawLongBits(x);
final int xRawExp = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;
if (yRawExp == 2047) {
// special values for y
if (yRawMantissa == 0) {
// y = -infinity or y = +infinity
if ((xRawExp == 2047 && xRawMantissa != 0) || (xRawExp == 1023 && xRawMantissa == 0)) {
// x = NaN or x = -1.0 or x = +1.0
return Double.NaN;
} else {
// on both sides of the limit points abs(x) = 1,
// infinite power asymtotically reaches different limits
if ((y > 0) ^ (xRawExp < 1023)) {
// either y = +infinity and abs(x) > 1.0
// or y = -infinity and abs(x) < 1.0
return Double.POSITIVE_INFINITY;
} else {
// either y = +infinity and abs(x) < 1.0
// or y = -infinity and abs(x) > 1.0
return +0.0;
}
}
} else {
// NaN
return Double.NaN;
}
} else {
// y is a regular non-zero number
if (yRawExp >= 1023) {
// y may be an integral value, which should be handled specifically
final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
if (yRawExp < 1075) {
// normal number with negative shift that may have a fractional part
final long integralMask = (-1L) << (1075 - yRawExp);
if ((yFullMantissa & integralMask) == yFullMantissa) {
// all fractional bits are 0, the number is really integral
final long l = yFullMantissa >> (1075 - yRawExp);
return FastMath.pow(x, (y < 0) ? -l : l);
}
} else {
// normal number with positive shift, always an integral value
double shiftedX = (y < 0) ? 1.0 / x : x;
int shiftedExp = yRawExp;
if (shiftedExp > 1085) {
// the integral value is too large to fit in a primitive long
// we want an intermediate odd power (in order not to lose the sign)
// and such that the remaining operation will not have power 0 (so result is not forced to 1.0)
final int intermediatePower = ((shiftedExp & 0x1) == 0) ? shiftedExp - 1077 : shiftedExp - 1076;
shiftedX = FastMath.pow(shiftedX, intermediatePower);
shiftedExp -= intermediatePower;
}
// the integral value fits in a primitive long
return FastMath.pow(shiftedX, yFullMantissa << (shiftedExp - 1075));
}
}
// y is a non-integral value
if (x == 0) {
// x = -0 or x = +0
// the integer powers have already been handled above
return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
} else if (xRawExp == 2047) {
if (xRawMantissa == 0) {
// x = -infinity or x = +infinity
return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
} else {
// NaN
return Double.NaN;
}
} else if (x < 0) {
// the integer powers have already been handled above
return Double.NaN;
} else {
// this is the general case, for regular fractional numbers x and y
// Split y into ya and yb such that y = ya+yb
final double tmp = y * HEX_40000000;
final double ya = (y + tmp) - tmp;
final double yb = y - ya;
/* Compute ln(x) */
final double lns[] = new double[2];
final double lores = log(x, lns);
if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
return lores;
}
double lna = lns[0];
double lnb = lns[1];
/* resplit lns */
final double tmp1 = lna * HEX_40000000;
final double tmp2 = (lna + tmp1) - tmp1;
lnb += lna - tmp2;
lna = tmp2;
// y*ln(x) = (aa+ab)
final double aa = lna * ya;
final double ab = lna * yb + lnb * ya + lnb * yb;
lna = aa+ab;
lnb = -(lna - aa - ab);
double z = 1.0 / 120.0;
z = z * lnb + (1.0 / 24.0);
z = z * lnb + (1.0 / 6.0);
z = z * lnb + 0.5;
z = z * lnb + 1.0;
z *= lnb;
final double result = exp(lna, z, null);
//result = result + result * z;
return result;
}
}
}
/* Compute ln(x) */
final double lores = log(x, lns);
if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
return lores;
}
double lna = lns[0];
double lnb = lns[1];
/* resplit lns */
double tmp1 = lna * HEX_40000000;
double tmp2 = lna + tmp1 - tmp1;
lnb += lna - tmp2;
lna = tmp2;
// y*ln(x) = (aa+ab)
final double aa = lna * ya;
final double ab = lna * yb + lnb * ya + lnb * yb;
lna = aa+ab;
lnb = -(lna - aa - ab);
double z = 1.0 / 120.0;
z = z * lnb + (1.0 / 24.0);
z = z * lnb + (1.0 / 6.0);
z = z * lnb + 0.5;
z = z * lnb + 1.0;
z *= lnb;
final double result = exp(lna, z, null);
//result = result + result * z;
return result;
}
/**
* Raise a double to an int power.
*
@ -1602,62 +1609,152 @@ public class FastMath {
* @since 3.1
*/
public static double pow(double d, int e) {
return pow(d, (long) e);
}
/**
* Raise a double to a long power.
*
* @param d Number to raise.
* @param e Exponent.
* @return d<sup>e</sup>
* @since 4.0
*/
public static double pow(double d, long e) {
if (e == 0) {
return 1.0;
} else if (e < 0) {
e = -e;
d = 1.0 / d;
} else if (e > 0) {
return new Split(d).pow(e).full;
} else {
return new Split(d).reciprocal().pow(-e).full;
}
}
/** Class operator on double numbers split into one 26 bits number and one 27 bits number. */
private static class Split {
/** Split version of NaN. */
public static final Split NAN = new Split(Double.NaN, 0);
/** Split version of positive infinity. */
public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
/** Split version of negative infinity. */
public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
/** Full number. */
private final double full;
/** High order bits. */
private final double high;
/** Low order bits. */
private final double low;
/** Simple constructor.
* @param x number to split
*/
public Split(final double x) {
full = x;
high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
low = x - high;
}
// split d as two 26 bits numbers
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final int splitFactor = 0x8000001;
final double cd = splitFactor * d;
final double d1High = cd - (cd - d);
final double d1Low = d - d1High;
/** Simple constructor.
* @param high high order bits
* @param low low order bits
*/
public Split(final double high, final double low) {
this(high + low, high, low);
}
// prepare result
double resultHigh = 1;
double resultLow = 0;
/** Simple constructor.
* @param full full number
* @param high high order bits
* @param low low order bits
*/
public Split(final double full, final double high, final double low) {
this.full = full;
this.high = high;
this.low = low;
}
// d^(2p)
double d2p = d;
double d2pHigh = d1High;
double d2pLow = d1Low;
/** Multiply the instance by another one.
* @param b other instance to multiply by
* @return product
*/
public Split multiply(final Split b) {
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final Split mulBasic = new Split(full * b.full);
final double mulError = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
return new Split(mulBasic.high, mulBasic.low + mulError);
}
while (e != 0) {
/** Compute the reciprocal of the instance.
* @return reciprocal of the instance
*/
public Split reciprocal() {
final double approximateInv = 1.0 / full;
final Split splitInv = new Split(approximateInv);
// if 1.0/d were computed perfectly, remultiplying it by d should give 1.0
// we want to estimate the error so we can fix the low order bits of approximateInvLow
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final Split product = multiply(splitInv);
final double error = (product.high - 1) + product.low;
// better accuracy estimate of reciprocal
return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
}
/** Computes this^e.
* @param e exponent (beware, here it MUST be > 0)
* @return d^e, split in high and low bits
* @since 4.0
*/
private Split pow(final long e) {
// prepare result
Split result = new Split(1);
// d^(2p)
Split d2p = new Split(full, high, low);
for (long p = e; p != 0; p >>= 1) {
if ((p & 0x1) != 0) {
// accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
result = result.multiply(d2p);
}
// accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
d2p = d2p.multiply(d2p);
if ((e & 0x1) != 0) {
// accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final double tmpHigh = resultHigh * d2p;
final double cRH = splitFactor * resultHigh;
final double rHH = cRH - (cRH - resultHigh);
final double rHL = resultHigh - rHH;
final double tmpLow = rHL * d2pLow - (((tmpHigh - rHH * d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
resultHigh = tmpHigh;
resultLow = resultLow * d2p + tmpLow;
}
// accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
final double tmpHigh = d2pHigh * d2p;
final double cD2pH = splitFactor * d2pHigh;
final double d2pHH = cD2pH - (cD2pH - d2pHigh);
final double d2pHL = d2pHigh - d2pHH;
final double tmpLow = d2pHL * d2pLow - (((tmpHigh - d2pHH * d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
final double cTmpH = splitFactor * tmpHigh;
d2pHigh = cTmpH - (cTmpH - tmpHigh);
d2pLow = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
d2p = d2pHigh + d2pLow;
e >>= 1;
if (Double.isNaN(result.full)) {
if (Double.isNaN(full)) {
return Split.NAN;
} else {
// some intermediate numbers exceeded capacity,
// and the low order bits became NaN (because infinity - infinity = NaN)
if (FastMath.abs(full) < 1) {
return new Split(FastMath.copySign(0.0, full), 0.0);
} else if (full < 0 && (e & 0x1) == 1) {
return Split.NEGATIVE_INFINITY;
} else {
return Split.POSITIVE_INFINITY;
}
}
} else {
return result;
}
}
return resultHigh + resultLow;
}
/**
@ -2576,7 +2673,7 @@ public class FastMath {
* @return phase angle of point (x,y) between {@code -PI} and {@code PI}
*/
public static double atan2(double y, double x) {
if (x != x || y != y) {
if (Double.isNaN(x) || Double.isNaN(y)) {
return Double.NaN;
}
@ -2697,7 +2794,7 @@ public class FastMath {
* @return arc sine of x
*/
public static double asin(double x) {
if (x != x) {
if (Double.isNaN(x)) {
return Double.NaN;
}
@ -2773,7 +2870,7 @@ public class FastMath {
* @return arc cosine of x
*/
public static double acos(double x) {
if (x != x) {
if (Double.isNaN(x)) {
return Double.NaN;
}
@ -3333,7 +3430,7 @@ public class FastMath {
public static double floor(double x) {
long y;
if (x != x) { // NaN
if (Double.isNaN(x)) {
return x;
}
@ -3360,7 +3457,7 @@ public class FastMath {
public static double ceil(double x) {
double y;
if (x != x) { // NaN
if (Double.isNaN(x)) {
return x;
}

91
src/main/java/org/apache/commons/math4/util/MathArrays.java
View File

@ -47,8 +47,6 @@ import org.apache.commons.math4.random.Well19937c;
* @since 3.0
*/
public class MathArrays {
/** Factor used for splitting double numbers: n = 2^27 + 1 (i.e. {@value}). */
private static final int SPLIT_FACTOR = 0x8000001;
/**
* Private constructor.
@ -864,15 +862,13 @@ public class MathArrays {
double prodLowSum = 0;
for (int i = 0; i < len; i++) {
final double ai = a[i];
final double ca = SPLIT_FACTOR * ai;
final double aHigh = ca - (ca - ai);
final double aLow = ai - aHigh;
final double ai = a[i];
final double aHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(ai) & ((-1L) << 27));
final double aLow = ai - aHigh;
final double bi = b[i];
final double cb = SPLIT_FACTOR * bi;
final double bHigh = cb - (cb - bi);
final double bLow = bi - bHigh;
final double bi = b[i];
final double bHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(bi) & ((-1L) << 27));
final double bLow = bi - bHigh;
prodHigh[i] = ai * bi;
final double prodLow = aLow * bLow - (((prodHigh[i] -
aHigh * bHigh) -
@ -938,7 +934,6 @@ public class MathArrays {
// the code below is split in many additions/subtractions that may
// appear redundant. However, they should NOT be simplified, as they
// use IEEE754 floating point arithmetic rounding properties.
// as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
// The variable naming conventions are that xyzHigh contains the most significant
// bits of xyz and xyzLow contains its least significant bits. So theoretically
// xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
@ -946,24 +941,20 @@ public class MathArrays {
// to hold it as long as we can, combining the high and low order bits together
// only at the end, after cancellation may have occurred on high order bits
// split a1 and b1 as two 26 bits numbers
final double ca1 = SPLIT_FACTOR * a1;
final double a1High = ca1 - (ca1 - a1);
// split a1 and b1 as one 26 bits number and one 27 bits number
final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
final double a1Low = a1 - a1High;
final double cb1 = SPLIT_FACTOR * b1;
final double b1High = cb1 - (cb1 - b1);
final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
final double b1Low = b1 - b1High;
// accurate multiplication a1 * b1
final double prod1High = a1 * b1;
final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
// split a2 and b2 as two 26 bits numbers
final double ca2 = SPLIT_FACTOR * a2;
final double a2High = ca2 - (ca2 - a2);
// split a2 and b2 as one 26 bits number and one 27 bits number
final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
final double a2Low = a2 - a2High;
final double cb2 = SPLIT_FACTOR * b2;
final double b2High = cb2 - (cb2 - b2);
final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
final double b2Low = b2 - b2High;
// accurate multiplication a2 * b2
@ -1018,7 +1009,6 @@ public class MathArrays {
// the code below is split in many additions/subtractions that may
// appear redundant. However, they should NOT be simplified, as they
// do use IEEE754 floating point arithmetic rounding properties.
// as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
// The variables naming conventions are that xyzHigh contains the most significant
// bits of xyz and xyzLow contains its least significant bits. So theoretically
// xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
@ -1026,36 +1016,30 @@ public class MathArrays {
// to hold it as long as we can, combining the high and low order bits together
// only at the end, after cancellation may have occurred on high order bits
// split a1 and b1 as two 26 bits numbers
final double ca1 = SPLIT_FACTOR * a1;
final double a1High = ca1 - (ca1 - a1);
// split a1 and b1 as one 26 bits number and one 27 bits number
final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
final double a1Low = a1 - a1High;
final double cb1 = SPLIT_FACTOR * b1;
final double b1High = cb1 - (cb1 - b1);
final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
final double b1Low = b1 - b1High;
// accurate multiplication a1 * b1
final double prod1High = a1 * b1;
final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
// split a2 and b2 as two 26 bits numbers
final double ca2 = SPLIT_FACTOR * a2;
final double a2High = ca2 - (ca2 - a2);
// split a2 and b2 as one 26 bits number and one 27 bits number
final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
final double a2Low = a2 - a2High;
final double cb2 = SPLIT_FACTOR * b2;
final double b2High = cb2 - (cb2 - b2);
final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
final double b2Low = b2 - b2High;
// accurate multiplication a2 * b2
final double prod2High = a2 * b2;
final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
// split a3 and b3 as two 26 bits numbers
final double ca3 = SPLIT_FACTOR * a3;
final double a3High = ca3 - (ca3 - a3);
// split a3 and b3 as one 26 bits number and one 27 bits number
final double a3High = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27));
final double a3Low = a3 - a3High;
final double cb3 = SPLIT_FACTOR * b3;
final double b3High = cb3 - (cb3 - b3);
final double b3High = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27));
final double b3Low = b3 - b3High;
// accurate multiplication a3 * b3
@ -1120,7 +1104,6 @@ public class MathArrays {
// the code below is split in many additions/subtractions that may
// appear redundant. However, they should NOT be simplified, as they
// do use IEEE754 floating point arithmetic rounding properties.
// as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
// The variables naming conventions are that xyzHigh contains the most significant
// bits of xyz and xyzLow contains its least significant bits. So theoretically
// xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
@ -1128,48 +1111,40 @@ public class MathArrays {
// to hold it as long as we can, combining the high and low order bits together
// only at the end, after cancellation may have occurred on high order bits
// split a1 and b1 as two 26 bits numbers
final double ca1 = SPLIT_FACTOR * a1;
final double a1High = ca1 - (ca1 - a1);
// split a1 and b1 as one 26 bits number and one 27 bits number
final double a1High = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
final double a1Low = a1 - a1High;
final double cb1 = SPLIT_FACTOR * b1;
final double b1High = cb1 - (cb1 - b1);
final double b1High = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
final double b1Low = b1 - b1High;
// accurate multiplication a1 * b1
final double prod1High = a1 * b1;
final double prod1Low = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
// split a2 and b2 as two 26 bits numbers
final double ca2 = SPLIT_FACTOR * a2;
final double a2High = ca2 - (ca2 - a2);
// split a2 and b2 as one 26 bits number and one 27 bits number
final double a2High = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
final double a2Low = a2 - a2High;
final double cb2 = SPLIT_FACTOR * b2;
final double b2High = cb2 - (cb2 - b2);
final double b2High = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
final double b2Low = b2 - b2High;
// accurate multiplication a2 * b2
final double prod2High = a2 * b2;
final double prod2Low = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
// split a3 and b3 as two 26 bits numbers
final double ca3 = SPLIT_FACTOR * a3;
final double a3High = ca3 - (ca3 - a3);
// split a3 and b3 as one 26 bits number and one 27 bits number
final double a3High = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27));
final double a3Low = a3 - a3High;
final double cb3 = SPLIT_FACTOR * b3;
final double b3High = cb3 - (cb3 - b3);
final double b3High = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27));
final double b3Low = b3 - b3High;
// accurate multiplication a3 * b3
final double prod3High = a3 * b3;
final double prod3Low = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
// split a4 and b4 as two 26 bits numbers
final double ca4 = SPLIT_FACTOR * a4;
final double a4High = ca4 - (ca4 - a4);
// split a4 and b4 as one 26 bits number and one 27 bits number
final double a4High = Double.longBitsToDouble(Double.doubleToRawLongBits(a4) & ((-1L) << 27));
final double a4Low = a4 - a4High;
final double cb4 = SPLIT_FACTOR * b4;
final double b4High = cb4 - (cb4 - b4);
final double b4High = Double.longBitsToDouble(Double.doubleToRawLongBits(b4) & ((-1L) << 27));
final double b4Low = b4 - b4High;
// accurate multiplication a4 * b4

10
src/test/java/org/apache/commons/math4/dfp/DfpTest.java
View File

@ -1652,4 +1652,14 @@ public class DfpTest extends ExtendedFieldElementAbstractTest<Dfp> {
Assert.assertTrue(field.newDfp("NaN").isNaN());
}
@Test
public void testEqualsHashcodeContract() {
DfpField var1 = new DfpField(1);
Dfp var6 = var1.newDfp(-0.0d);
Dfp var5 = var1.newDfp(0L);
// Checks the contract: equals-hashcode on var5 and var6
Assert.assertTrue(var5.equals(var6) ? var5.hashCode() == var6.hashCode() : true);
}
}

2
src/test/java/org/apache/commons/math4/distribution/GammaDistributionTest.java
View File

@ -350,7 +350,7 @@ public class GammaDistributionTest extends RealDistributionAbstractTest {
@Test
public void testMath753Shape142() throws IOException {
doTestMath753(142.0, 0.5, 1.5, 40.0, 40.0, "gamma-distribution-shape-142.csv");
doTestMath753(142.0, 3.3, 1.6, 40.0, 40.0, "gamma-distribution-shape-142.csv");
}
@Test

2
src/test/java/org/apache/commons/math4/geometry/euclidean/threed/FieldVector3DTest.java
View File

@ -585,7 +585,7 @@ public class FieldVector3DTest {
DerivativeStructure sNaive = u1.getX().multiply(u2.getX()).add(u1.getY().multiply(u2.getY())).add(u1.getZ().multiply(u2.getZ()));
DerivativeStructure sAccurate = FieldVector3D.dotProduct(u1, u2);
Assert.assertEquals(0.0, sNaive.getReal(), 1.0e-30);
Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate.getReal(), 1.0e-16);
Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate.getReal(), 1.0e-15);
}
@Test

2
src/test/java/org/apache/commons/math4/geometry/euclidean/threed/Vector3DTest.java
View File

@ -343,7 +343,7 @@ public class Vector3DTest {
double sNaive = u1.getX() * u2.getX() + u1.getY() * u2.getY() + u1.getZ() * u2.getZ();
double sAccurate = u1.dotProduct(u2);
Assert.assertEquals(0.0, sNaive, 1.0e-30);
Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate, 1.0e-16);
Assert.assertEquals(-2088690039198397.0 / 1125899906842624.0, sAccurate, 1.0e-15);
}
@Test

2
src/test/java/org/apache/commons/math4/optim/nonlinear/scalar/noderiv/BOBYQAOptimizerTest.java
View File

@ -189,7 +189,7 @@ public class BOBYQAOptimizerTest {
new PointValuePair(point(DIM/2,0.0),0.0);
doTest(new DiffPow(), startPoint, boundaries,
GoalType.MINIMIZE,
1e-8, 1e-1, 12000, expected);
1e-8, 1e-1, 21000, expected);
}
@Test

64
src/test/java/org/apache/commons/math4/util/FastMathTest.java
View File

@ -338,9 +338,9 @@ public class FastMathTest {
@Test
public void testLogSpecialCases() {
Assert.assertTrue("Log of zero should be -Inf", Double.isInfinite(FastMath.log(0.0)));
Assert.assertEquals("Log of zero should be -Inf", Double.NEGATIVE_INFINITY, FastMath.log(0.0), 1.0);
Assert.assertTrue("Log of -zero should be -Inf", Double.isInfinite(FastMath.log(-0.0)));
Assert.assertEquals("Log of -zero should be -Inf", Double.NEGATIVE_INFINITY, FastMath.log(-0.0), 1.0);
Assert.assertTrue("Log of NaN should be NaN", Double.isNaN(FastMath.log(Double.NaN)));
@ -348,8 +348,9 @@ public class FastMathTest {
Assert.assertEquals("Log of Double.MIN_VALUE should be -744.4400719213812", -744.4400719213812, FastMath.log(Double.MIN_VALUE), Precision.EPSILON);
Assert.assertTrue("Log of infinity should be infinity", Double.isInfinite(FastMath.log(Double.POSITIVE_INFINITY)));
Assert.assertEquals("Log of infinity should be infinity", Double.POSITIVE_INFINITY, FastMath.log(Double.POSITIVE_INFINITY), 1.0);
}
@Test
public void testExpSpecialCases() {
@ -360,7 +361,7 @@ public class FastMathTest {
Assert.assertTrue("exp of NaN should be NaN", Double.isNaN(FastMath.exp(Double.NaN)));
Assert.assertTrue("exp of infinity should be infinity", Double.isInfinite(FastMath.exp(Double.POSITIVE_INFINITY)));
Assert.assertEquals("exp of infinity should be infinity", Double.POSITIVE_INFINITY, FastMath.exp(Double.POSITIVE_INFINITY), 1.0);
Assert.assertEquals("exp of -infinity should be 0.0", 0.0, FastMath.exp(Double.NEGATIVE_INFINITY), Precision.EPSILON);
@ -383,9 +384,11 @@ public class FastMathTest {
Assert.assertTrue("pow(NaN, PI) should be NaN", Double.isNaN(FastMath.pow(Double.NaN, Math.PI)));
Assert.assertTrue("pow(2.0, Infinity) should be Infinity", Double.isInfinite(FastMath.pow(2.0, Double.POSITIVE_INFINITY)));
Assert.assertEquals("pow(2.0, Infinity) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(2.0, Double.POSITIVE_INFINITY), 1.0);
Assert.assertTrue("pow(0.5, -Infinity) should be Infinity", Double.isInfinite(FastMath.pow(0.5, Double.NEGATIVE_INFINITY)));
Assert.assertEquals("pow(0.5, -Infinity) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(0.5, Double.NEGATIVE_INFINITY), 1.0);
Assert.assertEquals("pow(0.5, Infinity) should be 0.0", 0.0, FastMath.pow(0.5, Double.POSITIVE_INFINITY), Precision.EPSILON);
Assert.assertEquals("pow(2.0, -Infinity) should be 0.0", 0.0, FastMath.pow(2.0, Double.NEGATIVE_INFINITY), Precision.EPSILON);
@ -393,22 +396,50 @@ public class FastMathTest {
Assert.assertEquals("pow(Infinity, -0.5) should be 0.0", 0.0, FastMath.pow(Double.POSITIVE_INFINITY, -0.5), Precision.EPSILON);
Assert.assertTrue("pow(0.0, -0.5) should be Inf", Double.isInfinite(FastMath.pow(0.0, -0.5)));
Assert.assertEquals("pow(0.0, -0.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(0.0, -0.5), 1.0);
Assert.assertTrue("pow(Inf, 0.5) should be Inf", Double.isInfinite(FastMath.pow(Double.POSITIVE_INFINITY, 0.5)));
Assert.assertEquals("pow(Inf, 0.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(Double.POSITIVE_INFINITY, 0.5), 1.0);
Assert.assertTrue("pow(-0.0, -3.0) should be -Inf", Double.isInfinite(FastMath.pow(-0.0, -3.0)));
Assert.assertEquals("pow(-0.0, -3.0) should be -Inf", Double.NEGATIVE_INFINITY, FastMath.pow(-0.0, -3.0), 1.0);
Assert.assertTrue("pow(-Inf, -3.0) should be -Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 3.0)));
Assert.assertEquals("pow(-0.0, Infinity) should be 0.0", 0.0, FastMath.pow(-0.0, Double.POSITIVE_INFINITY), Precision.EPSILON);
Assert.assertTrue("pow(-0.0, -3.5) should be Inf", Double.isInfinite(FastMath.pow(-0.0, -3.5)));
Assert.assertTrue("pow(-0.0, NaN) should be NaN", Double.isNaN(FastMath.pow(-0.0, Double.NaN)));
Assert.assertTrue("pow(Inf, 3.5) should be Inf", Double.isInfinite(FastMath.pow(Double.POSITIVE_INFINITY, 3.5)));
Assert.assertEquals("pow(-0.0, -tiny) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(-0.0, -Double.MIN_VALUE), 1.0);
Assert.assertEquals("pow(-0.0, -huge) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(-0.0, -Double.MAX_VALUE), 1.0);
Assert.assertEquals("pow(-Inf, -3.0) should be -Inf", Double.NEGATIVE_INFINITY, FastMath.pow(Double.NEGATIVE_INFINITY, 3.0), 1.0);
Assert.assertEquals("pow(-0.0, -3.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(-0.0, -3.5), 1.0);
Assert.assertEquals("pow(Inf, 3.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(Double.POSITIVE_INFINITY, 3.5), 1.0);
Assert.assertEquals("pow(-2.0, 3.0) should be -8.0", -8.0, FastMath.pow(-2.0, 3.0), Precision.EPSILON);
Assert.assertTrue("pow(-2.0, 3.5) should be NaN", Double.isNaN(FastMath.pow(-2.0, 3.5)));
Assert.assertTrue("pow(NaN, -Infinity) should be NaN", Double.isNaN(FastMath.pow(Double.NaN, Double.NEGATIVE_INFINITY)));
Assert.assertEquals("pow(NaN, 0.0) should be 1.0", 1.0, FastMath.pow(Double.NaN, 0.0), Precision.EPSILON);
Assert.assertEquals("pow(-Infinity, -Infinity) should be 0.0", 0.0, FastMath.pow(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY), Precision.EPSILON);
Assert.assertEquals("pow(-huge, -huge) should be 0.0", 0.0, FastMath.pow(-Double.MAX_VALUE, -Double.MAX_VALUE), Precision.EPSILON);
Assert.assertTrue("pow(-huge, huge) should be +Inf", Double.isInfinite(FastMath.pow(-Double.MAX_VALUE, Double.MAX_VALUE)));
Assert.assertTrue("pow(NaN, -Infinity) should be NaN", Double.isNaN(FastMath.pow(Double.NaN, Double.NEGATIVE_INFINITY)));
Assert.assertEquals("pow(NaN, 0.0) should be 1.0", 1.0, FastMath.pow(Double.NaN, 0.0), Precision.EPSILON);
Assert.assertEquals("pow(-Infinity, -Infinity) should be 0.0", 0.0, FastMath.pow(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY), Precision.EPSILON);
Assert.assertEquals("pow(-huge, -huge) should be 0.0", 0.0, FastMath.pow(-Double.MAX_VALUE, -Double.MAX_VALUE), Precision.EPSILON);
Assert.assertEquals("pow(-huge, huge) should be +Inf", Double.POSITIVE_INFINITY, FastMath.pow(-Double.MAX_VALUE, Double.MAX_VALUE), 1.0);
// Added tests for a 100% coverage
Assert.assertTrue("pow(+Inf, NaN) should be NaN", Double.isNaN(FastMath.pow(Double.POSITIVE_INFINITY, Double.NaN)));
@ -421,9 +452,9 @@ public class FastMathTest {
Assert.assertEquals("pow(-Inf, -2.0) should be 0.0", 0.0, FastMath.pow(Double.NEGATIVE_INFINITY, -2.0), Precision.EPSILON);
Assert.assertTrue("pow(-Inf, 1.0) should be -Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 1.0)));
Assert.assertEquals("pow(-Inf, 1.0) should be -Inf", Double.NEGATIVE_INFINITY, FastMath.pow(Double.NEGATIVE_INFINITY, 1.0), 1.0);
Assert.assertTrue("pow(-Inf, 2.0) should be +Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 2.0)));
Assert.assertEquals("pow(-Inf, 2.0) should be +Inf", Double.POSITIVE_INFINITY, FastMath.pow(Double.NEGATIVE_INFINITY, 2.0), 1.0);
Assert.assertTrue("pow(1.0, -Inf) should be NaN", Double.isNaN(FastMath.pow(1.0, Double.NEGATIVE_INFINITY)));
@ -1203,6 +1234,11 @@ public class FastMathTest {
}
}
@Test
public void testIntPowHuge() {
Assert.assertTrue(Double.isInfinite(FastMath.pow(FastMath.scalb(1.0, 500), 4)));
}
@Test
public void testIncrementExactInt() {
int[] specialValues = new int[] {

33
src/test/java/org/apache/commons/math4/util/MathArraysTest.java
View File

@ -711,6 +711,39 @@ public class MathArraysTest {
}
}
@Test
public void testLinearCombinationHuge() {
int scale = 971;
final double[] a = new double[] {
-1321008684645961.0 / 268435456.0,
-5774608829631843.0 / 268435456.0,
-7645843051051357.0 / 8589934592.0
};
final double[] b = new double[] {
-5712344449280879.0 / 2097152.0,
-4550117129121957.0 / 2097152.0,
8846951984510141.0 / 131072.0
};
double[] scaledA = new double[a.length];
double[] scaledB = new double[b.length];
for (int i = 0; i < scaledA.length; ++i) {
scaledA[i] = FastMath.scalb(a[i], -scale);
scaledB[i] = FastMath.scalb(b[i], +scale);
}
final double abSumInline = MathArrays.linearCombination(scaledA[0], scaledB[0],
scaledA[1], scaledB[1],
scaledA[2], scaledB[2]);
final double abSumArray = MathArrays.linearCombination(scaledA, scaledB);
Assert.assertEquals(abSumInline, abSumArray, 0);
Assert.assertEquals(-1.8551294182586248737720779899, abSumInline, 1.0e-15);
final double naive = scaledA[0] * scaledB[0] + scaledA[1] * scaledB[1] + scaledA[2] * scaledB[2];
Assert.assertTrue(FastMath.abs(naive - abSumInline) > 1.5);
}
@Test
public void testLinearCombinationInfinite() {
final double[][] a = new double[][] {