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/09/11 18:57:30 UTC
svn commit: r996172 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/util/FastMath.java
test/java/org/apache/commons/math/util/FastMathTest.java
Author: luc
Date: Sat Sep 11 16:57:30 2010
New Revision: 996172
URL: http://svn.apache.org/viewvc?rev=996172&view=rev
Log:
added fast cubic root computation
JIRA: MATH-375
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/FastMathTest.java
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java?rev=996172&r1=996171&r2=996172&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/FastMath.java Sat Sep 11 16:57:30 2010
@@ -154,6 +154,13 @@ public class FastMath {
*/
private static final double EIGHTHES[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
+ /* Table of 2^((n+2)/3) */
+ private static final double CBRTTWO[] = { 0.6299605249474366,
+ 0.7937005259840998,
+ 1.0,
+ 1.2599210498948732,
+ 1.5874010519681994 };
+
// Initialize tables
static {
int i;
@@ -212,14 +219,6 @@ public class FastMath {
return Math.sqrt(a);
}
- /** Compute the cubic root of a number.
- * @param a number on which evaluation is done
- * @return cubic root of a
- */
- public static double cbrt(final double a) {
- return Math.cbrt(a);
- }
-
/** Compute the hyperbolic cosine of a number.
* @param a number on which evaluation is done
* @return hyperbolic cosine of a
@@ -1020,11 +1019,6 @@ public class FastMath {
ya = aa + tmp - tmp;
yb = aa - ya + ab;
- if (hiPrec != null) {
- hiPrec[0] = ya;
- hiPrec[1] = yb;
- }
-
return ya + yb;
}
}
@@ -2886,6 +2880,90 @@ public class FastMath {
return atan(ra, rb, x<0);
}
+ /** Compute the cubic root of a number.
+ * @param a number on which evaluation is done
+ * @return cubic root of a
+ */
+ public static double cbrt(double x) {
+ /* Convert input double to bits */
+ long inbits = Double.doubleToLongBits(x);
+ int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
+ boolean subnormal = false;
+
+ if (exponent == -1023) {
+ if (x == 0) {
+ return x;
+ }
+
+ /* Subnormal, so normalize */
+ subnormal = true;
+ x *= 1.8014398509481984E16; // 2^54
+ inbits = Double.doubleToLongBits(x);
+ exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
+ }
+
+ if (exponent == 1024) {
+ // Nan or infinity. Don't care which.
+ return x;
+ }
+
+ /* Divide the exponent by 3 */
+ int exp3 = exponent / 3;
+
+ /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
+ double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
+ (long)(((exp3 + 1023) & 0x7ff)) << 52);
+
+ /* This will be a number between 1 and 2 */
+ final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
+
+ /* Estimate the cube root of mant by polynomial */
+ double est = -0.010714690733195933;
+ est = est * mant + 0.0875862700108075;
+ est = est * mant + -0.3058015757857271;
+ est = est * mant + 0.7249995199969751;
+ est = est * mant + 0.5039018405998233;
+
+ est *= CBRTTWO[exponent % 3 + 2];
+
+ // est should now be good to about 15 bits of precision. Do 2 rounds of
+ // Newton's method to get closer, this should get us full double precision
+ // Scale down x for the purpose of doing newtons method. This avoids over/under flows.
+ final double xs = x / (p2*p2*p2);
+ est += (xs - est*est*est) / (3*est*est);
+ est += (xs - est*est*est) / (3*est*est);
+
+ // Do one round of Newton's method in extended precision to get the last bit right.
+ double temp = est * 1073741824.0;
+ double ya = est + temp - temp;
+ double yb = est - ya;
+
+ double za = ya * ya;
+ double zb = ya * yb * 2.0 + yb * yb;
+ temp = za * 1073741824.0;
+ double temp2 = za + temp - temp;
+ zb += (za - temp2);
+ za = temp2;
+
+ zb = za * yb + ya * zb + zb * yb;
+ za = za * ya;
+
+ double na = xs - za;
+ double nb = -(na - xs + za);
+ nb -= zb;
+
+ est += (na+nb)/(3*est*est);
+
+ /* Scale by a power of two, so this is exact. */
+ est *= p2;
+
+ if (subnormal) {
+ est *= 3.814697265625E-6; // 2^-18
+ }
+
+ return est;
+ }
+
/**
* Convert degrees to radians, with error of less than 0.5 ULP
* @param x angle in degrees
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/FastMathTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/FastMathTest.java?rev=996172&r1=996171&r2=996172&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/FastMathTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/FastMathTest.java Sat Sep 11 16:57:30 2010
@@ -788,6 +788,45 @@ public class FastMathTest {
}
@Test
+ public void testCbrtAccuracy() {
+ double maxerrulp = 0.0;
+
+ for (int i=0; i<10000; i++) {
+ double x = ((generator.nextDouble() * 200.0) - 100.0) * generator.nextDouble();
+
+ double tst = FastMath.cbrt(x);
+ double ref = cbrt(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(cbrt(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("cbrt() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
+ }
+
+ private Dfp cbrt(Dfp x) {
+ boolean negative=false;
+
+ if (x.lessThan(field.getZero())) {
+ negative = true;
+ x = x.negate();
+ }
+
+ Dfp y = DfpMath.pow(x, field.getOne().divide(3));
+
+ if (negative) {
+ y = y.negate();
+ }
+
+ return y;
+ }
+
+ @Test
public void testToDegrees() {
double maxerrulp = 0.0;
for (int i = 0; i < NUMBER_OF_TRIALS; i++) {
@@ -966,6 +1005,13 @@ public class FastMathTest {
x = 0;
time = System.currentTimeMillis();
for (int i = 0; i < numberOfRuns; i++)
+ x += StrictMath.cbrt(i / 1000000.0);
+ time = System.currentTimeMillis() - time;
+ System.out.print("StrictMath.cbrt " + time + "\t" + x + "\t");
+
+ x = 0;
+ time = System.currentTimeMillis();
+ for (int i = 0; i < numberOfRuns; i++)
x += FastMath.cbrt(i / 1000000.0);
time = System.currentTimeMillis() - time;
System.out.println("FastMath.cbrt " + time + "\t" + x);