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
|
/* s_cbrtl.c -- long double version of s_cbrt.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* 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.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
#include "math.h"
#include "math_private.h"
/* cbrtl(x)
* Return cube root of x
*/
#ifdef __STDC__
static const u_int32_t
#else
static u_int32_t
#endif
B1_EXP = 10921, /* = Int(B1) */
B1_MANT = 0x7bc4b064, /* = Int(1.0-0.03306235651)*2**31 */
B2_EXP = 10900,
B2_MANT = 0x7bc4b064; /* = Int(1.0-0.03306235651)*2**31 */
#ifdef __STDC__
static const long double
#else
static long double
#endif
C = 5.42857142857142815906e-01L, /* 19/35 */
D = -7.05306122448979611050e-01L, /* -864/1225 */
E = 1.41428571428571436819e+00L, /* 99/70 */
F = 1.60714285714285720630e+00L, /* 45/28 */
G = 3.57142857142857150787e-01L; /* 5/14 */
#ifdef __STDC__
long double __cbrtl(long double x)
#else
long double __cbrtl(x)
long double x;
#endif
{
long double r,s,t=0.0,w;
u_int32_t sign, se, x0, x1;
GET_LDOUBLE_WORDS(se,x0,x1,x);
sign=se&0x8000; /* sign= sign(x) */
se ^= sign;
if(se==0x7fff) return(x+x); /* cbrt(NaN,INF) is itself */
if((se|x0|x1)==0)
return(x); /* cbrt(0) is itself */
SET_LDOUBLE_EXP(x,se); /* x <- |x| */
/* XXX I don't know whether the numbers for correct are correct. The
precalculation is extended from 20 bits to 32 bits. This hopefully
gives us the needed bits to get us still along with one iteration
step. */
/* rough cbrt to 5 bits */
if(se==0) /* subnormal number */
{
u_int64_t xxl;
u_int32_t set,t0,t1;
SET_LDOUBLE_EXP(t,0x4035); /* set t= 2**54 */
SET_LDOUBLE_MSW(t,0x80000000);
t*=x;
GET_LDOUBLE_WORDS(set,t0,t1,t);
xxl = ((u_int64_t) set) << 32 | t0;
xxl /= 3;
xxl += B2_EXP << 16 | B2_MANT;
t0 = xxl & 0xffffffffu;
set = xxl >> 32;
SET_LDOUBLE_WORDS(t,set,t0,t1);
}
else
{
u_int64_t xxl = ((u_int64_t) se) << 32 | x0;
xxl /= 3;
xxl += ((u_int64_t) B1_EXP) << 32 | B1_MANT;
SET_LDOUBLE_MSW(t,xxl&0xffffffffu);
xxl >>= 32;
SET_LDOUBLE_EXP(t,xxl);
}
/* new cbrt to 23 bits, may be implemented in single precision */
r=t*t/x;
s=C+r*t;
t*=G+F/(s+E+D/s);
/* chopped to 32 bits and make it larger than cbrt(x) */
GET_LDOUBLE_WORDS(se,x0,x1,t);
SET_LDOUBLE_WORDS(t,se,x0+1,0);
/* one step newton iteration to 53 bits with error less than 0.667 ulps */
s=t*t; /* t*t is exact */
r=x/s;
w=t+t;
r=(r-t)/(w+r); /* r-s is exact */
t=t+t*r;
/* retore the sign bit */
GET_LDOUBLE_EXP(se,t);
SET_LDOUBLE_EXP(t,se|sign);
return(t);
}
weak_alias (__cbrtl, cbrtl)
|
