986
|
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 |
}
|