4.3BSD - /usr/src/usr.lib/libm/IEEE/cbrt.c

   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

Defined functions

Defined variables