OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_lgammaf_r.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_lgammaf_r.c */ |
2 /* | 2 /* |
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 */ | 4 */ |
5 /* | 5 /* |
6 * ==================================================== | 6 * ==================================================== |
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 * | 8 * |
9 * Developed at SunPro, a Sun Microsystems, Inc. business. | 9 * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 * Permission to use, copy, modify, and distribute this | 10 * Permission to use, copy, modify, and distribute this |
11 * software is freely granted, provided that this notice | 11 * software is freely granted, provided that this notice |
12 * is preserved. | 12 * is preserved. |
13 * ==================================================== | 13 * ==================================================== |
14 */ | 14 */ |
15 | 15 |
16 #include "libm.h" | 16 #include "libm.h" |
17 #include "libc.h" | 17 #include "libc.h" |
18 | 18 |
19 static const float | 19 static const float pi = 3.1415927410e+00, /* 0x40490fdb */ |
20 pi = 3.1415927410e+00, /* 0x40490fdb */ | 20 a0 = 7.7215664089e-02, /* 0x3d9e233f */ |
21 a0 = 7.7215664089e-02, /* 0x3d9e233f */ | 21 a1 = 3.2246702909e-01, /* 0x3ea51a66 */ |
22 a1 = 3.2246702909e-01, /* 0x3ea51a66 */ | 22 a2 = 6.7352302372e-02, /* 0x3d89f001 */ |
23 a2 = 6.7352302372e-02, /* 0x3d89f001 */ | 23 a3 = 2.0580807701e-02, /* 0x3ca89915 */ |
24 a3 = 2.0580807701e-02, /* 0x3ca89915 */ | 24 a4 = 7.3855509982e-03, /* 0x3bf2027e */ |
25 a4 = 7.3855509982e-03, /* 0x3bf2027e */ | 25 a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ |
26 a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ | 26 a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ |
27 a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ | 27 a7 = 5.1006977446e-04, /* 0x3a05b634 */ |
28 a7 = 5.1006977446e-04, /* 0x3a05b634 */ | 28 a8 = 2.2086278477e-04, /* 0x39679767 */ |
29 a8 = 2.2086278477e-04, /* 0x39679767 */ | 29 a9 = 1.0801156895e-04, /* 0x38e28445 */ |
30 a9 = 1.0801156895e-04, /* 0x38e28445 */ | 30 a10 = 2.5214456400e-05, /* 0x37d383a2 */ |
31 a10 = 2.5214456400e-05, /* 0x37d383a2 */ | 31 a11 = 4.4864096708e-05, /* 0x383c2c75 */ |
32 a11 = 4.4864096708e-05, /* 0x383c2c75 */ | 32 tc = 1.4616321325e+00, /* 0x3fbb16c3 */ |
33 tc = 1.4616321325e+00, /* 0x3fbb16c3 */ | 33 tf = -1.2148628384e-01, /* 0xbdf8cdcd */ |
34 tf = -1.2148628384e-01, /* 0xbdf8cdcd */ | 34 /* tt = -(tail of tf) */ |
35 /* tt = -(tail of tf) */ | 35 tt = 6.6971006518e-09, /* 0x31e61c52 */ |
36 tt = 6.6971006518e-09, /* 0x31e61c52 */ | 36 t0 = 4.8383611441e-01, /* 0x3ef7b95e */ |
37 t0 = 4.8383611441e-01, /* 0x3ef7b95e */ | 37 t1 = -1.4758771658e-01, /* 0xbe17213c */ |
38 t1 = -1.4758771658e-01, /* 0xbe17213c */ | 38 t2 = 6.4624942839e-02, /* 0x3d845a15 */ |
39 t2 = 6.4624942839e-02, /* 0x3d845a15 */ | 39 t3 = -3.2788541168e-02, /* 0xbd064d47 */ |
40 t3 = -3.2788541168e-02, /* 0xbd064d47 */ | 40 t4 = 1.7970675603e-02, /* 0x3c93373d */ |
41 t4 = 1.7970675603e-02, /* 0x3c93373d */ | 41 t5 = -1.0314224288e-02, /* 0xbc28fcfe */ |
42 t5 = -1.0314224288e-02, /* 0xbc28fcfe */ | 42 t6 = 6.1005386524e-03, /* 0x3bc7e707 */ |
43 t6 = 6.1005386524e-03, /* 0x3bc7e707 */ | 43 t7 = -3.6845202558e-03, /* 0xbb7177fe */ |
44 t7 = -3.6845202558e-03, /* 0xbb7177fe */ | 44 t8 = 2.2596477065e-03, /* 0x3b141699 */ |
45 t8 = 2.2596477065e-03, /* 0x3b141699 */ | 45 t9 = -1.4034647029e-03, /* 0xbab7f476 */ |
46 t9 = -1.4034647029e-03, /* 0xbab7f476 */ | 46 t10 = 8.8108185446e-04, /* 0x3a66f867 */ |
47 t10 = 8.8108185446e-04, /* 0x3a66f867 */ | 47 t11 = -5.3859531181e-04, /* 0xba0d3085 */ |
48 t11 = -5.3859531181e-04, /* 0xba0d3085 */ | 48 t12 = 3.1563205994e-04, /* 0x39a57b6b */ |
49 t12 = 3.1563205994e-04, /* 0x39a57b6b */ | 49 t13 = -3.1275415677e-04, /* 0xb9a3f927 */ |
50 t13 = -3.1275415677e-04, /* 0xb9a3f927 */ | 50 t14 = 3.3552918467e-04, /* 0x39afe9f7 */ |
51 t14 = 3.3552918467e-04, /* 0x39afe9f7 */ | 51 u0 = -7.7215664089e-02, /* 0xbd9e233f */ |
52 u0 = -7.7215664089e-02, /* 0xbd9e233f */ | 52 u1 = 6.3282704353e-01, /* 0x3f2200f4 */ |
53 u1 = 6.3282704353e-01, /* 0x3f2200f4 */ | 53 u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ |
54 u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ | 54 u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ |
55 u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ | 55 u4 = 2.2896373272e-01, /* 0x3e6a7578 */ |
56 u4 = 2.2896373272e-01, /* 0x3e6a7578 */ | 56 u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ |
57 u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ | 57 v1 = 2.4559779167e+00, /* 0x401d2ebe */ |
58 v1 = 2.4559779167e+00, /* 0x401d2ebe */ | 58 v2 = 2.1284897327e+00, /* 0x4008392d */ |
59 v2 = 2.1284897327e+00, /* 0x4008392d */ | 59 v3 = 7.6928514242e-01, /* 0x3f44efdf */ |
60 v3 = 7.6928514242e-01, /* 0x3f44efdf */ | 60 v4 = 1.0422264785e-01, /* 0x3dd572af */ |
61 v4 = 1.0422264785e-01, /* 0x3dd572af */ | 61 v5 = 3.2170924824e-03, /* 0x3b52d5db */ |
62 v5 = 3.2170924824e-03, /* 0x3b52d5db */ | 62 s0 = -7.7215664089e-02, /* 0xbd9e233f */ |
63 s0 = -7.7215664089e-02, /* 0xbd9e233f */ | 63 s1 = 2.1498242021e-01, /* 0x3e5c245a */ |
64 s1 = 2.1498242021e-01, /* 0x3e5c245a */ | 64 s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ |
65 s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ | 65 s3 = 1.4635047317e-01, /* 0x3e15dce6 */ |
66 s3 = 1.4635047317e-01, /* 0x3e15dce6 */ | 66 s4 = 2.6642270386e-02, /* 0x3cda40e4 */ |
67 s4 = 2.6642270386e-02, /* 0x3cda40e4 */ | 67 s5 = 1.8402845599e-03, /* 0x3af135b4 */ |
68 s5 = 1.8402845599e-03, /* 0x3af135b4 */ | 68 s6 = 3.1947532989e-05, /* 0x3805ff67 */ |
69 s6 = 3.1947532989e-05, /* 0x3805ff67 */ | 69 r1 = 1.3920053244e+00, /* 0x3fb22d3b */ |
70 r1 = 1.3920053244e+00, /* 0x3fb22d3b */ | 70 r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ |
71 r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ | 71 r3 = 1.7193385959e-01, /* 0x3e300f6e */ |
72 r3 = 1.7193385959e-01, /* 0x3e300f6e */ | 72 r4 = 1.8645919859e-02, /* 0x3c98bf54 */ |
73 r4 = 1.8645919859e-02, /* 0x3c98bf54 */ | 73 r5 = 7.7794247773e-04, /* 0x3a4beed6 */ |
74 r5 = 7.7794247773e-04, /* 0x3a4beed6 */ | 74 r6 = 7.3266842264e-06, /* 0x36f5d7bd */ |
75 r6 = 7.3266842264e-06, /* 0x36f5d7bd */ | 75 w0 = 4.1893854737e-01, /* 0x3ed67f1d */ |
76 w0 = 4.1893854737e-01, /* 0x3ed67f1d */ | 76 w1 = 8.3333335817e-02, /* 0x3daaaaab */ |
77 w1 = 8.3333335817e-02, /* 0x3daaaaab */ | 77 w2 = -2.7777778450e-03, /* 0xbb360b61 */ |
78 w2 = -2.7777778450e-03, /* 0xbb360b61 */ | 78 w3 = 7.9365057172e-04, /* 0x3a500cfd */ |
79 w3 = 7.9365057172e-04, /* 0x3a500cfd */ | 79 w4 = -5.9518753551e-04, /* 0xba1c065c */ |
80 w4 = -5.9518753551e-04, /* 0xba1c065c */ | 80 w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ |
81 w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ | 81 w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ |
82 w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ | |
83 | 82 |
84 /* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ | 83 /* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ |
85 static float sin_pi(float x) | 84 static float sin_pi(float x) { |
86 { | 85 double_t y; |
87 » double_t y; | 86 int n; |
88 » int n; | 87 |
89 | 88 /* spurious inexact if odd int */ |
90 » /* spurious inexact if odd int */ | 89 x = 2 * (x * 0.5f - floorf(x * 0.5f)); /* x mod 2.0 */ |
91 » x = 2*(x*0.5f - floorf(x*0.5f)); /* x mod 2.0 */ | 90 |
92 | 91 n = (int)(x * 4); |
93 » n = (int)(x*4); | 92 n = (n + 1) / 2; |
94 » n = (n+1)/2; | 93 y = x - n * 0.5f; |
95 » y = x - n*0.5f; | 94 y *= 3.14159265358979323846; |
96 » y *= 3.14159265358979323846; | 95 switch (n) { |
97 » switch (n) { | 96 default: /* case 4: */ |
98 » default: /* case 4: */ | 97 case 0: |
99 » case 0: return __sindf(y); | 98 return __sindf(y); |
100 » case 1: return __cosdf(y); | 99 case 1: |
101 » case 2: return __sindf(-y); | 100 return __cosdf(y); |
102 » case 3: return -__cosdf(y); | 101 case 2: |
103 » } | 102 return __sindf(-y); |
| 103 case 3: |
| 104 return -__cosdf(y); |
| 105 } |
104 } | 106 } |
105 | 107 |
106 float __lgammaf_r(float x, int *signgamp) | 108 float __lgammaf_r(float x, int* signgamp) { |
107 { | 109 union { |
108 » union {float f; uint32_t i;} u = {x}; | 110 float f; |
109 » float t,y,z,nadj,p,p1,p2,p3,q,r,w; | 111 uint32_t i; |
110 » uint32_t ix; | 112 } u = {x}; |
111 » int i,sign; | 113 float t, y, z, nadj, p, p1, p2, p3, q, r, w; |
112 | 114 uint32_t ix; |
113 » /* purge off +-inf, NaN, +-0, tiny and negative arguments */ | 115 int i, sign; |
114 » *signgamp = 1; | 116 |
115 » sign = u.i>>31; | 117 /* purge off +-inf, NaN, +-0, tiny and negative arguments */ |
116 » ix = u.i & 0x7fffffff; | 118 *signgamp = 1; |
117 » if (ix >= 0x7f800000) | 119 sign = u.i >> 31; |
118 » » return x*x; | 120 ix = u.i & 0x7fffffff; |
119 » if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ | 121 if (ix >= 0x7f800000) |
120 » » if (sign) { | 122 return x * x; |
121 » » » *signgamp = -1; | 123 if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ |
122 » » » x = -x; | 124 if (sign) { |
123 » » } | 125 *signgamp = -1; |
124 » » return -logf(x); | 126 x = -x; |
125 » } | 127 } |
126 » if (sign) { | 128 return -logf(x); |
127 » » x = -x; | 129 } |
128 » » t = sin_pi(x); | 130 if (sign) { |
129 » » if (t == 0.0f) /* -integer */ | 131 x = -x; |
130 » » » return 1.0f/(x-x); | 132 t = sin_pi(x); |
131 » » if (t > 0.0f) | 133 if (t == 0.0f) /* -integer */ |
132 » » » *signgamp = -1; | 134 return 1.0f / (x - x); |
133 » » else | 135 if (t > 0.0f) |
134 » » » t = -t; | 136 *signgamp = -1; |
135 » » nadj = logf(pi/(t*x)); | 137 else |
136 » } | 138 t = -t; |
137 | 139 nadj = logf(pi / (t * x)); |
138 » /* purge off 1 and 2 */ | 140 } |
139 » if (ix == 0x3f800000 || ix == 0x40000000) | 141 |
140 » » r = 0; | 142 /* purge off 1 and 2 */ |
141 » /* for x < 2.0 */ | 143 if (ix == 0x3f800000 || ix == 0x40000000) |
142 » else if (ix < 0x40000000) { | 144 r = 0; |
143 » » if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ | 145 /* for x < 2.0 */ |
144 » » » r = -logf(x); | 146 else if (ix < 0x40000000) { |
145 » » » if (ix >= 0x3f3b4a20) { | 147 if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ |
146 » » » » y = 1.0f - x; | 148 r = -logf(x); |
147 » » » » i = 0; | 149 if (ix >= 0x3f3b4a20) { |
148 » » » } else if (ix >= 0x3e6d3308) { | 150 y = 1.0f - x; |
149 » » » » y = x - (tc-1.0f); | 151 i = 0; |
150 » » » » i = 1; | 152 } else if (ix >= 0x3e6d3308) { |
151 » » » } else { | 153 y = x - (tc - 1.0f); |
152 » » » » y = x; | 154 i = 1; |
153 » » » » i = 2; | 155 } else { |
154 » » » } | 156 y = x; |
155 » » } else { | 157 i = 2; |
156 » » » r = 0.0f; | 158 } |
157 » » » if (ix >= 0x3fdda618) { /* [1.7316,2] */ | 159 } else { |
158 » » » » y = 2.0f - x; | 160 r = 0.0f; |
159 » » » » i = 0; | 161 if (ix >= 0x3fdda618) { /* [1.7316,2] */ |
160 » » » } else if (ix >= 0x3F9da620) { /* [1.23,1.73] */ | 162 y = 2.0f - x; |
161 » » » » y = x - tc; | 163 i = 0; |
162 » » » » i = 1; | 164 } else if (ix >= 0x3F9da620) { /* [1.23,1.73] */ |
163 » » » } else { | 165 y = x - tc; |
164 » » » » y = x - 1.0f; | 166 i = 1; |
165 » » » » i = 2; | 167 } else { |
166 » » » } | 168 y = x - 1.0f; |
167 » » } | 169 i = 2; |
168 » » switch(i) { | 170 } |
169 » » case 0: | 171 } |
170 » » » z = y*y; | 172 switch (i) { |
171 » » » p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); | 173 case 0: |
172 » » » p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); | 174 z = y * y; |
173 » » » p = y*p1+p2; | 175 p1 = a0 + z * (a2 + z * (a4 + z * (a6 + z * (a8 + z * a10)))); |
174 » » » r += p - 0.5f*y; | 176 p2 = z * (a1 + z * (a3 + z * (a5 + z * (a7 + z * (a9 + z * a11))))); |
175 » » » break; | 177 p = y * p1 + p2; |
176 » » case 1: | 178 r += p - 0.5f * y; |
177 » » » z = y*y; | 179 break; |
178 » » » w = z*y; | 180 case 1: |
179 » » » p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp
*/ | 181 z = y * y; |
180 » » » p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); | 182 w = z * y; |
181 » » » p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); | 183 p1 = t0 + w * (t3 + w * (t6 + w * (t9 + w * t12))); /* parallel comp */ |
182 » » » p = z*p1-(tt-w*(p2+y*p3)); | 184 p2 = t1 + w * (t4 + w * (t7 + w * (t10 + w * t13))); |
183 » » » r += (tf + p); | 185 p3 = t2 + w * (t5 + w * (t8 + w * (t11 + w * t14))); |
184 » » » break; | 186 p = z * p1 - (tt - w * (p2 + y * p3)); |
185 » » case 2: | 187 r += (tf + p); |
186 » » » p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); | 188 break; |
187 » » » p2 = 1.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); | 189 case 2: |
188 » » » r += -0.5f*y + p1/p2; | 190 p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * u5))))); |
189 » » } | 191 p2 = 1.0f + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5)))); |
190 » } else if (ix < 0x41000000) { /* x < 8.0 */ | 192 r += -0.5f * y + p1 / p2; |
191 » » i = (int)x; | 193 } |
192 » » y = x - (float)i; | 194 } else if (ix < 0x41000000) { /* x < 8.0 */ |
193 » » p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); | 195 i = (int)x; |
194 » » q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); | 196 y = x - (float)i; |
195 » » r = 0.5f*y+p/q; | 197 p = y * |
196 » » z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */ | 198 (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))); |
197 » » switch (i) { | 199 q = 1.0f + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))))); |
198 » » case 7: z *= y + 6.0f; /* FALLTHRU */ | 200 r = 0.5f * y + p / q; |
199 » » case 6: z *= y + 5.0f; /* FALLTHRU */ | 201 z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */ |
200 » » case 5: z *= y + 4.0f; /* FALLTHRU */ | 202 switch (i) { |
201 » » case 4: z *= y + 3.0f; /* FALLTHRU */ | 203 case 7: |
202 » » case 3: z *= y + 2.0f; /* FALLTHRU */ | 204 z *= y + 6.0f; /* FALLTHRU */ |
203 » » » r += logf(z); | 205 case 6: |
204 » » » break; | 206 z *= y + 5.0f; /* FALLTHRU */ |
205 » » } | 207 case 5: |
206 » } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */ | 208 z *= y + 4.0f; /* FALLTHRU */ |
207 » » t = logf(x); | 209 case 4: |
208 » » z = 1.0f/x; | 210 z *= y + 3.0f; /* FALLTHRU */ |
209 » » y = z*z; | 211 case 3: |
210 » » w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); | 212 z *= y + 2.0f; /* FALLTHRU */ |
211 » » r = (x-0.5f)*(t-1.0f)+w; | 213 r += logf(z); |
212 » } else /* 2**58 <= x <= inf */ | 214 break; |
213 » » r = x*(logf(x)-1.0f); | 215 } |
214 » if (sign) | 216 } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */ |
215 » » r = nadj - r; | 217 t = logf(x); |
216 » return r; | 218 z = 1.0f / x; |
| 219 y = z * z; |
| 220 w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6))))); |
| 221 r = (x - 0.5f) * (t - 1.0f) + w; |
| 222 } else /* 2**58 <= x <= inf */ |
| 223 r = x * (logf(x) - 1.0f); |
| 224 if (sign) |
| 225 r = nadj - r; |
| 226 return r; |
217 } | 227 } |
218 | 228 |
219 weak_alias(__lgammaf_r, lgammaf_r); | 229 weak_alias(__lgammaf_r, lgammaf_r); |
OLD | NEW |