1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtl.c */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
* Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* The argument reduction and testing for exceptional cases was
* written by Steven G. Kargl with input from Bruce D. Evans
* and David A. Schultz.
*/
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double cbrtl(long double x)
{
return cbrt(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#define BIAS (LDBL_MAX_EXP - 1)
static const unsigned
B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
long double cbrtl(long double x)
{
union IEEEl2bits u, v;
long double r, s, t, w;
double dr, dt, dx;
float ft, fx;
uint32_t hx;
uint16_t expsign;
int k;
u.e = x;
expsign = u.xbits.expsign;
k = expsign & 0x7fff;
/*
* If x = +-Inf, then cbrt(x) = +-Inf.
* If x = NaN, then cbrt(x) = NaN.
*/
if (k == BIAS + LDBL_MAX_EXP)
return x + x;
// FIXME: extended precision is default on linux..
#undef __i386__
#ifdef __i386__
fp_prec_t oprec;
oprec = fpgetprec();
if (oprec != FP_PE)
fpsetprec(FP_PE);
#endif
if (k == 0) {
/* If x = +-0, then cbrt(x) = +-0. */
if ((u.bits.manh | u.bits.manl) == 0) {
#ifdef __i386__
if (oprec != FP_PE)
fpsetprec(oprec);
#endif
return (x);
}
/* Adjust subnormal numbers. */
u.e *= 0x1.0p514;
k = u.bits.exp;
k -= BIAS + 514;
} else
k -= BIAS;
u.xbits.expsign = BIAS;
v.e = 1;
x = u.e;
switch (k % 3) {
case 1:
case -2:
x = 2*x;
k--;
break;
case 2:
case -1:
x = 4*x;
k -= 2;
break;
}
v.xbits.expsign = (expsign & 0x8000) | (BIAS + k / 3);
/*
* The following is the guts of s_cbrtf, with the handling of
* special values removed and extra care for accuracy not taken,
* but with most of the extra accuracy not discarded.
*/
/* ~5-bit estimate: */
fx = x;
GET_FLOAT_WORD(hx, fx);
SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
/* ~16-bit estimate: */
dx = x;
dt = ft;
dr = dt * dt * dt;
dt = dt * (dx + dx + dr) / (dx + dr + dr);
/* ~47-bit estimate: */
dr = dt * dt * dt;
dt = dt * (dx + dx + dr) / (dx + dr + dr);
#if LDBL_MANT_DIG == 64
/*
* dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
* Round it away from zero to 32 bits (32 so that t*t is exact, and
* away from zero for technical reasons).
*/
t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32;
#elif LDBL_MANT_DIG == 113
/*
* Round dt away from zero to 47 bits. Since we don't trust the 47,
* add 2 47-bit ulps instead of 1 to round up. Rounding is slow and
* might be avoidable in this case, since on most machines dt will
* have been evaluated in 53-bit precision and the technical reasons
* for rounding up might not apply to either case in cbrtl() since
* dt is much more accurate than needed.
*/
t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
#else
#error "Unsupported long double format"
#endif
/*
* Final step Newton iteration to 64 or 113 bits with
* error < 0.667 ulps
*/
s = t*t; /* t*t is exact */
r = x/s; /* error <= 0.5 ulps; |r| < |t| */
w = t+t; /* t+t is exact */
r = (r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
t *= v.e;
#ifdef __i386__
if (oprec != FP_PE)
fpsetprec(oprec);
#endif
return t;
}
#endif
|
