Return to cbrt.c CVS log

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

Annotation of 43BSDTahoe/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:  * All rights reserved.
                      4:  *
                      5:  * Redistribution and use in source and binary forms are permitted
                      6:  * provided that the above copyright notice and this paragraph are
                      7:  * duplicated in all such forms and that any documentation,
                      8:  * advertising materials, and other materials related to such
                      9:  * distribution and use acknowledge that the software was developed
                     10:  * by the University of California, Berkeley.  The name of the
                     11:  * University may not be used to endorse or promote products derived
                     12:  * from this software without specific prior written permission.
                     13:  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
                     14:  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
                     15:  * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
                     16:  *
                     17:  * All recipients should regard themselves as participants in an ongoing
                     18:  * research project and hence should feel obligated to report their
                     19:  * experiences (good or bad) with these elementary function codes, using
                     20:  * the sendbug(8) program, to the authors.
                     21:  */
                     22: 
                     23: #ifndef lint
                     24: static char sccsid[] = "@(#)cbrt.c     5.4 (Berkeley) 6/30/88";
                     25: #endif /* not lint */
                     26: 
                     27: /* kahan's cube root (53 bits IEEE double precision)
                     28:  * for IEEE machines only
                     29:  * coded in C by K.C. Ng, 4/30/85
                     30:  *
                     31:  * Accuracy:
                     32:  *     better than 0.667 ulps according to an error analysis. Maximum
                     33:  * error observed was 0.666 ulps in an 1,000,000 random arguments test.
                     34:  *
                     35:  * Warning: this code is semi machine dependent; the ordering of words in
                     36:  * a floating point number must be known in advance. I assume that the
                     37:  * long interger at the address of a floating point number will be the
                     38:  * leading 32 bits of that floating point number (i.e., sign, exponent,
                     39:  * and the 20 most significant bits).
                     40:  * On a National machine, it has different ordering; therefore, this code 
                     41:  * must be compiled with flag -DNATIONAL. 
                     42:  */
                     43: #if !defined(vax)&&!defined(tahoe)
                     44: 
                     45: static unsigned long B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
                     46:                     B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
                     47: static double
                     48:            C= 19./35.,
                     49:            D= -864./1225.,
                     50:            E= 99./70.,
                     51:            F= 45./28.,
                     52:            G= 5./14.;
                     53: 
                     54: double cbrt(x) 
                     55: double x;
                     56: {
                     57:        double r,s,t=0.0,w;
                     58:        unsigned long *px = (unsigned long *) &x,
                     59:                      *pt = (unsigned long *) &t,
                     60:                      mexp,sign;
                     61: 
                     62: #ifdef national /* ordering of words in a floating points number */
                     63:        int n0=1,n1=0;
                     64: #else  /* national */
                     65:        int n0=0,n1=1;
                     66: #endif /* national */
                     67: 
                     68:        mexp=px[n0]&0x7ff00000;
                     69:        if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
                     70:        if(x==0.0) return(x);           /* cbrt(0) is itself */
                     71: 
                     72:        sign=px[n0]&0x80000000; /* sign= sign(x) */
                     73:        px[n0] ^= sign;         /* x=|x| */
                     74: 
                     75: 
                     76:     /* rough cbrt to 5 bits */
                     77:        if(mexp==0)             /* subnormal number */
                     78:          {pt[n0]=0x43500000;   /* set t= 2**54 */
                     79:           t*=x; pt[n0]=pt[n0]/3+B2;
                     80:          }
                     81:        else
                     82:          pt[n0]=px[n0]/3+B1;   
                     83: 
                     84: 
                     85:     /* new cbrt to 23 bits, may be implemented in single precision */
                     86:        r=t*t/x;
                     87:        s=C+r*t;
                     88:        t*=G+F/(s+E+D/s);       
                     89: 
                     90:     /* chopped to 20 bits and make it larger than cbrt(x) */ 
                     91:        pt[n1]=0; pt[n0]+=0x00000001;
                     92: 
                     93: 
                     94:     /* one step newton iteration to 53 bits with error less than 0.667 ulps */
                     95:        s=t*t;          /* t*t is exact */
                     96:        r=x/s;
                     97:        w=t+t;
                     98:        r=(r-t)/(w+r);  /* r-t is exact */
                     99:        t=t+t*r;
                    100: 
                    101: 
                    102:     /* retore the sign bit */
                    103:        pt[n0] |= sign;
                    104:        return(t);
                    105: }
                    106: #endif

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