Return to cbrt.c CVS log

Up to [CSRG BSD Unix] / 43BSD / usr.lib / libm / IEEE

Annotation of 43BSD/usr.lib/libm/IEEE/cbrt.c, revision 1.1.1.1

1.1       root        1: /* 
                      2:  * Copyright (c) 1985 Regents of the University of California.
                      3:  * 
                      4:  * Use and reproduction of this software are granted  in  accordance  with
                      5:  * the terms and conditions specified in  the  Berkeley  Software  License
                      6:  * Agreement (in particular, this entails acknowledgement of the programs'
                      7:  * source, and inclusion of this notice) with the additional understanding
                      8:  * that  all  recipients  should regard themselves as participants  in  an
                      9:  * ongoing  research  project and hence should  feel  obligated  to report
                     10:  * their  experiences (good or bad) with these elementary function  codes,
                     11:  * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
                     12:  */
                     13: 
                     14: #ifndef lint
                     15: static char sccsid[] = "@(#)cbrt.c     1.1 (Berkeley) 5/23/85";
                     16: #endif not lint
                     17: 
                     18: /* kahan's cube root (53 bits IEEE double precision)
                     19:  * for IEEE machines only
                     20:  * coded in C by K.C. Ng, 4/30/85
                     21:  *
                     22:  * Accuracy:
                     23:  *     better than 0.667 ulps according to an error analysis. Maximum
                     24:  * error observed was 0.666 ulps in an 1,000,000 random arguments test.
                     25:  *
                     26:  * Warning: this code is semi machine dependent; the ordering of words in
                     27:  * a floating point number must be known in advance. I assume that the
                     28:  * long interger at the address of a floating point number will be the
                     29:  * leading 32 bits of that floating point number (i.e., sign, exponent,
                     30:  * and the 20 most significant bits).
                     31:  * On a National machine, it has different ordering; therefore, this code 
                     32:  * must be compiled with flag -DNATIONAL. 
                     33:  */
                     34: #ifndef VAX
                     35: 
                     36: static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
                     37:                     B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
                     38: static double
                     39:            C= 19./35.,
                     40:            D= -864./1225.,
                     41:            E= 99./70.,
                     42:            F= 45./28.,
                     43:            G= 5./14.;
                     44: 
                     45: double cbrt(x) 
                     46: double x;
                     47: {
                     48:        double r,s,t=0.0,w;
                     49:        unsigned long *px = (unsigned long *) &x,
                     50:                      *pt = (unsigned long *) &t,
                     51:                      mexp,sign;
                     52: 
                     53: #ifdef NATIONAL /* ordering of words in a floating points number */
                     54:        int n0=1,n1=0;
                     55: #else
                     56:        int n0=0,n1=1;
                     57: #endif
                     58: 
                     59:        mexp=px[n0]&0x7ff00000;
                     60:        if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
                     61:        if(x==0.0) return(x);           /* cbrt(0) is itself */
                     62: 
                     63:        sign=px[n0]&0x80000000; /* sign= sign(x) */
                     64:        px[n0] ^= sign;         /* x=|x| */
                     65: 
                     66: 
                     67:     /* rough cbrt to 5 bits */
                     68:        if(mexp==0)             /* subnormal number */
                     69:          {pt[n0]=0x43500000;   /* set t= 2**54 */
                     70:           t*=x; pt[n0]=pt[n0]/3+B2;
                     71:          }
                     72:        else
                     73:          pt[n0]=px[n0]/3+B1;   
                     74: 
                     75: 
                     76:     /* new cbrt to 23 bits, may be implemented in single precision */
                     77:        r=t*t/x;
                     78:        s=C+r*t;
                     79:        t*=G+F/(s+E+D/s);       
                     80: 
                     81:     /* chopped to 20 bits and make it larger than cbrt(x) */ 
                     82:        pt[n1]=0; pt[n0]+=0x00000001;
                     83: 
                     84: 
                     85:     /* one step newton iteration to 53 bits with error less than 0.667 ulps */
                     86:        s=t*t;          /* t*t is exact */
                     87:        r=x/s;
                     88:        w=t+t;
                     89:        r=(r-t)/(w+r);  /* r-s is exact */
                     90:        t=t+t*r;
                     91: 
                     92: 
                     93:     /* retore the sign bit */
                     94:        pt[n0] |= sign;
                     95:        return(t);
                     96: }
                     97: #endif

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.