Added comments about the magic behind

    <cbrt(x) in bits> ~= <x in bits>/3 + BIAS.
Keep the large comments only in the double version as usual.

Fixed some style bugs (mainly grammar and spelling errors in comments).

2 files changed, 25 insertions, 14 deletions

diff --git a/lib/msun/src/s_cbrt.c b/lib/msun/src/s_cbrt.c
index 4d70825..4b140c7 100644
--- a/lib/msun/src/s_cbrt.c
+++ b/lib/msun/src/s_cbrt.c
@@ -21,8 +21,8 @@ static char rcsid[] = "$FreeBSD$";
  * Return cube root of x
  */
 static const u_int32_t
-	B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
-	B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
+	B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
+	B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
 
 static const double
 C =  5.42857142857142815906e-01, /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
@@ -48,7 +48,21 @@ cbrt(double x)
 	    return(x);		/* cbrt(0) is itself */
 
 	SET_HIGH_WORD(x,hx);	/* x <- |x| */
-    /* rough cbrt to 5 bits */
+    /*
+     * Rough cbrt to 5 bits:
+     *    cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
+     * where e is integral and >= 0, m is real and in [0, 1), and "/" and
+     * "%" are integer division and modulus with rounding towards minus
+     * infinity.  The RHS is always >= the LHS and has a maximum relative
+     * error of about 1 in 16.  Adding a bias of -0.03306235651 to the
+     * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
+     * floating point representation, for finite positive normal values,
+     * ordinary integer divison of the value in bits magically gives
+     * almost exactly the RHS of the above provided we first subtract the
+     * exponent bias (1023 for doubles) and later add it back.  We do the
+     * subtraction virtually to keep e >= 0 so that ordinary integer
+     * division rounds towards minus infinity; this is also efficient.
+     */
 	if(hx<0x00100000) 		/* subnormal number */
 	  {SET_HIGH_WORD(t,0x43500000);	/* set t= 2**54 */
 	   t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);
@@ -56,25 +70,23 @@ cbrt(double x)
 	else
 	  SET_HIGH_WORD(t,hx/3+B1);
 
-
-    /* new cbrt to 23 bits, may be implemented in single precision */
+    /* new cbrt to 23 bits; may be implemented in single precision */
 	r=t*t/x;
 	s=C+r*t;
 	t*=G+F/(s+E+D/s);
 
-    /* chopped to 20 bits and make it larger than cbrt(x) */
+    /* chop t to 20 bits and make it larger than cbrt(x) */
 	GET_HIGH_WORD(high,t);
 	INSERT_WORDS(t,high+0x00000001,0);
 
-
-    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
+    /* 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 */
+	r=(r-t)/(w+r);	/* r-t is exact */
 	t=t+t*r;
 
-    /* retore the sign bit */
+    /* restore the sign bit */
 	GET_HIGH_WORD(high,t);
 	SET_HIGH_WORD(t,high|sign);
 	return(t);
diff --git a/lib/msun/src/s_cbrtf.c b/lib/msun/src/s_cbrtf.c
index 1453ac2..6de11e5 100644
--- a/lib/msun/src/s_cbrtf.c
+++ b/lib/msun/src/s_cbrtf.c
@@ -25,8 +25,8 @@ static char rcsid[] = "$FreeBSD$";
  * Return cube root of x
  */
 static const unsigned
-	B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */
-	B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */
+	B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
+	B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
 
 static const float
 C =  5.4285717010e-01, /* 19/35     = 0x3f0af8b0 */
@@ -59,7 +59,6 @@ cbrtf(float x)
 	else
 	  SET_FLOAT_WORD(t,hx/3+B1);
 
-
     /* new cbrt to 23 bits */
 	r=t*t/x;
 	s=C+r*t;
@@ -76,7 +75,7 @@ cbrtf(float x)
 	r=(r-t)/(w+r);	/* r-t is exact */
 	t=t+t*r;
 
-    /* retore the sign bit */
+    /* restore the sign bit */
 	GET_FLOAT_WORD(high,t);
 	SET_FLOAT_WORD(t,high|sign);
 	return(t);