160484Sobrien/* 260484Sobrien * ==================================================== 360484Sobrien * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 460484Sobrien * 560484Sobrien * Developed at SunPro, a Sun Microsystems, Inc. business. 660484Sobrien * Permission to use, copy, modify, and distribute this 760484Sobrien * software is freely granted, provided that this notice 860484Sobrien * is preserved. 960484Sobrien * ==================================================== 1077298Sobrien*/ 1177298Sobrien 1260484Sobrien#ifndef __vax__ 1360484Sobrienstatic const unsigned long 1460484Sobrien B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ 1560484Sobrien B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ 1660484Sobrien 1760484Sobrienstatic const double 1860484Sobrien C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */ 1960484Sobrien D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */ 2060484Sobrien E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */ 2160484Sobrien F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */ 2260484Sobrien G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */ 2360484Sobrien 2460484Sobriendouble 2560484Sobriencbrtl (double x) 2660484Sobrien{ 2760484Sobrien long hx; 2860484Sobrien double r,s,w; 2960484Sobrien double lt; 3060484Sobrien unsigned sign; 3160484Sobrien typedef unsigned unsigned32 __attribute__((mode(SI))); 3260484Sobrien union { 3360484Sobrien double t; 3460484Sobrien unsigned32 pt[2]; 3577298Sobrien } ut, ux; 3677298Sobrien int n0; 37218822Sdim 38218822Sdim ut.t = 1.0; 39218822Sdim n0 = (ut.pt[0] == 0); 40218822Sdim 4177298Sobrien ut.t = 0.0; 4277298Sobrien ux.t = x; 43218822Sdim 44218822Sdim hx = ux.pt[n0]; /* high word of x */ 4560484Sobrien sign=hx&0x80000000; /* sign= sign(x) */ 4660484Sobrien hx ^=sign; 4760484Sobrien if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */ 4860484Sobrien if((hx| ux.pt[1-n0])==0) 49104834Sobrien return(ux.t); /* cbrt(0) is itself */ 50104834Sobrien 5160484Sobrien ux.pt[n0] = hx; 5260484Sobrien /* rough cbrt to 5 bits */ 53104834Sobrien if(hx<0x00100000) /* subnormal number */ 54104834Sobrien {ut.pt[n0]=0x43500000; /* set t= 2**54 */ 5560484Sobrien ut.t*=x; ut.pt[n0]=ut.pt[n0]/3+B2; 5660484Sobrien } 5760484Sobrien else 5860484Sobrien ut.pt[n0]=hx/3+B1; 5960484Sobrien 6060484Sobrien /* new cbrt to 23 bits, may be implemented in single precision */ 6160484Sobrien r=ut.t*ut.t/ux.t; 6260484Sobrien s=C+r*ut.t; 6360484Sobrien ut.t*=G+F/(s+E+D/s); 6460484Sobrien 6560484Sobrien /* chopped to 20 bits and make it larger than cbrt(x) */ 6660484Sobrien ut.pt[1-n0]=0; ut.pt[n0]+=0x00000001; 6777298Sobrien 6877298Sobrien /* one step newton iteration to 53 bits with error less than 0.667 ulps */ 69130561Sobrien s=ut.t*ut.t; /* t*t is exact */ 70130561Sobrien r=ux.t/s; 71218822Sdim w=ut.t+ut.t; 72218822Sdim r=(r-ut.t)/(w+r); /* r-s is exact */ 7360484Sobrien ut.t=ut.t+ut.t*r; 7460484Sobrien 7577298Sobrien /* restore the sign bit */ 7677298Sobrien ut.pt[n0] |= sign; 7760484Sobrien 7860484Sobrien lt = ut.t; 7960484Sobrien lt -= (lt - (x/(lt*lt))) * 0.333333333333333333333; 8060484Sobrien return lt; 81218822Sdim} 82218822Sdim 8360484Sobrienmain () 8460484Sobrien{ 8589857Sobrien if ((int) (cbrtl (27.0) + 0.5) != 3) 8689857Sobrien abort (); 8760484Sobrien 8860484Sobrien exit (0); 8960484Sobrien} 9060484Sobrien#else 91130561Sobrienmain () { exit (0); } 92130561Sobrien#endif 9360484Sobrien