home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / math / s_cbrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-16  |  2.0 KB  |  88 lines
  1.  
  2. /* @(#)s_cbrt.c 1.3 95/01/18 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  8.  * Permission to use, copy, modify, and distribute this
  9.  * software is freely granted, provided that this notice 
  10.  * is preserved.
  11.  * ====================================================
  12.  *
  13.  */
  14.  
  15. #include "fdlibm.h"
  16.  
  17. /* cbrt(x)
  18.  * Return cube root of x
  19.  */
  20. #ifdef __STDC__
  21. static const unsigned 
  22. #else
  23. static unsigned 
  24. #endif
  25.     B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
  26.     B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
  27.  
  28. #ifdef __STDC__
  29. static const double
  30. #else
  31. static double
  32. #endif
  33. C =  5.42857142857142815906e-01, /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
  34. D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
  35. E =  1.41428571428571436819e+00, /* 99/70     = 0x3FF6A0EA, 0x0EA0EA0F */
  36. F =  1.60714285714285720630e+00, /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
  37. G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
  38.  
  39. #ifdef __STDC__
  40.     double cbrt(double x) 
  41. #else
  42.     double cbrt(x) 
  43.     double x;
  44. #endif
  45. {
  46.     int    hx;
  47.     double r,s,t=0.0,w;
  48.     unsigned sign;
  49.  
  50.  
  51.     hx = __HI(x);        /* high word of x */
  52.     sign=hx&0x80000000;         /* sign= sign(x) */
  53.     hx  ^=sign;
  54.     if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
  55.     if((hx|__LO(x))==0) 
  56.         return(x);        /* cbrt(0) is itself */
  57.  
  58.     __HI(x) = hx;    /* x <- |x| */
  59.     /* rough cbrt to 5 bits */
  60.     if(hx<0x00100000)         /* subnormal number */
  61.       {__HI(t)=0x43500000;         /* set t= 2**54 */
  62.        t*=x; __HI(t)=__HI(t)/3+B2;
  63.       }
  64.     else
  65.       __HI(t)=hx/3+B1;    
  66.  
  67.  
  68.     /* new cbrt to 23 bits, may be implemented in single precision */
  69.     r=t*t/x;
  70.     s=C+r*t;
  71.     t*=G+F/(s+E+D/s);    
  72.  
  73.     /* chopped to 20 bits and make it larger than cbrt(x) */ 
  74.     __LO(t)=0; __HI(t)+=0x00000001;
  75.  
  76.  
  77.     /* one step newton iteration to 53 bits with error less than 0.667 ulps */
  78.     s=t*t;        /* t*t is exact */
  79.     r=x/s;
  80.     w=t+t;
  81.     r=(r-t)/(w+r);    /* r-s is exact */
  82.     t=t+t*r;
  83.  
  84.     /* retore the sign bit */
  85.     __HI(t) |= sign;
  86.     return(t);
  87. }
  88.