[Commit] nickle ChangeLog,1.139,1.140 math.5c,1.42,1.43

[Keith Packard]

Committed by: keithp

Update of /local/src/CVS/nickle
In directory home.keithp.com:/tmp/cvs-serv15686

Modified Files:
	ChangeLog math.5c 
Log Message:
2005-12-01  Keith Packard  <keithp at keithp.com>

	* math.5c:
	cbrt must use more intermediate precision to hit the
	specified error bound


Index: ChangeLog
===================================================================
RCS file: /local/src/CVS/nickle/ChangeLog,v
retrieving revision 1.139
retrieving revision 1.140
diff -u -d -r1.139 -r1.140
--- ChangeLog	20 Nov 2005 08:49:33 -0000	1.139
+++ ChangeLog	1 Dec 2005 16:15:31 -0000	1.140
@@ -1,3 +1,9 @@
+2005-12-01  Keith Packard  <keithp at keithp.com>
+
+	* math.5c:
+	cbrt must use more intermediate precision to hit the
+	specified error bound
+
 2005-11-20  Bart Massey <bart at cs.pdx.edu>
 
         * process.5c:

Index: math.5c
===================================================================
RCS file: /local/src/CVS/nickle/math.5c,v
retrieving revision 1.42
retrieving revision 1.43
diff -u -d -r1.42 -r1.43
--- math.5c	5 Jun 2005 01:03:42 -0000	1.42
+++ math.5c	1 Dec 2005 16:15:31 -0000	1.43
@@ -64,23 +64,22 @@
     {
 	real real_cbrt(real v)
 	{
-	    real	err;
 	    real	prev, cur;
-	    int	n;
-	    int	s;
-
-	    s = sign (v);
+	    
+	    int s = sign (v);
 	    v = imprecise (abs (v));
-	    err = 1/(2**(precision (v)+3));
-	    prev = imprecise (1 / (0.75 * 2**(floor (exponent(v)/3))), precision(v));
-	    for (;;)
-	    {
-		cur = 1/3 * (2 * prev + v / prev**2);
-		if (abs (cur - prev) < err)
-		    break;
+	    int result_bits = precision (v);
+	    int intermediate_bits = result_bits + 10;
+	    v = imprecise (v, intermediate_bits);
+	    real err = imprecise (1/(2**(result_bits+3)), intermediate_bits);
+
+	    cur = imprecise (1 / (0.75 * 2**(floor (exponent(v)/3))),
+			     intermediate_bits);
+	    do {
 		prev = cur;
-	    }
-	    return s * abs (cur);
+		cur = 1/3 * (2 * prev + v / prev**2);
+	    } while (abs (cur - prev) > err);
+	    return s * imprecise (abs (cur), result_bits);
 	}
 
 	if (is_rational (v))