1 | /* sf_cbrt.c -- float version of s_cbrt.c. |
---|
2 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
---|
3 | */ |
---|
4 | |
---|
5 | /* |
---|
6 | * ==================================================== |
---|
7 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
---|
8 | * |
---|
9 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
---|
10 | * Permission to use, copy, modify, and distribute this |
---|
11 | * software is freely granted, provided that this notice |
---|
12 | * is preserved. |
---|
13 | * ==================================================== |
---|
14 | * |
---|
15 | */ |
---|
16 | |
---|
17 | #include "fdlibm.h" |
---|
18 | |
---|
19 | /* cbrtf(x) |
---|
20 | * Return cube root of x |
---|
21 | */ |
---|
22 | #ifdef __STDC__ |
---|
23 | static const __uint32_t |
---|
24 | #else |
---|
25 | static __uint32_t |
---|
26 | #endif |
---|
27 | B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */ |
---|
28 | B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */ |
---|
29 | |
---|
30 | #ifdef __STDC__ |
---|
31 | static const float |
---|
32 | #else |
---|
33 | static float |
---|
34 | #endif |
---|
35 | C = 5.4285717010e-01, /* 19/35 = 0x3f0af8b0 */ |
---|
36 | D = -7.0530611277e-01, /* -864/1225 = 0xbf348ef1 */ |
---|
37 | E = 1.4142856598e+00, /* 99/70 = 0x3fb50750 */ |
---|
38 | F = 1.6071428061e+00, /* 45/28 = 0x3fcdb6db */ |
---|
39 | G = 3.5714286566e-01; /* 5/14 = 0x3eb6db6e */ |
---|
40 | |
---|
41 | #ifdef __STDC__ |
---|
42 | float cbrtf(float x) |
---|
43 | #else |
---|
44 | float cbrtf(x) |
---|
45 | float x; |
---|
46 | #endif |
---|
47 | { |
---|
48 | __int32_t hx; |
---|
49 | float r,s,t; |
---|
50 | __uint32_t sign; |
---|
51 | __uint32_t high; |
---|
52 | |
---|
53 | GET_FLOAT_WORD(hx,x); |
---|
54 | sign=hx&0x80000000; /* sign= sign(x) */ |
---|
55 | hx ^=sign; |
---|
56 | if(!FLT_UWORD_IS_FINITE(hx)) |
---|
57 | return(x+x); /* cbrt(NaN,INF) is itself */ |
---|
58 | if(FLT_UWORD_IS_ZERO(hx)) |
---|
59 | return(x); /* cbrt(0) is itself */ |
---|
60 | |
---|
61 | SET_FLOAT_WORD(x,hx); /* x <- |x| */ |
---|
62 | /* rough cbrt to 5 bits */ |
---|
63 | if(FLT_UWORD_IS_SUBNORMAL(hx)) /* subnormal number */ |
---|
64 | {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */ |
---|
65 | t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2); |
---|
66 | } |
---|
67 | else |
---|
68 | SET_FLOAT_WORD(t,hx/3+B1); |
---|
69 | |
---|
70 | |
---|
71 | /* new cbrt to 23 bits */ |
---|
72 | r=t*t/x; |
---|
73 | s=C+r*t; |
---|
74 | t*=G+F/(s+E+D/s); |
---|
75 | |
---|
76 | /* retore the sign bit */ |
---|
77 | GET_FLOAT_WORD(high,t); |
---|
78 | SET_FLOAT_WORD(t,high|sign); |
---|
79 | return(t); |
---|
80 | } |
---|
81 | |
---|
82 | #ifdef _DOUBLE_IS_32BITS |
---|
83 | |
---|
84 | #ifdef __STDC__ |
---|
85 | double cbrt(double x) |
---|
86 | #else |
---|
87 | double cbrt(x) |
---|
88 | double x; |
---|
89 | #endif |
---|
90 | { |
---|
91 | return (double) cbrtf((float) x); |
---|
92 | } |
---|
93 | |
---|
94 | #endif /* defined(_DOUBLE_IS_32BITS) */ |
---|