1 /* @(#)s_cbrt.c 5.1 93/09/24 */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 #ifndef lint 14 static char rcsid[] = "$Id: s_cbrt.c,v 1.4 1994/03/03 17:04:28 jtc Exp $"; 15 #endif 16 17 #include <math.h> 18 #include <machine/endian.h> 19 20 #if BYTE_ORDER == LITTLE_ENDIAN 21 #define n0 1 22 #else 23 #define n0 0 24 #endif 25 26 /* cbrt(x) 27 * Return cube root of x 28 */ 29 #ifdef __STDC__ 30 static const unsigned 31 #else 32 static unsigned 33 #endif 34 B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ 35 B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ 36 37 #ifdef __STDC__ 38 static const double 39 #else 40 static double 41 #endif 42 C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */ 43 D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */ 44 E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */ 45 F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */ 46 G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */ 47 48 #ifdef __STDC__ 49 double cbrt(double x) 50 #else 51 double cbrt(x) 52 double x; 53 #endif 54 { 55 int hx; 56 double r,s,t=0.0,w; 57 unsigned *pt = (unsigned *) &t, sign; 58 59 hx = *( n0 + (int*)&x); /* high word of x */ 60 sign=hx&0x80000000; /* sign= sign(x) */ 61 hx ^=sign; 62 if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */ 63 if((hx|*(1-n0+(int*)&x))==0) 64 return(x); /* cbrt(0) is itself */ 65 66 *(n0+(int*)&x) = hx; /* x <- |x| */ 67 /* rough cbrt to 5 bits */ 68 if(hx<0x00100000) /* 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]=hx/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[1-n0]=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 /* retore the sign bit */ 93 pt[n0] |= sign; 94 return(t); 95 } 96