Return to cbrt.c CVS log

Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / ieee

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

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