date and time created 85/05/23 17:12:33 by miriam

SCCS-vsn: old/libm/libm/IEEE/cbrt.c 1.1

diff --git a/usr/src/old/libm/libm/IEEE/cbrt.c b/usr/src/old/libm/libm/IEEE/cbrt.c
new file mode 100644 (file)
index 0000000..5901b89
--- /dev/null
+++ b/usr/src/old/libm/libm/IEEE/cbrt.c
@@ -0,0 +1,98 @@
+From Prof. Kahan at UC at Berkeley
+/* 
+ * Copyright (c) 1985 Regents of the University of California.
+ * 
+ * Use and reproduction of this software are granted  in  accordance  with
+ * the terms and conditions specified in  the  Berkeley  Software  License
+ * Agreement (in particular, this entails acknowledgement of the programs'
+ * source, and inclusion of this notice) with the additional understanding
+ * that  all  recipients  should regard themselves as participants  in  an
+ * ongoing  research  project and hence should  feel  obligated  to report
+ * their  experiences (good or bad) with these elementary function  codes,
+ * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
+ */
+
+#ifndef lint
+static char sccsid[] = "@(#)cbrt.c     1.1 (Berkeley) %G%";
+#endif not lint
+
+/* kahan's cube root (53 bits IEEE double precision)
+ * for IEEE machines only
+ * coded in C by K.C. Ng, 4/30/85
+ *
+ * Accuracy:
+ *     better than 0.667 ulps according to an error analysis. Maximum
+ * error observed was 0.666 ulps in an 1,000,000 random arguments test.
+ *
+ * Warning: this code is semi machine dependent; the ordering of words in
+ * a floating point number must be known in advance. I assume that the
+ * long interger at the address of a floating point number will be the
+ * leading 32 bits of that floating point number (i.e., sign, exponent,
+ * and the 20 most significant bits).
+ * On a National machine, it has different ordering; therefore, this code 
+ * must be compiled with flag -DNATIONAL. 
+ */
+#ifndef VAX
+
+static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
+                    B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
+static double
+           C= 19./35.,
+           D= -864./1225.,
+           E= 99./70.,
+           F= 45./28.,
+           G= 5./14.;
+
+double cbrt(x) 
+double x;
+{
+       double r,s,t=0.0,w;
+       unsigned long *px = (unsigned long *) &x,
+                     *pt = (unsigned long *) &t,
+                     mexp,sign;
+
+#ifdef NATIONAL /* ordering of words in a floating points number */
+       int n0=1,n1=0;
+#else
+       int n0=0,n1=1;
+#endif
+
+       mexp=px[n0]&0x7ff00000;
+       if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
+       if(x==0.0) return(x);           /* cbrt(0) is itself */
+
+       sign=px[n0]&0x80000000; /* sign= sign(x) */
+       px[n0] ^= sign;         /* x=|x| */
+
+
+    /* rough cbrt to 5 bits */
+       if(mexp==0)             /* subnormal number */
+         {pt[n0]=0x43500000;   /* set t= 2**54 */
+          t*=x; pt[n0]=pt[n0]/3+B2;
+         }
+       else
+         pt[n0]=px[n0]/3+B1;   
+
+
+    /* 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) */ 
+       pt[n1]=0; pt[n0]+=0x00000001;
+
+
+    /* 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;
+
+
+    /* retore the sign bit */
+       pt[n0] |= sign;
+       return(t);
+}
+#endif