OLD | NEW |
1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */ |
2 /* | 2 /* |
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> | 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
4 * | 4 * |
5 * Permission to use, copy, modify, and distribute this software for any | 5 * Permission to use, copy, modify, and distribute this software for any |
6 * purpose with or without fee is hereby granted, provided that the above | 6 * purpose with or without fee is hereby granted, provided that the above |
7 * copyright notice and this permission notice appear in all copies. | 7 * copyright notice and this permission notice appear in all copies. |
8 * | 8 * |
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES | 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF | 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
(...skipping 33 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
44 * IEEE -40,+40 10000 3.6e-19 7.9e-20 | 44 * IEEE -40,+40 10000 3.6e-19 7.9e-20 |
45 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 | 45 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 |
46 * | 46 * |
47 * Accuracy for large arguments is dominated by error in powl(). | 47 * Accuracy for large arguments is dominated by error in powl(). |
48 * | 48 * |
49 */ | 49 */ |
50 | 50 |
51 #include "libm.h" | 51 #include "libm.h" |
52 | 52 |
53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
54 long double tgammal(long double x) | 54 long double tgammal(long double x) { |
55 { | 55 return tgamma(x); |
56 » return tgamma(x); | |
57 } | 56 } |
58 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 | 57 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
59 /* | 58 /* |
60 tgamma(x+2) = tgamma(x+2) P(x)/Q(x) | 59 tgamma(x+2) = tgamma(x+2) P(x)/Q(x) |
61 0 <= x <= 1 | 60 0 <= x <= 1 |
62 Relative error | 61 Relative error |
63 n=7, d=8 | 62 n=7, d=8 |
64 Peak error = 1.83e-20 | 63 Peak error = 1.83e-20 |
65 Relative error spread = 8.4e-23 | 64 Relative error spread = 8.4e-23 |
66 */ | 65 */ |
67 static const long double P[8] = { | 66 static const long double P[8] = { |
68 4.212760487471622013093E-5L, | 67 4.212760487471622013093E-5L, 4.542931960608009155600E-4L, |
69 4.542931960608009155600E-4L, | 68 4.092666828394035500949E-3L, 2.385363243461108252554E-2L, |
70 4.092666828394035500949E-3L, | 69 1.113062816019361559013E-1L, 3.629515436640239168939E-1L, |
71 2.385363243461108252554E-2L, | 70 8.378004301573126728826E-1L, 1.000000000000000000009E0L, |
72 1.113062816019361559013E-1L, | |
73 3.629515436640239168939E-1L, | |
74 8.378004301573126728826E-1L, | |
75 1.000000000000000000009E0L, | |
76 }; | 71 }; |
77 static const long double Q[9] = { | 72 static const long double Q[9] = { |
78 -1.397148517476170440917E-5L, | 73 -1.397148517476170440917E-5L, 2.346584059160635244282E-4L, |
79 2.346584059160635244282E-4L, | 74 -1.237799246653152231188E-3L, -7.955933682494738320586E-4L, |
80 -1.237799246653152231188E-3L, | 75 2.773706565840072979165E-2L, -4.633887671244534213831E-2L, |
81 -7.955933682494738320586E-4L, | 76 -2.243510905670329164562E-1L, 4.150160950588455434583E-1L, |
82 2.773706565840072979165E-2L, | 77 9.999999999999999999908E-1L, |
83 -4.633887671244534213831E-2L, | |
84 -2.243510905670329164562E-1L, | |
85 4.150160950588455434583E-1L, | |
86 9.999999999999999999908E-1L, | |
87 }; | 78 }; |
88 | 79 |
89 /* | 80 /* |
90 static const long double P[] = { | 81 static const long double P[] = { |
91 -3.01525602666895735709e0L, | 82 -3.01525602666895735709e0L, |
92 -3.25157411956062339893e1L, | 83 -3.25157411956062339893e1L, |
93 -2.92929976820724030353e2L, | 84 -2.92929976820724030353e2L, |
94 -1.70730828800510297666e3L, | 85 -1.70730828800510297666e3L, |
95 -7.96667499622741999770e3L, | 86 -7.96667499622741999770e3L, |
96 -2.59780216007146401957e4L, | 87 -2.59780216007146401957e4L, |
(...skipping 18 matching lines...) Expand all Loading... |
115 /* Stirling's formula for the gamma function | 106 /* Stirling's formula for the gamma function |
116 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) | 107 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) |
117 z(x) = x | 108 z(x) = x |
118 13 <= x <= 1024 | 109 13 <= x <= 1024 |
119 Relative error | 110 Relative error |
120 n=8, d=0 | 111 n=8, d=0 |
121 Peak error = 9.44e-21 | 112 Peak error = 9.44e-21 |
122 Relative error spread = 8.8e-4 | 113 Relative error spread = 8.8e-4 |
123 */ | 114 */ |
124 static const long double STIR[9] = { | 115 static const long double STIR[9] = { |
125 7.147391378143610789273E-4L, | 116 7.147391378143610789273E-4L, -2.363848809501759061727E-5L, |
126 -2.363848809501759061727E-5L, | 117 -5.950237554056330156018E-4L, 6.989332260623193171870E-5L, |
127 -5.950237554056330156018E-4L, | 118 7.840334842744753003862E-4L, -2.294719747873185405699E-4L, |
128 6.989332260623193171870E-5L, | 119 -2.681327161876304418288E-3L, 3.472222222230075327854E-3L, |
129 7.840334842744753003862E-4L, | 120 8.333333333333331800504E-2L, |
130 -2.294719747873185405699E-4L, | |
131 -2.681327161876304418288E-3L, | |
132 3.472222222230075327854E-3L, | |
133 8.333333333333331800504E-2L, | |
134 }; | 121 }; |
135 | 122 |
136 #define MAXSTIR 1024.0L | 123 #define MAXSTIR 1024.0L |
137 static const long double SQTPI = 2.50662827463100050242E0L; | 124 static const long double SQTPI = 2.50662827463100050242E0L; |
138 | 125 |
139 /* 1/tgamma(x) = z P(z) | 126 /* 1/tgamma(x) = z P(z) |
140 * z(x) = 1/x | 127 * z(x) = 1/x |
141 * 0 < x < 0.03125 | 128 * 0 < x < 0.03125 |
142 * Peak relative error 4.2e-23 | 129 * Peak relative error 4.2e-23 |
143 */ | 130 */ |
144 static const long double S[9] = { | 131 static const long double S[9] = { |
145 -1.193945051381510095614E-3L, | 132 -1.193945051381510095614E-3L, 7.220599478036909672331E-3L, |
146 7.220599478036909672331E-3L, | 133 -9.622023360406271645744E-3L, -4.219773360705915470089E-2L, |
147 -9.622023360406271645744E-3L, | 134 1.665386113720805206758E-1L, -4.200263503403344054473E-2L, |
148 -4.219773360705915470089E-2L, | 135 -6.558780715202540684668E-1L, 5.772156649015328608253E-1L, |
149 1.665386113720805206758E-1L, | 136 1.000000000000000000000E0L, |
150 -4.200263503403344054473E-2L, | |
151 -6.558780715202540684668E-1L, | |
152 5.772156649015328608253E-1L, | |
153 1.000000000000000000000E0L, | |
154 }; | 137 }; |
155 | 138 |
156 /* 1/tgamma(-x) = z P(z) | 139 /* 1/tgamma(-x) = z P(z) |
157 * z(x) = 1/x | 140 * z(x) = 1/x |
158 * 0 < x < 0.03125 | 141 * 0 < x < 0.03125 |
159 * Peak relative error 5.16e-23 | 142 * Peak relative error 5.16e-23 |
160 * Relative error spread = 2.5e-24 | 143 * Relative error spread = 2.5e-24 |
161 */ | 144 */ |
162 static const long double SN[9] = { | 145 static const long double SN[9] = { |
163 1.133374167243894382010E-3L, | 146 1.133374167243894382010E-3L, 7.220837261893170325704E-3L, |
164 7.220837261893170325704E-3L, | 147 9.621911155035976733706E-3L, -4.219773343731191721664E-2L, |
165 9.621911155035976733706E-3L, | 148 -1.665386113944413519335E-1L, -4.200263503402112910504E-2L, |
166 -4.219773343731191721664E-2L, | 149 6.558780715202536547116E-1L, 5.772156649015328608727E-1L, |
167 -1.665386113944413519335E-1L, | 150 -1.000000000000000000000E0L, |
168 -4.200263503402112910504E-2L, | |
169 6.558780715202536547116E-1L, | |
170 5.772156649015328608727E-1L, | |
171 -1.000000000000000000000E0L, | |
172 }; | 151 }; |
173 | 152 |
174 static const long double PIL = 3.1415926535897932384626L; | 153 static const long double PIL = 3.1415926535897932384626L; |
175 | 154 |
176 /* Gamma function computed by Stirling's formula. | 155 /* Gamma function computed by Stirling's formula. |
177 */ | 156 */ |
178 static long double stirf(long double x) | 157 static long double stirf(long double x) { |
179 { | 158 long double y, w, v; |
180 » long double y, w, v; | |
181 | 159 |
182 » w = 1.0/x; | 160 w = 1.0 / x; |
183 » /* For large x, use rational coefficients from the analytical expansion.
*/ | 161 /* For large x, use rational coefficients from the analytical expansion. */ |
184 » if (x > 1024.0) | 162 if (x > 1024.0) |
185 » » w = (((((6.97281375836585777429E-5L * w | 163 w = (((((6.97281375836585777429E-5L * w + 7.84039221720066627474E-4L) * w - |
186 » » + 7.84039221720066627474E-4L) * w | 164 2.29472093621399176955E-4L) * |
187 » » - 2.29472093621399176955E-4L) * w | 165 w - |
188 » » - 2.68132716049382716049E-3L) * w | 166 2.68132716049382716049E-3L) * |
189 » » + 3.47222222222222222222E-3L) * w | 167 w + |
190 » » + 8.33333333333333333333E-2L) * w | 168 3.47222222222222222222E-3L) * |
191 » » + 1.0; | 169 w + |
192 » else | 170 8.33333333333333333333E-2L) * |
193 » » w = 1.0 + w * __polevll(w, STIR, 8); | 171 w + |
194 » y = expl(x); | 172 1.0; |
195 » if (x > MAXSTIR) { /* Avoid overflow in pow() */ | 173 else |
196 » » v = powl(x, 0.5L * x - 0.25L); | 174 w = 1.0 + w * __polevll(w, STIR, 8); |
197 » » y = v * (v / y); | 175 y = expl(x); |
198 » } else { | 176 if (x > MAXSTIR) { /* Avoid overflow in pow() */ |
199 » » y = powl(x, x - 0.5L) / y; | 177 v = powl(x, 0.5L * x - 0.25L); |
200 » } | 178 y = v * (v / y); |
201 » y = SQTPI * y * w; | 179 } else { |
202 » return y; | 180 y = powl(x, x - 0.5L) / y; |
| 181 } |
| 182 y = SQTPI * y * w; |
| 183 return y; |
203 } | 184 } |
204 | 185 |
205 long double tgammal(long double x) | 186 long double tgammal(long double x) { |
206 { | 187 long double p, q, z; |
207 » long double p, q, z; | |
208 | 188 |
209 » if (!isfinite(x)) | 189 if (!isfinite(x)) |
210 » » return x + INFINITY; | 190 return x + INFINITY; |
211 | 191 |
212 » q = fabsl(x); | 192 q = fabsl(x); |
213 » if (q > 13.0) { | 193 if (q > 13.0) { |
214 » » if (x < 0.0) { | 194 if (x < 0.0) { |
215 » » » p = floorl(q); | 195 p = floorl(q); |
216 » » » z = q - p; | 196 z = q - p; |
217 » » » if (z == 0) | 197 if (z == 0) |
218 » » » » return 0 / z; | 198 return 0 / z; |
219 » » » if (q > MAXGAML) { | 199 if (q > MAXGAML) { |
220 » » » » z = 0; | 200 z = 0; |
221 » » » } else { | 201 } else { |
222 » » » » if (z > 0.5) { | 202 if (z > 0.5) { |
223 » » » » » p += 1.0; | 203 p += 1.0; |
224 » » » » » z = q - p; | 204 z = q - p; |
225 » » » » } | 205 } |
226 » » » » z = q * sinl(PIL * z); | 206 z = q * sinl(PIL * z); |
227 » » » » z = fabsl(z) * stirf(q); | 207 z = fabsl(z) * stirf(q); |
228 » » » » z = PIL/z; | 208 z = PIL / z; |
229 » » » } | 209 } |
230 » » » if (0.5 * p == floorl(q * 0.5)) | 210 if (0.5 * p == floorl(q * 0.5)) |
231 » » » » z = -z; | 211 z = -z; |
232 » » } else if (x > MAXGAML) { | 212 } else if (x > MAXGAML) { |
233 » » » z = x * 0x1p16383L; | 213 z = x * 0x1p16383L; |
234 » » } else { | 214 } else { |
235 » » » z = stirf(x); | 215 z = stirf(x); |
236 » » } | 216 } |
237 » » return z; | 217 return z; |
238 » } | 218 } |
239 | 219 |
240 » z = 1.0; | 220 z = 1.0; |
241 » while (x >= 3.0) { | 221 while (x >= 3.0) { |
242 » » x -= 1.0; | 222 x -= 1.0; |
243 » » z *= x; | 223 z *= x; |
244 » } | 224 } |
245 » while (x < -0.03125L) { | 225 while (x < -0.03125L) { |
246 » » z /= x; | 226 z /= x; |
247 » » x += 1.0; | 227 x += 1.0; |
248 » } | 228 } |
249 » if (x <= 0.03125L) | 229 if (x <= 0.03125L) |
250 » » goto small; | 230 goto small; |
251 » while (x < 2.0) { | 231 while (x < 2.0) { |
252 » » z /= x; | 232 z /= x; |
253 » » x += 1.0; | 233 x += 1.0; |
254 » } | 234 } |
255 » if (x == 2.0) | 235 if (x == 2.0) |
256 » » return z; | 236 return z; |
257 | 237 |
258 » x -= 2.0; | 238 x -= 2.0; |
259 » p = __polevll(x, P, 7); | 239 p = __polevll(x, P, 7); |
260 » q = __polevll(x, Q, 8); | 240 q = __polevll(x, Q, 8); |
261 » z = z * p / q; | 241 z = z * p / q; |
262 » return z; | 242 return z; |
263 | 243 |
264 small: | 244 small: |
265 » /* z==1 if x was originally +-0 */ | 245 /* z==1 if x was originally +-0 */ |
266 » if (x == 0 && z != 1) | 246 if (x == 0 && z != 1) |
267 » » return x / x; | 247 return x / x; |
268 » if (x < 0.0) { | 248 if (x < 0.0) { |
269 » » x = -x; | 249 x = -x; |
270 » » q = z / (x * __polevll(x, SN, 8)); | 250 q = z / (x * __polevll(x, SN, 8)); |
271 » } else | 251 } else |
272 » » q = z / (x * __polevll(x, S, 8)); | 252 q = z / (x * __polevll(x, S, 8)); |
273 » return q; | 253 return q; |
274 } | 254 } |
275 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 | 255 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
276 // TODO: broken implementation to make things compile | 256 // TODO: broken implementation to make things compile |
277 long double tgammal(long double x) | 257 long double tgammal(long double x) { |
278 { | 258 return tgamma(x); |
279 » return tgamma(x); | |
280 } | 259 } |
281 #endif | 260 #endif |
OLD | NEW |