Hi Tom,
Tom Tromey wrote:
>>>>>> "Carsten" == Carsten Neumann <[EMAIL PROTECTED]> writes:
>
> Carsten> + int[] words = new int[2];
>
> I think cbrt should not require allocation. It shouldn't be too hard
> to rewrite it without this... could you do that?
>
Here is an updated version.
Cheers,
Carsten
PS: If accepted, please commit it for me, thanks.
2006-07-13 Carsten Neumann <[EMAIL PROTECTED]>
* StrictMath.java (cbrt): New method.
(CBRT_B1): New field.
(CBRT_B2): Likewise.
(CBRT_C): Likewise.
(CBRT_D): Likewise.
(CBRT_E): Likewise.
(CBRT_F): Likewise.
(CBRT_G): 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 19:31:28 -0000
@@ -633,6 +633,111 @@
}
/**
+ * 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;
+
+ long bits;
+ int lowWord;
+ int highWord;
+
+ 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;
+ }
+
+ bits = Double.doubleToLongBits(x);
+
+ if (bits < 0x0010000000000000L) // subnormal number
+ {
+ t = TWO_54;
+ t *= x;
+
+ // __HI(t)=__HI(t)/3+B2;
+ bits = Double.doubleToLongBits(t);
+ highWord = (int) (bits >> 32);
+ lowWord = (int) (bits & 0x00000000ffffffffL);
+
+ highWord = highWord / 3 + CBRT_B2;
+ bits = ((long) highWord << 32) | ((long) lowWord);
+ t = Double.longBitsToDouble(bits);
+ }
+ else
+ {
+ // __HI(t)=__HI(x)/3+B1;
+ bits = Double.doubleToLongBits(x);
+ highWord = (int) (bits >> 32);
+ lowWord = (int) (bits & 0x00000000ffffffffL);
+
+ highWord = highWord / 3 + CBRT_B1;
+ bits = ((long) highWord << 32) | ((long) lowWord);
+ t = Double.longBitsToDouble(bits);
+ }
+
+ // 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)
+ bits = Double.doubleToLongBits(x);
+ highWord = (int) (bits >> 32);
+
+ // __LO(t)=0;
+ // __HI(t)+=0x00000001;
+ lowWord = 0;
+ highWord += 1;
+
+ bits = ((long) highWord << 32) | ((long) lowWord);
+ t = Double.longBitsToDouble(bits);
+
+ // 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 +1534,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].
*
2006-07-13 Carsten Neumann <[EMAIL PROTECTED]>
* StrictMath.java (cbrt): New method.
(CBRT_B1): New field.
(CBRT_B2): Likewise.
(CBRT_C): Likewise.
(CBRT_D): Likewise.
(CBRT_E): Likewise.
(CBRT_F): Likewise.
(CBRT_G): Likewise.
