Hi,
the attached patch implements java.lang.StrictMath#cbrt, but I'm
actually quite ashamed to post it: I'm not very knowledgeable about
floating point arithmetic, so I just followed the implementation in
fdlibm (www.netlib.org/fdlibm). Unfortunately this library does some
weird manipulations of the bits of a double, which I simply mimicked in
java.
Cheers,
Carsten
2006-07-13 Carsten Neumann <[EMAIL PROTECTED]>
* StrictMath.java (cbrt): New public method.
(doubleToLowWord): New helper method.
(doubleToHighWord): Likewise.
(doubleToWords): Likewise.
(wordsToDouble): Likewise.
(wordsToDouble): Likewise.
Index: java/lang/StrictMath.java
===================================================================
RCS file: /sources/classpath/classpath/java/lang/StrictMath.java,v
retrieving revision 1.8
diff -u -r1.8 StrictMath.java
--- java/lang/StrictMath.java 15 Mar 2006 21:57:44 -0000 1.8
+++ java/lang/StrictMath.java 13 Jul 2006 17:29:15 -0000
@@ -633,6 +633,147 @@
}
/**
+ * Helper method to obtain the lower 32 bits of a double in IEEE 754
+ * representation.
+ */
+ private static int doubleToLowWord(double x)
+ {
+ return (int) (Double.doubleToLongBits(x) & 0x00000000ffffffffL);
+ }
+
+ /**
+ * Helper method to obtain the upper 32 bits of a double in IEEE 754
+ * representation.
+ */
+ private static int doubleToHighWord(double x)
+ {
+ return (int) ((Double.doubleToLongBits(x) & 0xffffffff00000000L ) >> 32);
+ }
+
+ /**
+ * Helper method to obtain the upper and lower 32 bits of a double
+ * in IEEE 754 representation.
+ */
+ private static void doubleToWords(double x, int[] words)
+ {
+ long longBits = Double.doubleToLongBits(x);
+
+ words[0] = (int) (longBits & 0x00000000ffffffffL);
+ words[1] = (int) ((longBits & 0xffffffff00000000L) >> 32);
+ }
+
+ /**
+ * Helper method to construct a double from the upper and lower 32
+ * bits in words[1] and words[0] resprectively.
+ */
+ private static double wordsToDouble(int[] words)
+ {
+ return wordsToDouble(words[0], words[1]);
+ }
+
+ /**
+ * Helper method to construct a double from the upper and lower 32
+ * bits in highWord and lowWord resprectively.
+ */
+ private static double wordsToDouble(int lowWord, int highWord)
+ {
+ return Double.longBitsToDouble((((long) highWord) << 32) | ((long) lowWord));
+ }
+
+ /**
+ * Returns the cube root of <code>x</code>. The sign of the cube root
+ * is equal to the sign of the sign of <code>x</code>.
+ * Special cases:
+ * <ul>
+ * <li>If the argument is NaN, the result is NaN</li>
+ * <li>If the argument is positive infinity, the result is positive
+ * infinity.</li>
+ * <li>If the argument is negative infinity, the result is negative
+ * infinity.</li>
+ * <li>If the argument is zero, the result is zero with the same
+ * sign as the argument.</li>
+ * </ul>
+ *
+ * @param x the number to take the cube root of
+ * @return the cube root of <code>x</code>
+ * @see #sqrt(double)
+ * @since 1.5
+ */
+ public static double cbrt(double x)
+ {
+ boolean negative = (x < 0);
+ double r;
+ double s;
+ double t = 0.0;
+ double w;
+ int hx;
+ int[] words = new int[2];
+
+ if (negative)
+ {
+ x = -x;
+ }
+
+ // handle the special cases
+ if (x != x)
+ {
+ return Double.NaN;
+ }
+ else if (x == Double.POSITIVE_INFINITY)
+ {
+ return negative ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
+ }
+ else if (x == 0)
+ {
+ return negative ? - 0.0 : 0.0;
+ }
+
+ hx = doubleToHighWord(x);
+
+ if (hx < 0x00100000) // subnormal number
+ {
+ t = TWO_54;
+ t *= x;
+
+ // __HI(t)=__HI(t)/3+B2;
+ doubleToWords(t, words);
+ words[1] = words[1] / 3 + CBRT_B2;
+ t = wordsToDouble(words);
+ }
+ else
+ {
+ // __HI(t)=hx/3+B1;
+ doubleToWords(t, words);
+ words[1] = hx / 3 + CBRT_B1;
+ t = wordsToDouble(words);
+ }
+
+ // new cbrt to 23 bits
+ r = t * t / x;
+ s = CBRT_C + r * t;
+ t *= CBRT_G + CBRT_F / (s + CBRT_E + CBRT_D / s);
+
+ // chopped to 20 bits and make it larger than cbrt(x)
+ doubleToWords(t, words);
+
+ // __LO(t)=0;
+ // __HI(t)+=0x00000001;
+ words[0] = 0;
+ words[1] += 1;
+
+ t = wordsToDouble(words);
+
+ // one step newton iteration to 53 bits with error less than 0.667 ulps
+ s = t * t; // t * t is exact
+ r = x / s;
+ w = t + t;
+ r = (r - t) / (w + r); // r - s is exact
+ t = t + t * r;
+
+ return negative ? -t : t;
+ }
+
+ /**
* Take <em>e</em><sup>a</sup>. The opposite of <code>log()</code>. If the
* argument is NaN, the result is NaN; if the argument is positive infinity,
* the result is positive infinity; and if the argument is negative
@@ -1429,6 +1570,23 @@
AT10 = 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L.
/**
+ * Constants for computing [EMAIL PROTECTED] #cbrt(double)}.
+ */
+ private static final int
+ CBRT_B1 = 715094163, // B1 = (682-0.03306235651)*2**20
+ CBRT_B2 = 696219795; // B2 = (664-0.03306235651)*2**20
+
+ /**
+ * Constants for computing [EMAIL PROTECTED] #cbrt(double)}.
+ */
+ private static final double
+ CBRT_C = 5.42857142857142815906e-01, // Long bits 0x3fe15f15f15f15f1L
+ CBRT_D = -7.05306122448979611050e-01, // Long bits 0xbfe691de2532c834L
+ CBRT_E = 1.41428571428571436819e+00, // Long bits 0x3ff6a0ea0ea0ea0fL
+ CBRT_F = 1.60714285714285720630e+00, // Long bits 0x3ff9b6db6db6db6eL
+ CBRT_G = 3.57142857142857150787e-01; // Long bits 0x3fd6db6db6db6db7L
+
+ /**
* Helper function for reducing an angle to a multiple of pi/2 within
* [-pi/4, pi/4].
*
