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 /* ieee_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 ieee_cbrt(double x) 41 #else 42 double ieee_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); /* ieee_cbrt(NaN,INF) is itself */ 55 if((hx|__LO(x))==0) 56 return(x); /* ieee_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 ieee_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
