You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2010/10/03 22:52:54 UTC
svn commit: r1004045 - in /commons/proper/math/branches/MATH_2_X/src:
main/java/org/apache/commons/math/util/FastMath.java
test/java/org/apache/commons/math/util/FastMathTest.java
Author: luc
Date: Sun Oct 3 20:52:53 2010
New Revision: 1004045
URL: http://svn.apache.org/viewvc?rev=1004045&view=rev
Log:
applied Bill's patch adding sinh, cosh and tanh functions to FastMath
JIRA: MATH-375
Modified:
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java?rev=1004045&r1=1004044&r2=1004045&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java Sun Oct 3 20:52:53 2010
@@ -220,68 +220,286 @@ public class FastMath {
}
/** Compute the hyperbolic cosine of a number.
- * @param a number on which evaluation is done
- * @return hyperbolic cosine of a
+ * @param x number on which evaluation is done
+ * @return hyperbolic cosine of x
*/
- public static double cosh(final double a) {
- return 0.5 * (FastMath.exp(a) + FastMath.exp(-a));
+ public static double cosh(double x) {
+ if (x != x) {
+ return x;
+ }
+
+ if (x > 20.0) {
+ return exp(x)/2.0;
+ }
+
+ if (x < -20) {
+ return exp(-x)/2.0;
+ }
+
+ double hiPrec[] = new double[2];
+ if (x < 0.0) {
+ x = -x;
+ }
+ exp(x, 0.0, hiPrec);
+
+ double ya = hiPrec[0] + hiPrec[1];
+ double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+ double temp = ya * 1073741824.0;
+ double yaa = ya + temp - temp;
+ double yab = ya - yaa;
+
+ // recip = 1/y
+ double recip = 1.0/ya;
+ temp = recip * 1073741824.0;
+ double recipa = recip + temp - temp;
+ double recipb = recip - recipa;
+
+ // Correct for rounding in division
+ recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
+ // Account for yb
+ recipb += -yb * recip * recip;
+
+ // y = y + 1/y
+ temp = ya + recipa;
+ yb += -(temp - ya - recipa);
+ ya = temp;
+ temp = ya + recipb;
+ yb += -(temp - ya - recipb);
+ ya = temp;
+
+ double result = ya + yb;
+ result *= 0.5;
+ return result;
}
/** Compute the hyperbolic sine of a number.
- * @param a number on which evaluation is done
- * @return hyperbolic sine of a
+ * @param x number on which evaluation is done
+ * @return hyperbolic sine of x
*/
- public static double sinh(double a) {
+ public static double sinh(double x) {
+ boolean negate = false;
+ if (x != x) {
+ return x;
+ }
- boolean negative = false;
- if (a < 0) {
- negative = true;
- a = -a;
- }
+ if (x > 20.0) {
+ return exp(x)/2.0;
+ }
- double absSinh;
- if (a > 0.3) {
- absSinh = 0.5 * (FastMath.exp(a) - FastMath.exp(-a));
- } else {
- final double a2 = a * a;
- if (a > 0.05) {
- absSinh = a * (1 + a2 * (1 + a2 * (1 + a2 * (1 + a2 * (1 + a2 / 110) / 72) / 42) / 20) / 6);
- } else {
- absSinh = a * (1 + a2 * (1 + a2 * (1 + a2 / 42) / 20) / 6);
- }
- }
+ if (x < -20) {
+ return -exp(-x)/2.0;
+ }
- return negative ? -absSinh : absSinh;
+ if (x == 0) {
+ return x;
+ }
+
+ if (x < 0.0) {
+ x = -x;
+ negate = true;
+ }
+ double result;
+
+ if (x > 0.25) {
+ double hiPrec[] = new double[2];
+ exp(x, 0.0, hiPrec);
+
+ double ya = hiPrec[0] + hiPrec[1];
+ double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+ double temp = ya * 1073741824.0;
+ double yaa = ya + temp - temp;
+ double yab = ya - yaa;
+
+ // recip = 1/y
+ double recip = 1.0/ya;
+ temp = recip * 1073741824.0;
+ double recipa = recip + temp - temp;
+ double recipb = recip - recipa;
+
+ // Correct for rounding in division
+ recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
+ // Account for yb
+ recipb += -yb * recip * recip;
+
+ recipa = -recipa;
+ recipb = -recipb;
+
+ // y = y + 1/y
+ temp = ya + recipa;
+ yb += -(temp - ya - recipa);
+ ya = temp;
+ temp = ya + recipb;
+ yb += -(temp - ya - recipb);
+ ya = temp;
+
+ result = ya + yb;
+ result *= 0.5;
+ }
+ else {
+ double hiPrec[] = new double[2];
+ expm1(x, hiPrec);
+
+ double ya = hiPrec[0] + hiPrec[1];
+ double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+ /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
+ double denom = 1.0+ya;
+ double denomr = 1.0/denom;
+ double denomb = -(denom - 1.0 - ya) + yb;
+ double ratio = ya*denomr;
+ double ra, rb;
+ double temp = ratio * 1073741824.0;
+ ra = ratio + temp - temp;
+ rb = ratio - ra;
+
+ temp = denom * 1073741824.0;
+ double za = denom + temp - temp;
+ double zb = denom - za;
+
+ rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
+
+ // Adjust for yb
+ rb += yb*denomr; // numerator
+ rb += -ya * denomb * denomr * denomr; // denominator
+
+ // y = y - 1/y
+ temp = ya + ra;
+ yb += -(temp - ya - ra);
+ ya = temp;
+ temp = ya + rb;
+ yb += -(temp - ya - rb);
+ ya = temp;
+
+ result = ya + yb;
+ result *= 0.5;
+ }
+
+ if (negate) {
+ result = -result;
+ }
+
+ return result;
}
/** Compute the hyperbolic tangent of a number.
- * @param a number on which evaluation is done
- * @return hyperbolic tangent of a
+ * @param x number on which evaluation is done
+ * @return hyperbolic tangent of x
*/
- public static double tanh(double a) {
+ public static double tanh(double x) {
+ boolean negate = false;
- boolean negative = false;
- if (a < 0) {
- negative = true;
- a = -a;
- }
+ if (x != x) {
+ return x;
+ }
- double absTanh;
- if (a > 0.074) {
- final double twoA = 2 * a;
- absTanh = FastMath.expm1(twoA) / (FastMath.exp(twoA) + 1);
- } else {
- final double a2 = a * a;
- if (a > 0.016) {
- absTanh = a * (1 - a2 * (1 - a2 * (2 - a2 * (17 - a2 * (62 - a2 * 1382 / 55 ) / 9) / 21) / 5) / 3);
- } else {
- absTanh = a * (1 - a2 * (1 - a2 * (2 - a2 * 17 / 21) / 5) / 3);
- }
- }
+ if (x > 20.0) {
+ return 1.0;
+ }
- return negative ? -absTanh : absTanh;
+ if (x < -20) {
+ return -1.0;
+ }
+ if (x == 0) {
+ return x;
+ }
+
+ if (x < 0.0) {
+ x = -x;
+ negate = true;
+ }
+
+ double result;
+ if (x >= 0.5) {
+ double hiPrec[] = new double[2];
+ // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
+ exp(x*2.0, 0.0, hiPrec);
+
+ double ya = hiPrec[0] + hiPrec[1];
+ double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+ /* Numerator */
+ double na = -1.0 + ya;
+ double nb = -(na + 1.0 - ya);
+ double temp = na + yb;
+ nb += -(temp - na - yb);
+ na = temp;
+
+ /* Denominator */
+ double da = 1.0 + ya;
+ double db = -(da - 1.0 - ya);
+ temp = da + yb;
+ db += -(temp - da - yb);
+ da = temp;
+
+ temp = da * 1073741824.0;
+ double daa = da + temp - temp;
+ double dab = da - daa;
+
+ // ratio = na/da
+ double ratio = na/da;
+ temp = ratio * 1073741824.0;
+ double ratioa = ratio + temp - temp;
+ double ratiob = ratio - ratioa;
+
+ // Correct for rounding in division
+ ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
+
+ // Account for nb
+ ratiob += nb / da;
+ // Account for db
+ ratiob += -db * na / da / da;
+
+ result = ratioa + ratiob;
+ }
+ else {
+ double hiPrec[] = new double[2];
+ // tanh(x) = expm1(2x) / (expm1(2x) + 2)
+ expm1(x*2.0, hiPrec);
+
+ double ya = hiPrec[0] + hiPrec[1];
+ double yb = -(ya - hiPrec[0] - hiPrec[1]);
+
+ /* Numerator */
+ double na = ya;
+ double nb = yb;
+
+ /* Denominator */
+ double da = 2.0 + ya;
+ double db = -(da - 2.0 - ya);
+ double temp = da + yb;
+ db += -(temp - da - yb);
+ da = temp;
+
+ temp = da * 1073741824.0;
+ double daa = da + temp - temp;
+ double dab = da - daa;
+
+ // ratio = na/da
+ double ratio = na/da;
+ temp = ratio * 1073741824.0;
+ double ratioa = ratio + temp - temp;
+ double ratiob = ratio - ratioa;
+
+ // Correct for rounding in division
+ ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
+
+ // Account for nb
+ ratiob += nb / da;
+ // Account for db
+ ratiob += -db * na / da / da;
+
+ result = ratioa + ratiob;
+ }
+
+ if (negate) {
+ result = -result;
+ }
+
+ return result;
}
/** Compute the inverse hyperbolic cosine of a number.
@@ -533,6 +751,15 @@ public class FastMath {
* @return exp(x) - 1
*/
public static double expm1(double x) {
+ return expm1(x, null);
+ }
+
+ /** Internal helper method for expm1
+ * @param x number to compute shifted exponential
+ * @param hiPrecOut[] receive high precision result for -1.0 < x < 1.0
+ * @return exp(x) - 1
+ */
+ private static double expm1(double x, double hiPrecOut[]) {
if (x != x || x == 0.0) { // NaN or zero
return x;
}
@@ -666,6 +893,11 @@ public class FastMath {
yb = -rb;
}
+ if (hiPrecOut != null) {
+ hiPrecOut[0] = ya;
+ hiPrecOut[1] = yb;
+ }
+
return ya + yb;
}
Modified: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java?rev=1004045&r1=1004044&r2=1004045&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/util/FastMathTest.java Sun Oct 3 20:52:53 2010
@@ -788,6 +788,84 @@ public class FastMathTest {
Assert.assertTrue("acos() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
}
+ private Dfp cosh(Dfp x) {
+ return DfpMath.exp(x).add(DfpMath.exp(x.negate())).divide(2);
+ }
+
+ private Dfp sinh(Dfp x) {
+ return DfpMath.exp(x).subtract(DfpMath.exp(x.negate())).divide(2);
+ }
+
+ private Dfp tanh(Dfp x) {
+ return sinh(x).divide(cosh(x));
+ }
+
+ @Test
+ public void testSinhAccuracy() {
+ double maxerrulp = 0.0;
+
+ for (int i=0; i<10000; i++) {
+ double x = ((generator.nextDouble() * 16.0) - 8.0) * generator.nextDouble();
+
+ double tst = FastMath.sinh(x);
+ double ref = sinh(field.newDfp(x)).toDouble();
+ double err = (tst - ref) / ref;
+
+ if (err != 0) {
+ double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref) ^ 1)));
+ double errulp = field.newDfp(tst).subtract(sinh(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble();
+ //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp);
+ maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+ }
+ }
+
+ Assert.assertTrue("sinh() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
+ }
+
+ @Test
+ public void testCoshAccuracy() {
+ double maxerrulp = 0.0;
+
+ for (int i=0; i<10000; i++) {
+ double x = ((generator.nextDouble() * 16.0) - 8.0) * generator.nextDouble();
+
+ double tst = FastMath.cosh(x);
+ double ref = cosh(field.newDfp(x)).toDouble();
+ double err = (tst - ref) / ref;
+
+ if (err != 0) {
+ double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref) ^ 1)));
+ double errulp = field.newDfp(tst).subtract(cosh(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble();
+ //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp);
+ maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+ }
+ }
+
+ Assert.assertTrue("cosh() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
+ }
+
+ @Test
+ public void testTanhAccuracy() {
+ double maxerrulp = 0.0;
+
+ for (int i=0; i<10000; i++) {
+ double x = ((generator.nextDouble() * 16.0) - 8.0) * generator.nextDouble();
+
+ double tst = FastMath.tanh(x);
+ double ref = tanh(field.newDfp(x)).toDouble();
+ double err = (tst - ref) / ref;
+
+ if (err != 0) {
+ double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref) ^ 1)));
+ double errulp = field.newDfp(tst).subtract(tanh(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble();
+ //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp);
+ maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+ }
+ }
+
+ Assert.assertTrue("tanh() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
+ }
+
@Test
public void testCbrtAccuracy() {
double maxerrulp = 0.0;
@@ -1017,7 +1095,49 @@ public class FastMathTest {
time = System.currentTimeMillis() - time;
System.out.println("FastMath.cbrt " + time + "\t" + x);
- x = 0;
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
+ x += StrictMath.cosh(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.print("StrictMath.cosh " + time + "\t" + x + "\t");
+
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
+ x += FastMath.cosh(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.println("FastMath.cosh " + time + "\t" + x);
+
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
+ x += StrictMath.sinh(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.print("StrictMath.sinh " + time + "\t" + x + "\t");
+
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
+ x += FastMath.sinh(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.println("FastMath.sinh " + time + "\t" + x);
+
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
+ x += StrictMath.tanh(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.print("StrictMath.tanh " + time + "\t" + x + "\t");
+
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
+ x += FastMath.tanh(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.println("FastMath.tanh " + time + "\t" + x);
+
+ x = 0;
time = System.currentTimeMillis();
for (int i = 0; i < numberOfRuns; i++)
x += StrictMath.expm1(-i / 100000.0);
@@ -1037,6 +1157,7 @@ public class FastMathTest {
x += FastMath.expm1(-i / 100000.0);
time = System.currentTimeMillis() - time;
System.out.println("FastMath.expm1 " + time + "\t" + x);
+
}
}