1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtl.c */
|
---|
2 | /*-
|
---|
3 | * ====================================================
|
---|
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
---|
5 | * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
|
---|
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 | * The argument reduction and testing for exceptional cases was
|
---|
14 | * written by Steven G. Kargl with input from Bruce D. Evans
|
---|
15 | * and David A. Schultz.
|
---|
16 | */
|
---|
17 |
|
---|
18 | #include "libm.h"
|
---|
19 |
|
---|
20 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
---|
21 | long double cbrtl(long double x)
|
---|
22 | {
|
---|
23 | return cbrt(x);
|
---|
24 | }
|
---|
25 | #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
|
---|
26 | static const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
|
---|
27 |
|
---|
28 | long double cbrtl(long double x)
|
---|
29 | {
|
---|
30 | union ldshape u = {x}, v;
|
---|
31 | union {float f; uint32_t i;} uft;
|
---|
32 | long double r, s, t, w;
|
---|
33 | double_t dr, dt, dx;
|
---|
34 | float_t ft;
|
---|
35 | int e = u.i.se & 0x7fff;
|
---|
36 | int sign = u.i.se & 0x8000;
|
---|
37 |
|
---|
38 | /*
|
---|
39 | * If x = +-Inf, then cbrt(x) = +-Inf.
|
---|
40 | * If x = NaN, then cbrt(x) = NaN.
|
---|
41 | */
|
---|
42 | if (e == 0x7fff)
|
---|
43 | return x + x;
|
---|
44 | if (e == 0) {
|
---|
45 | /* Adjust subnormal numbers. */
|
---|
46 | u.f *= 0x1p120;
|
---|
47 | e = u.i.se & 0x7fff;
|
---|
48 | /* If x = +-0, then cbrt(x) = +-0. */
|
---|
49 | if (e == 0)
|
---|
50 | return x;
|
---|
51 | e -= 120;
|
---|
52 | }
|
---|
53 | e -= 0x3fff;
|
---|
54 | u.i.se = 0x3fff;
|
---|
55 | x = u.f;
|
---|
56 | switch (e % 3) {
|
---|
57 | case 1:
|
---|
58 | case -2:
|
---|
59 | x *= 2;
|
---|
60 | e--;
|
---|
61 | break;
|
---|
62 | case 2:
|
---|
63 | case -1:
|
---|
64 | x *= 4;
|
---|
65 | e -= 2;
|
---|
66 | break;
|
---|
67 | }
|
---|
68 | v.f = 1.0;
|
---|
69 | v.i.se = sign | (0x3fff + e/3);
|
---|
70 |
|
---|
71 | /*
|
---|
72 | * The following is the guts of s_cbrtf, with the handling of
|
---|
73 | * special values removed and extra care for accuracy not taken,
|
---|
74 | * but with most of the extra accuracy not discarded.
|
---|
75 | */
|
---|
76 |
|
---|
77 | /* ~5-bit estimate: */
|
---|
78 | uft.f = x;
|
---|
79 | uft.i = (uft.i & 0x7fffffff)/3 + B1;
|
---|
80 | ft = uft.f;
|
---|
81 |
|
---|
82 | /* ~16-bit estimate: */
|
---|
83 | dx = x;
|
---|
84 | dt = ft;
|
---|
85 | dr = dt * dt * dt;
|
---|
86 | dt = dt * (dx + dx + dr) / (dx + dr + dr);
|
---|
87 |
|
---|
88 | /* ~47-bit estimate: */
|
---|
89 | dr = dt * dt * dt;
|
---|
90 | dt = dt * (dx + dx + dr) / (dx + dr + dr);
|
---|
91 |
|
---|
92 | #if LDBL_MANT_DIG == 64
|
---|
93 | /*
|
---|
94 | * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
|
---|
95 | * Round it away from zero to 32 bits (32 so that t*t is exact, and
|
---|
96 | * away from zero for technical reasons).
|
---|
97 | */
|
---|
98 | t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32;
|
---|
99 | #elif LDBL_MANT_DIG == 113
|
---|
100 | /*
|
---|
101 | * Round dt away from zero to 47 bits. Since we don't trust the 47,
|
---|
102 | * add 2 47-bit ulps instead of 1 to round up. Rounding is slow and
|
---|
103 | * might be avoidable in this case, since on most machines dt will
|
---|
104 | * have been evaluated in 53-bit precision and the technical reasons
|
---|
105 | * for rounding up might not apply to either case in cbrtl() since
|
---|
106 | * dt is much more accurate than needed.
|
---|
107 | */
|
---|
108 | t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
|
---|
109 | #endif
|
---|
110 |
|
---|
111 | /*
|
---|
112 | * Final step Newton iteration to 64 or 113 bits with
|
---|
113 | * error < 0.667 ulps
|
---|
114 | */
|
---|
115 | s = t*t; /* t*t is exact */
|
---|
116 | r = x/s; /* error <= 0.5 ulps; |r| < |t| */
|
---|
117 | w = t+t; /* t+t is exact */
|
---|
118 | r = (r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
|
---|
119 | t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
|
---|
120 |
|
---|
121 | t *= v.f;
|
---|
122 | return t;
|
---|
123 | }
|
---|
124 | #endif
|
---|