diff --git a/src/main/java/org/apache/commons/math/util/FastMath.java b/src/main/java/org/apache/commons/math/util/FastMath.java index 09c7cd0bf..9c471e636 100644 --- a/src/main/java/org/apache/commons/math/util/FastMath.java +++ b/src/main/java/org/apache/commons/math/util/FastMath.java @@ -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 diff --git a/src/test/java/org/apache/commons/math/util/FastMathTest.java b/src/test/java/org/apache/commons/math/util/FastMathTest.java index c2d8c54a2..46a2ec005 100644 --- a/src/test/java/org/apache/commons/math/util/FastMathTest.java +++ b/src/test/java/org/apache/commons/math/util/FastMathTest.java @@ -787,6 +787,28 @@ public class FastMathTest { Assert.assertTrue("acos() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP); } + @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; @@ -938,6 +960,7 @@ public class FastMathTest { time = System.currentTimeMillis() - time; System.out.println("FastMath.cos " + time + "\t" + x); + x = 0; time = System.currentTimeMillis(); for (int i = 0; i < numberOfRuns; i++) x += StrictMath.acos(i / 10000000.0); @@ -980,6 +1003,20 @@ public class FastMathTest { System.out.println("FastMath.atan " + time + "\t" + x); 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); + + x = 0; time = System.currentTimeMillis(); for (int i = 0; i < numberOfRuns; i++) x += StrictMath.expm1(-i / 100000.0);