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
15 /*
16 FUNCTION
17 <<cbrt>>, <<cbrtf>>---cube root
18
19 INDEX
20 cbrt
21 INDEX
22 cbrtf
23
24 SYNOPSIS
25 #include <math.h>
26 double cbrt(double <[x]>);
27 float cbrtf(float <[x]>);
28
29 DESCRIPTION
30 <<cbrt>> computes the cube root of the argument.
31
32 RETURNS
33 The cube root is returned.
34
35 PORTABILITY
36 <<cbrt>> is in System V release 4. <<cbrtf>> is an extension.
37 */
38
39 #include "fdlibm.h"
40
41 #ifdef _NEED_FLOAT64
42
43 /* cbrt(x)
44 * Return cube root of x
45 */
46 static const __uint32_t
47 B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
48 B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
49
50 static const __float64
51 C = _F_64(5.42857142857142815906e-01), /* 19/35 = 0x3FE15F15, 0xF15F15F1 */
52 D = _F_64(-7.05306122448979611050e-01), /* -864/1225 = 0xBFE691DE, 0x2532C834 */
53 E = _F_64(1.41428571428571436819e+00), /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */
54 F = _F_64(1.60714285714285720630e+00), /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */
55 G = _F_64(3.57142857142857150787e-01); /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
56
57 __float64
cbrt64(__float64 x)58 cbrt64(__float64 x)
59 {
60 __int32_t hx;
61 __float64 r,s,t=_F_64(0.0),w;
62 __uint32_t sign;
63 __uint32_t high,low;
64
65 GET_HIGH_WORD(hx,x);
66 sign=hx&0x80000000; /* sign= sign(x) */
67 hx ^=sign;
68 if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
69 GET_LOW_WORD(low,x);
70 if((hx|low)==0)
71 return(x); /* cbrt(0) is itself */
72
73 SET_HIGH_WORD(x,hx); /* x <- |x| */
74 /* rough cbrt to 5 bits */
75 if(hx<0x00100000) /* subnormal number */
76 {SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
77 t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);
78 }
79 else
80 SET_HIGH_WORD(t,hx/3+B1);
81
82
83 /* new cbrt to 23 bits, may be implemented in single precision */
84 r=t*t/x;
85 s=C+r*t;
86 t*=G+F/(s+E+D/s);
87
88 /* chopped to 20 bits and make it larger than cbrt(x) */
89 GET_HIGH_WORD(high,t);
90 INSERT_WORDS(t,high+0x00000001,0);
91
92
93 /* one step newton iteration to 53 bits with error less than 0.667 ulps */
94 s=t*t; /* t*t is exact */
95 r=x/s;
96 w=t+t;
97 r=(r-t)/(w+r); /* r-s is exact */
98 t=t+t*r;
99
100 /* retore the sign bit */
101 GET_HIGH_WORD(high,t);
102 SET_HIGH_WORD(t,high|sign);
103 return(t);
104 }
105
106 _MATH_ALIAS_d_d(cbrt)
107
108 #endif /* _NEED_FLOAT64 */
109
