Annotation of src/lib/libm/src/s_cbrt.c, Revision 1.2
1.1 jtc 1:
2: /* @(#)s_cbrt.c 5.1 93/09/24 */
3: /*
4: * ====================================================
5: * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6: *
7: * Developed at SunPro, 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:
1.2 ! jtc 15: #include <math.h>
1.1 jtc 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: double one = 1.0;
47: int n0,hx;
48: double r,s,t=0.0,w;
49: unsigned *pt = (unsigned *) &t, sign;
50:
51:
52: n0 = ((*(int*)&one)>>29)^1; /* index of high word */
53: hx = *( n0 + (int*)&x); /* high word of x */
54: sign=hx&0x80000000; /* sign= sign(x) */
55: hx ^=sign;
56: if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
57: if((hx|*(1-n0+(int*)&x))==0)
58: return(x); /* cbrt(0) is itself */
59:
60: *(n0+(int*)&x) = hx; /* x <- |x| */
61: /* rough cbrt to 5 bits */
62: if(hx<0x00100000) /* subnormal number */
63: {pt[n0]=0x43500000; /* set t= 2**54 */
64: t*=x; pt[n0]=pt[n0]/3+B2;
65: }
66: else
67: pt[n0]=hx/3+B1;
68:
69:
70: /* new cbrt to 23 bits, may be implemented in single precision */
71: r=t*t/x;
72: s=C+r*t;
73: t*=G+F/(s+E+D/s);
74:
75: /* chopped to 20 bits and make it larger than cbrt(x) */
76: pt[1-n0]=0; pt[n0]+=0x00000001;
77:
78:
79: /* one step newton iteration to 53 bits with error less than 0.667 ulps */
80: s=t*t; /* t*t is exact */
81: r=x/s;
82: w=t+t;
83: r=(r-t)/(w+r); /* r-s is exact */
84: t=t+t*r;
85:
86: /* retore the sign bit */
87: pt[n0] |= sign;
88: return(t);
89: }
CVSweb <webmaster@jp.NetBSD.org>