]>
Commit | Line | Data |
---|---|---|
a3c937ce AJ |
1 | /* lgammal |
2 | * | |
3 | * Natural logarithm of gamma function | |
4 | * | |
5 | * | |
6 | * | |
7 | * SYNOPSIS: | |
8 | * | |
9 | * long double x, y, lgammal(); | |
10 | * extern int sgngam; | |
11 | * | |
12 | * y = lgammal(x); | |
13 | * | |
14 | * | |
15 | * | |
16 | * DESCRIPTION: | |
17 | * | |
18 | * Returns the base e (2.718...) logarithm of the absolute | |
19 | * value of the gamma function of the argument. | |
20 | * The sign (+1 or -1) of the gamma function is returned in a | |
21 | * global (extern) variable named sgngam. | |
22 | * | |
23 | * The positive domain is partitioned into numerous segments for approximation. | |
24 | * For x > 10, | |
25 | * log gamma(x) = (x - 0.5) log(x) - x + log sqrt(2 pi) + 1/x R(1/x^2) | |
26 | * Near the minimum at x = x0 = 1.46... the approximation is | |
27 | * log gamma(x0 + z) = log gamma(x0) + z^2 P(z)/Q(z) | |
28 | * for small z. | |
29 | * Elsewhere between 0 and 10, | |
30 | * log gamma(n + z) = log gamma(n) + z P(z)/Q(z) | |
31 | * for various selected n and small z. | |
32 | * | |
33 | * The cosecant reflection formula is employed for negative arguments. | |
34 | * | |
35 | * | |
36 | * | |
37 | * ACCURACY: | |
38 | * | |
39 | * | |
40 | * arithmetic domain # trials peak rms | |
41 | * Relative error: | |
42 | * IEEE 10, 30 100000 3.9e-34 9.8e-35 | |
43 | * IEEE 0, 10 100000 3.8e-34 5.3e-35 | |
44 | * Absolute error: | |
45 | * IEEE -10, 0 100000 8.0e-34 8.0e-35 | |
46 | * IEEE -30, -10 100000 4.4e-34 1.0e-34 | |
47 | * IEEE -100, 100 100000 1.0e-34 | |
48 | * | |
49 | * The absolute error criterion is the same as relative error | |
50 | * when the function magnitude is greater than one but it is absolute | |
51 | * when the magnitude is less than one. | |
52 | * | |
53 | */ | |
54 | ||
cc7375ce RM |
55 | /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov> |
56 | ||
57 | This library is free software; you can redistribute it and/or | |
58 | modify it under the terms of the GNU Lesser General Public | |
59 | License as published by the Free Software Foundation; either | |
60 | version 2.1 of the License, or (at your option) any later version. | |
61 | ||
62 | This library is distributed in the hope that it will be useful, | |
63 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
64 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
65 | Lesser General Public License for more details. | |
66 | ||
67 | You should have received a copy of the GNU Lesser General Public | |
59ba27a6 PE |
68 | License along with this library; if not, see |
69 | <http://www.gnu.org/licenses/>. */ | |
a3c937ce | 70 | |
1ed0291c RH |
71 | #include <math.h> |
72 | #include <math_private.h> | |
9d969099 | 73 | #include <float.h> |
a3c937ce | 74 | |
1f5649f8 | 75 | static const long double PIL = 3.1415926535897932384626433832795028841972E0L; |
4b84e247 JM |
76 | #if LDBL_MANT_DIG == 106 |
77 | static const long double MAXLGM = 0x5.d53649e2d469dbc1f01e99fd66p+1012L; | |
78 | #else | |
1f5649f8 | 79 | static const long double MAXLGM = 1.0485738685148938358098967157129705071571E4928L; |
4b84e247 | 80 | #endif |
1f5649f8 | 81 | static const long double one = 1.0L; |
9d969099 | 82 | static const long double huge = LDBL_MAX; |
a3c937ce AJ |
83 | |
84 | /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2) | |
85 | 1/x <= 0.0741 (x >= 13.495...) | |
86 | Peak relative error 1.5e-36 */ | |
1f5649f8 | 87 | static const long double ls2pi = 9.1893853320467274178032973640561763986140E-1L; |
a3c937ce | 88 | #define NRASY 12 |
1f5649f8 | 89 | static const long double RASY[NRASY + 1] = |
a3c937ce AJ |
90 | { |
91 | 8.333333333333333333333333333310437112111E-2L, | |
92 | -2.777777777777777777777774789556228296902E-3L, | |
93 | 7.936507936507936507795933938448586499183E-4L, | |
94 | -5.952380952380952041799269756378148574045E-4L, | |
95 | 8.417508417507928904209891117498524452523E-4L, | |
96 | -1.917526917481263997778542329739806086290E-3L, | |
97 | 6.410256381217852504446848671499409919280E-3L, | |
98 | -2.955064066900961649768101034477363301626E-2L, | |
99 | 1.796402955865634243663453415388336954675E-1L, | |
100 | -1.391522089007758553455753477688592767741E0L, | |
101 | 1.326130089598399157988112385013829305510E1L, | |
102 | -1.420412699593782497803472576479997819149E2L, | |
103 | 1.218058922427762808938869872528846787020E3L | |
104 | }; | |
105 | ||
106 | ||
107 | /* log gamma(x+13) = log gamma(13) + x P(x)/Q(x) | |
108 | -0.5 <= x <= 0.5 | |
109 | 12.5 <= x+13 <= 13.5 | |
110 | Peak relative error 1.1e-36 */ | |
1f5649f8 UD |
111 | static const long double lgam13a = 1.9987213134765625E1L; |
112 | static const long double lgam13b = 1.3608962611495173623870550785125024484248E-6L; | |
a3c937ce | 113 | #define NRN13 7 |
1f5649f8 | 114 | static const long double RN13[NRN13 + 1] = |
a3c937ce AJ |
115 | { |
116 | 8.591478354823578150238226576156275285700E11L, | |
117 | 2.347931159756482741018258864137297157668E11L, | |
118 | 2.555408396679352028680662433943000804616E10L, | |
119 | 1.408581709264464345480765758902967123937E9L, | |
120 | 4.126759849752613822953004114044451046321E7L, | |
121 | 6.133298899622688505854211579222889943778E5L, | |
122 | 3.929248056293651597987893340755876578072E3L, | |
123 | 6.850783280018706668924952057996075215223E0L | |
124 | }; | |
125 | #define NRD13 6 | |
1f5649f8 | 126 | static const long double RD13[NRD13 + 1] = |
a3c937ce AJ |
127 | { |
128 | 3.401225382297342302296607039352935541669E11L, | |
129 | 8.756765276918037910363513243563234551784E10L, | |
130 | 8.873913342866613213078554180987647243903E9L, | |
131 | 4.483797255342763263361893016049310017973E8L, | |
132 | 1.178186288833066430952276702931512870676E7L, | |
133 | 1.519928623743264797939103740132278337476E5L, | |
134 | 7.989298844938119228411117593338850892311E2L | |
135 | /* 1.0E0L */ | |
136 | }; | |
137 | ||
138 | ||
139 | /* log gamma(x+12) = log gamma(12) + x P(x)/Q(x) | |
140 | -0.5 <= x <= 0.5 | |
141 | 11.5 <= x+12 <= 12.5 | |
142 | Peak relative error 4.1e-36 */ | |
1f5649f8 UD |
143 | static const long double lgam12a = 1.75023040771484375E1L; |
144 | static const long double lgam12b = 3.7687254483392876529072161996717039575982E-6L; | |
a3c937ce | 145 | #define NRN12 7 |
1f5649f8 | 146 | static const long double RN12[NRN12 + 1] = |
a3c937ce AJ |
147 | { |
148 | 4.709859662695606986110997348630997559137E11L, | |
149 | 1.398713878079497115037857470168777995230E11L, | |
150 | 1.654654931821564315970930093932954900867E10L, | |
151 | 9.916279414876676861193649489207282144036E8L, | |
152 | 3.159604070526036074112008954113411389879E7L, | |
153 | 5.109099197547205212294747623977502492861E5L, | |
154 | 3.563054878276102790183396740969279826988E3L, | |
155 | 6.769610657004672719224614163196946862747E0L | |
156 | }; | |
157 | #define NRD12 6 | |
1f5649f8 | 158 | static const long double RD12[NRD12 + 1] = |
a3c937ce AJ |
159 | { |
160 | 1.928167007860968063912467318985802726613E11L, | |
161 | 5.383198282277806237247492369072266389233E10L, | |
162 | 5.915693215338294477444809323037871058363E9L, | |
163 | 3.241438287570196713148310560147925781342E8L, | |
164 | 9.236680081763754597872713592701048455890E6L, | |
165 | 1.292246897881650919242713651166596478850E5L, | |
166 | 7.366532445427159272584194816076600211171E2L | |
167 | /* 1.0E0L */ | |
168 | }; | |
169 | ||
170 | ||
171 | /* log gamma(x+11) = log gamma(11) + x P(x)/Q(x) | |
172 | -0.5 <= x <= 0.5 | |
173 | 10.5 <= x+11 <= 11.5 | |
174 | Peak relative error 1.8e-35 */ | |
1f5649f8 UD |
175 | static const long double lgam11a = 1.5104400634765625E1L; |
176 | static const long double lgam11b = 1.1938309890295225709329251070371882250744E-5L; | |
a3c937ce | 177 | #define NRN11 7 |
1f5649f8 | 178 | static const long double RN11[NRN11 + 1] = |
a3c937ce AJ |
179 | { |
180 | 2.446960438029415837384622675816736622795E11L, | |
181 | 7.955444974446413315803799763901729640350E10L, | |
182 | 1.030555327949159293591618473447420338444E10L, | |
183 | 6.765022131195302709153994345470493334946E8L, | |
184 | 2.361892792609204855279723576041468347494E7L, | |
185 | 4.186623629779479136428005806072176490125E5L, | |
186 | 3.202506022088912768601325534149383594049E3L, | |
187 | 6.681356101133728289358838690666225691363E0L | |
188 | }; | |
189 | #define NRD11 6 | |
1f5649f8 | 190 | static const long double RD11[NRD11 + 1] = |
a3c937ce AJ |
191 | { |
192 | 1.040483786179428590683912396379079477432E11L, | |
193 | 3.172251138489229497223696648369823779729E10L, | |
194 | 3.806961885984850433709295832245848084614E9L, | |
195 | 2.278070344022934913730015420611609620171E8L, | |
196 | 7.089478198662651683977290023829391596481E6L, | |
197 | 1.083246385105903533237139380509590158658E5L, | |
198 | 6.744420991491385145885727942219463243597E2L | |
199 | /* 1.0E0L */ | |
200 | }; | |
201 | ||
202 | ||
203 | /* log gamma(x+10) = log gamma(10) + x P(x)/Q(x) | |
204 | -0.5 <= x <= 0.5 | |
205 | 9.5 <= x+10 <= 10.5 | |
206 | Peak relative error 5.4e-37 */ | |
1f5649f8 UD |
207 | static const long double lgam10a = 1.280181884765625E1L; |
208 | static const long double lgam10b = 8.6324252196112077178745667061642811492557E-6L; | |
a3c937ce | 209 | #define NRN10 7 |
1f5649f8 | 210 | static const long double RN10[NRN10 + 1] = |
a3c937ce AJ |
211 | { |
212 | -1.239059737177249934158597996648808363783E14L, | |
213 | -4.725899566371458992365624673357356908719E13L, | |
214 | -7.283906268647083312042059082837754850808E12L, | |
215 | -5.802855515464011422171165179767478794637E11L, | |
216 | -2.532349691157548788382820303182745897298E10L, | |
217 | -5.884260178023777312587193693477072061820E8L, | |
218 | -6.437774864512125749845840472131829114906E6L, | |
219 | -2.350975266781548931856017239843273049384E4L | |
220 | }; | |
221 | #define NRD10 7 | |
1f5649f8 | 222 | static const long double RD10[NRD10 + 1] = |
a3c937ce AJ |
223 | { |
224 | -5.502645997581822567468347817182347679552E13L, | |
225 | -1.970266640239849804162284805400136473801E13L, | |
226 | -2.819677689615038489384974042561531409392E12L, | |
227 | -2.056105863694742752589691183194061265094E11L, | |
228 | -8.053670086493258693186307810815819662078E9L, | |
229 | -1.632090155573373286153427982504851867131E8L, | |
230 | -1.483575879240631280658077826889223634921E6L, | |
231 | -4.002806669713232271615885826373550502510E3L | |
232 | /* 1.0E0L */ | |
233 | }; | |
234 | ||
235 | ||
236 | /* log gamma(x+9) = log gamma(9) + x P(x)/Q(x) | |
237 | -0.5 <= x <= 0.5 | |
238 | 8.5 <= x+9 <= 9.5 | |
239 | Peak relative error 3.6e-36 */ | |
1f5649f8 UD |
240 | static const long double lgam9a = 1.06045989990234375E1L; |
241 | static const long double lgam9b = 3.9037218127284172274007216547549861681400E-6L; | |
a3c937ce | 242 | #define NRN9 7 |
1f5649f8 | 243 | static const long double RN9[NRN9 + 1] = |
a3c937ce AJ |
244 | { |
245 | -4.936332264202687973364500998984608306189E13L, | |
246 | -2.101372682623700967335206138517766274855E13L, | |
247 | -3.615893404644823888655732817505129444195E12L, | |
248 | -3.217104993800878891194322691860075472926E11L, | |
249 | -1.568465330337375725685439173603032921399E10L, | |
250 | -4.073317518162025744377629219101510217761E8L, | |
251 | -4.983232096406156139324846656819246974500E6L, | |
252 | -2.036280038903695980912289722995505277253E4L | |
253 | }; | |
254 | #define NRD9 7 | |
1f5649f8 | 255 | static const long double RD9[NRD9 + 1] = |
a3c937ce AJ |
256 | { |
257 | -2.306006080437656357167128541231915480393E13L, | |
258 | -9.183606842453274924895648863832233799950E12L, | |
259 | -1.461857965935942962087907301194381010380E12L, | |
260 | -1.185728254682789754150068652663124298303E11L, | |
261 | -5.166285094703468567389566085480783070037E9L, | |
262 | -1.164573656694603024184768200787835094317E8L, | |
263 | -1.177343939483908678474886454113163527909E6L, | |
264 | -3.529391059783109732159524500029157638736E3L | |
265 | /* 1.0E0L */ | |
266 | }; | |
267 | ||
268 | ||
269 | /* log gamma(x+8) = log gamma(8) + x P(x)/Q(x) | |
270 | -0.5 <= x <= 0.5 | |
271 | 7.5 <= x+8 <= 8.5 | |
272 | Peak relative error 2.4e-37 */ | |
1f5649f8 UD |
273 | static const long double lgam8a = 8.525146484375E0L; |
274 | static const long double lgam8b = 1.4876690414300165531036347125050759667737E-5L; | |
a3c937ce | 275 | #define NRN8 8 |
1f5649f8 | 276 | static const long double RN8[NRN8 + 1] = |
a3c937ce AJ |
277 | { |
278 | 6.600775438203423546565361176829139703289E11L, | |
279 | 3.406361267593790705240802723914281025800E11L, | |
280 | 7.222460928505293914746983300555538432830E10L, | |
281 | 8.102984106025088123058747466840656458342E9L, | |
282 | 5.157620015986282905232150979772409345927E8L, | |
283 | 1.851445288272645829028129389609068641517E7L, | |
284 | 3.489261702223124354745894067468953756656E5L, | |
285 | 2.892095396706665774434217489775617756014E3L, | |
286 | 6.596977510622195827183948478627058738034E0L | |
287 | }; | |
288 | #define NRD8 7 | |
1f5649f8 | 289 | static const long double RD8[NRD8 + 1] = |
a3c937ce AJ |
290 | { |
291 | 3.274776546520735414638114828622673016920E11L, | |
292 | 1.581811207929065544043963828487733970107E11L, | |
293 | 3.108725655667825188135393076860104546416E10L, | |
294 | 3.193055010502912617128480163681842165730E9L, | |
295 | 1.830871482669835106357529710116211541839E8L, | |
296 | 5.790862854275238129848491555068073485086E6L, | |
297 | 9.305213264307921522842678835618803553589E4L, | |
298 | 6.216974105861848386918949336819572333622E2L | |
299 | /* 1.0E0L */ | |
300 | }; | |
301 | ||
302 | ||
303 | /* log gamma(x+7) = log gamma(7) + x P(x)/Q(x) | |
304 | -0.5 <= x <= 0.5 | |
305 | 6.5 <= x+7 <= 7.5 | |
306 | Peak relative error 3.2e-36 */ | |
1f5649f8 UD |
307 | static const long double lgam7a = 6.5792388916015625E0L; |
308 | static const long double lgam7b = 1.2320408538495060178292903945321122583007E-5L; | |
a3c937ce | 309 | #define NRN7 8 |
1f5649f8 | 310 | static const long double RN7[NRN7 + 1] = |
a3c937ce AJ |
311 | { |
312 | 2.065019306969459407636744543358209942213E11L, | |
313 | 1.226919919023736909889724951708796532847E11L, | |
314 | 2.996157990374348596472241776917953749106E10L, | |
315 | 3.873001919306801037344727168434909521030E9L, | |
316 | 2.841575255593761593270885753992732145094E8L, | |
317 | 1.176342515359431913664715324652399565551E7L, | |
318 | 2.558097039684188723597519300356028511547E5L, | |
319 | 2.448525238332609439023786244782810774702E3L, | |
320 | 6.460280377802030953041566617300902020435E0L | |
321 | }; | |
322 | #define NRD7 7 | |
1f5649f8 | 323 | static const long double RD7[NRD7 + 1] = |
a3c937ce AJ |
324 | { |
325 | 1.102646614598516998880874785339049304483E11L, | |
326 | 6.099297512712715445879759589407189290040E10L, | |
327 | 1.372898136289611312713283201112060238351E10L, | |
328 | 1.615306270420293159907951633566635172343E9L, | |
329 | 1.061114435798489135996614242842561967459E8L, | |
330 | 3.845638971184305248268608902030718674691E6L, | |
331 | 7.081730675423444975703917836972720495507E4L, | |
332 | 5.423122582741398226693137276201344096370E2L | |
333 | /* 1.0E0L */ | |
334 | }; | |
335 | ||
336 | ||
337 | /* log gamma(x+6) = log gamma(6) + x P(x)/Q(x) | |
338 | -0.5 <= x <= 0.5 | |
339 | 5.5 <= x+6 <= 6.5 | |
340 | Peak relative error 6.2e-37 */ | |
1f5649f8 UD |
341 | static const long double lgam6a = 4.7874908447265625E0L; |
342 | static const long double lgam6b = 8.9805548349424770093452324304839959231517E-7L; | |
a3c937ce | 343 | #define NRN6 8 |
1f5649f8 | 344 | static const long double RN6[NRN6 + 1] = |
a3c937ce AJ |
345 | { |
346 | -3.538412754670746879119162116819571823643E13L, | |
347 | -2.613432593406849155765698121483394257148E13L, | |
348 | -8.020670732770461579558867891923784753062E12L, | |
349 | -1.322227822931250045347591780332435433420E12L, | |
350 | -1.262809382777272476572558806855377129513E11L, | |
351 | -7.015006277027660872284922325741197022467E9L, | |
352 | -2.149320689089020841076532186783055727299E8L, | |
353 | -3.167210585700002703820077565539658995316E6L, | |
354 | -1.576834867378554185210279285358586385266E4L | |
355 | }; | |
356 | #define NRD6 8 | |
1f5649f8 | 357 | static const long double RD6[NRD6 + 1] = |
a3c937ce AJ |
358 | { |
359 | -2.073955870771283609792355579558899389085E13L, | |
360 | -1.421592856111673959642750863283919318175E13L, | |
361 | -4.012134994918353924219048850264207074949E12L, | |
362 | -6.013361045800992316498238470888523722431E11L, | |
363 | -5.145382510136622274784240527039643430628E10L, | |
364 | -2.510575820013409711678540476918249524123E9L, | |
365 | -6.564058379709759600836745035871373240904E7L, | |
366 | -7.861511116647120540275354855221373571536E5L, | |
367 | -2.821943442729620524365661338459579270561E3L | |
368 | /* 1.0E0L */ | |
369 | }; | |
370 | ||
371 | ||
372 | /* log gamma(x+5) = log gamma(5) + x P(x)/Q(x) | |
373 | -0.5 <= x <= 0.5 | |
374 | 4.5 <= x+5 <= 5.5 | |
375 | Peak relative error 3.4e-37 */ | |
1f5649f8 UD |
376 | static const long double lgam5a = 3.17803955078125E0L; |
377 | static const long double lgam5b = 1.4279566695619646941601297055408873990961E-5L; | |
a3c937ce | 378 | #define NRN5 9 |
1f5649f8 | 379 | static const long double RN5[NRN5 + 1] = |
a3c937ce AJ |
380 | { |
381 | 2.010952885441805899580403215533972172098E11L, | |
382 | 1.916132681242540921354921906708215338584E11L, | |
383 | 7.679102403710581712903937970163206882492E10L, | |
384 | 1.680514903671382470108010973615268125169E10L, | |
385 | 2.181011222911537259440775283277711588410E9L, | |
386 | 1.705361119398837808244780667539728356096E8L, | |
387 | 7.792391565652481864976147945997033946360E6L, | |
388 | 1.910741381027985291688667214472560023819E5L, | |
389 | 2.088138241893612679762260077783794329559E3L, | |
390 | 6.330318119566998299106803922739066556550E0L | |
391 | }; | |
392 | #define NRD5 8 | |
1f5649f8 | 393 | static const long double RD5[NRD5 + 1] = |
a3c937ce AJ |
394 | { |
395 | 1.335189758138651840605141370223112376176E11L, | |
396 | 1.174130445739492885895466097516530211283E11L, | |
397 | 4.308006619274572338118732154886328519910E10L, | |
398 | 8.547402888692578655814445003283720677468E9L, | |
399 | 9.934628078575618309542580800421370730906E8L, | |
400 | 6.847107420092173812998096295422311820672E7L, | |
401 | 2.698552646016599923609773122139463150403E6L, | |
402 | 5.526516251532464176412113632726150253215E4L, | |
403 | 4.772343321713697385780533022595450486932E2L | |
404 | /* 1.0E0L */ | |
405 | }; | |
406 | ||
407 | ||
408 | /* log gamma(x+4) = log gamma(4) + x P(x)/Q(x) | |
409 | -0.5 <= x <= 0.5 | |
410 | 3.5 <= x+4 <= 4.5 | |
411 | Peak relative error 6.7e-37 */ | |
1f5649f8 UD |
412 | static const long double lgam4a = 1.791748046875E0L; |
413 | static const long double lgam4b = 1.1422353055000812477358380702272722990692E-5L; | |
a3c937ce | 414 | #define NRN4 9 |
1f5649f8 | 415 | static const long double RN4[NRN4 + 1] = |
a3c937ce AJ |
416 | { |
417 | -1.026583408246155508572442242188887829208E13L, | |
418 | -1.306476685384622809290193031208776258809E13L, | |
419 | -7.051088602207062164232806511992978915508E12L, | |
420 | -2.100849457735620004967624442027793656108E12L, | |
421 | -3.767473790774546963588549871673843260569E11L, | |
422 | -4.156387497364909963498394522336575984206E10L, | |
423 | -2.764021460668011732047778992419118757746E9L, | |
424 | -1.036617204107109779944986471142938641399E8L, | |
425 | -1.895730886640349026257780896972598305443E6L, | |
426 | -1.180509051468390914200720003907727988201E4L | |
427 | }; | |
428 | #define NRD4 9 | |
1f5649f8 | 429 | static const long double RD4[NRD4 + 1] = |
a3c937ce AJ |
430 | { |
431 | -8.172669122056002077809119378047536240889E12L, | |
432 | -9.477592426087986751343695251801814226960E12L, | |
433 | -4.629448850139318158743900253637212801682E12L, | |
434 | -1.237965465892012573255370078308035272942E12L, | |
435 | -1.971624313506929845158062177061297598956E11L, | |
436 | -1.905434843346570533229942397763361493610E10L, | |
437 | -1.089409357680461419743730978512856675984E9L, | |
438 | -3.416703082301143192939774401370222822430E7L, | |
439 | -4.981791914177103793218433195857635265295E5L, | |
440 | -2.192507743896742751483055798411231453733E3L | |
441 | /* 1.0E0L */ | |
442 | }; | |
443 | ||
444 | ||
445 | /* log gamma(x+3) = log gamma(3) + x P(x)/Q(x) | |
446 | -0.25 <= x <= 0.5 | |
447 | 2.75 <= x+3 <= 3.5 | |
448 | Peak relative error 6.0e-37 */ | |
1f5649f8 UD |
449 | static const long double lgam3a = 6.93145751953125E-1L; |
450 | static const long double lgam3b = 1.4286068203094172321214581765680755001344E-6L; | |
a3c937ce AJ |
451 | |
452 | #define NRN3 9 | |
1f5649f8 | 453 | static const long double RN3[NRN3 + 1] = |
a3c937ce AJ |
454 | { |
455 | -4.813901815114776281494823863935820876670E11L, | |
456 | -8.425592975288250400493910291066881992620E11L, | |
457 | -6.228685507402467503655405482985516909157E11L, | |
458 | -2.531972054436786351403749276956707260499E11L, | |
459 | -6.170200796658926701311867484296426831687E10L, | |
460 | -9.211477458528156048231908798456365081135E9L, | |
461 | -8.251806236175037114064561038908691305583E8L, | |
462 | -4.147886355917831049939930101151160447495E7L, | |
463 | -1.010851868928346082547075956946476932162E6L, | |
464 | -8.333374463411801009783402800801201603736E3L | |
465 | }; | |
466 | #define NRD3 9 | |
1f5649f8 | 467 | static const long double RD3[NRD3 + 1] = |
a3c937ce AJ |
468 | { |
469 | -5.216713843111675050627304523368029262450E11L, | |
470 | -8.014292925418308759369583419234079164391E11L, | |
471 | -5.180106858220030014546267824392678611990E11L, | |
472 | -1.830406975497439003897734969120997840011E11L, | |
473 | -3.845274631904879621945745960119924118925E10L, | |
474 | -4.891033385370523863288908070309417710903E9L, | |
475 | -3.670172254411328640353855768698287474282E8L, | |
476 | -1.505316381525727713026364396635522516989E7L, | |
477 | -2.856327162923716881454613540575964890347E5L, | |
478 | -1.622140448015769906847567212766206894547E3L | |
479 | /* 1.0E0L */ | |
480 | }; | |
481 | ||
482 | ||
483 | /* log gamma(x+2.5) = log gamma(2.5) + x P(x)/Q(x) | |
484 | -0.125 <= x <= 0.25 | |
485 | 2.375 <= x+2.5 <= 2.75 */ | |
1f5649f8 UD |
486 | static const long double lgam2r5a = 2.8466796875E-1L; |
487 | static const long double lgam2r5b = 1.4901722919159632494669682701924320137696E-5L; | |
a3c937ce | 488 | #define NRN2r5 8 |
1f5649f8 | 489 | static const long double RN2r5[NRN2r5 + 1] = |
a3c937ce AJ |
490 | { |
491 | -4.676454313888335499356699817678862233205E9L, | |
492 | -9.361888347911187924389905984624216340639E9L, | |
493 | -7.695353600835685037920815799526540237703E9L, | |
494 | -3.364370100981509060441853085968900734521E9L, | |
495 | -8.449902011848163568670361316804900559863E8L, | |
496 | -1.225249050950801905108001246436783022179E8L, | |
497 | -9.732972931077110161639900388121650470926E6L, | |
498 | -3.695711763932153505623248207576425983573E5L, | |
499 | -4.717341584067827676530426007495274711306E3L | |
500 | }; | |
501 | #define NRD2r5 8 | |
1f5649f8 | 502 | static const long double RD2r5[NRD2r5 + 1] = |
a3c937ce AJ |
503 | { |
504 | -6.650657966618993679456019224416926875619E9L, | |
505 | -1.099511409330635807899718829033488771623E10L, | |
506 | -7.482546968307837168164311101447116903148E9L, | |
507 | -2.702967190056506495988922973755870557217E9L, | |
508 | -5.570008176482922704972943389590409280950E8L, | |
509 | -6.536934032192792470926310043166993233231E7L, | |
510 | -4.101991193844953082400035444146067511725E6L, | |
511 | -1.174082735875715802334430481065526664020E5L, | |
512 | -9.932840389994157592102947657277692978511E2L | |
513 | /* 1.0E0L */ | |
514 | }; | |
515 | ||
516 | ||
517 | /* log gamma(x+2) = x P(x)/Q(x) | |
518 | -0.125 <= x <= +0.375 | |
519 | 1.875 <= x+2 <= 2.375 | |
520 | Peak relative error 4.6e-36 */ | |
521 | #define NRN2 9 | |
1f5649f8 | 522 | static const long double RN2[NRN2 + 1] = |
a3c937ce AJ |
523 | { |
524 | -3.716661929737318153526921358113793421524E9L, | |
525 | -1.138816715030710406922819131397532331321E10L, | |
526 | -1.421017419363526524544402598734013569950E10L, | |
527 | -9.510432842542519665483662502132010331451E9L, | |
528 | -3.747528562099410197957514973274474767329E9L, | |
529 | -8.923565763363912474488712255317033616626E8L, | |
530 | -1.261396653700237624185350402781338231697E8L, | |
531 | -9.918402520255661797735331317081425749014E6L, | |
532 | -3.753996255897143855113273724233104768831E5L, | |
533 | -4.778761333044147141559311805999540765612E3L | |
534 | }; | |
535 | #define NRD2 9 | |
1f5649f8 | 536 | static const long double RD2[NRD2 + 1] = |
a3c937ce AJ |
537 | { |
538 | -8.790916836764308497770359421351673950111E9L, | |
539 | -2.023108608053212516399197678553737477486E10L, | |
540 | -1.958067901852022239294231785363504458367E10L, | |
541 | -1.035515043621003101254252481625188704529E10L, | |
542 | -3.253884432621336737640841276619272224476E9L, | |
543 | -6.186383531162456814954947669274235815544E8L, | |
544 | -6.932557847749518463038934953605969951466E7L, | |
545 | -4.240731768287359608773351626528479703758E6L, | |
546 | -1.197343995089189188078944689846348116630E5L, | |
547 | -1.004622911670588064824904487064114090920E3L | |
548 | /* 1.0E0 */ | |
549 | }; | |
550 | ||
551 | ||
552 | /* log gamma(x+1.75) = log gamma(1.75) + x P(x)/Q(x) | |
553 | -0.125 <= x <= +0.125 | |
554 | 1.625 <= x+1.75 <= 1.875 | |
555 | Peak relative error 9.2e-37 */ | |
1f5649f8 UD |
556 | static const long double lgam1r75a = -8.441162109375E-2L; |
557 | static const long double lgam1r75b = 1.0500073264444042213965868602268256157604E-5L; | |
a3c937ce | 558 | #define NRN1r75 8 |
1f5649f8 | 559 | static const long double RN1r75[NRN1r75 + 1] = |
a3c937ce AJ |
560 | { |
561 | -5.221061693929833937710891646275798251513E7L, | |
562 | -2.052466337474314812817883030472496436993E8L, | |
563 | -2.952718275974940270675670705084125640069E8L, | |
564 | -2.132294039648116684922965964126389017840E8L, | |
565 | -8.554103077186505960591321962207519908489E7L, | |
566 | -1.940250901348870867323943119132071960050E7L, | |
567 | -2.379394147112756860769336400290402208435E6L, | |
568 | -1.384060879999526222029386539622255797389E5L, | |
569 | -2.698453601378319296159355612094598695530E3L | |
570 | }; | |
571 | #define NRD1r75 8 | |
1f5649f8 | 572 | static const long double RD1r75[NRD1r75 + 1] = |
a3c937ce AJ |
573 | { |
574 | -2.109754689501705828789976311354395393605E8L, | |
575 | -5.036651829232895725959911504899241062286E8L, | |
576 | -4.954234699418689764943486770327295098084E8L, | |
577 | -2.589558042412676610775157783898195339410E8L, | |
578 | -7.731476117252958268044969614034776883031E7L, | |
579 | -1.316721702252481296030801191240867486965E7L, | |
580 | -1.201296501404876774861190604303728810836E6L, | |
581 | -5.007966406976106636109459072523610273928E4L, | |
582 | -6.155817990560743422008969155276229018209E2L | |
583 | /* 1.0E0L */ | |
584 | }; | |
585 | ||
586 | ||
587 | /* log gamma(x+x0) = y0 + x^2 P(x)/Q(x) | |
588 | -0.0867 <= x <= +0.1634 | |
589 | 1.374932... <= x+x0 <= 1.625032... | |
590 | Peak relative error 4.0e-36 */ | |
1f5649f8 UD |
591 | static const long double x0a = 1.4616241455078125L; |
592 | static const long double x0b = 7.9994605498412626595423257213002588621246E-6L; | |
593 | static const long double y0a = -1.21490478515625E-1L; | |
594 | static const long double y0b = 4.1879797753919044854428223084178486438269E-6L; | |
a3c937ce | 595 | #define NRN1r5 8 |
1f5649f8 | 596 | static const long double RN1r5[NRN1r5 + 1] = |
a3c937ce AJ |
597 | { |
598 | 6.827103657233705798067415468881313128066E5L, | |
599 | 1.910041815932269464714909706705242148108E6L, | |
600 | 2.194344176925978377083808566251427771951E6L, | |
601 | 1.332921400100891472195055269688876427962E6L, | |
602 | 4.589080973377307211815655093824787123508E5L, | |
603 | 8.900334161263456942727083580232613796141E4L, | |
604 | 9.053840838306019753209127312097612455236E3L, | |
605 | 4.053367147553353374151852319743594873771E2L, | |
606 | 5.040631576303952022968949605613514584950E0L | |
607 | }; | |
608 | #define NRD1r5 8 | |
1f5649f8 | 609 | static const long double RD1r5[NRD1r5 + 1] = |
a3c937ce AJ |
610 | { |
611 | 1.411036368843183477558773688484699813355E6L, | |
612 | 4.378121767236251950226362443134306184849E6L, | |
613 | 5.682322855631723455425929877581697918168E6L, | |
614 | 3.999065731556977782435009349967042222375E6L, | |
615 | 1.653651390456781293163585493620758410333E6L, | |
616 | 4.067774359067489605179546964969435858311E5L, | |
617 | 5.741463295366557346748361781768833633256E4L, | |
618 | 4.226404539738182992856094681115746692030E3L, | |
619 | 1.316980975410327975566999780608618774469E2L, | |
620 | /* 1.0E0L */ | |
621 | }; | |
622 | ||
623 | ||
624 | /* log gamma(x+1.25) = log gamma(1.25) + x P(x)/Q(x) | |
625 | -.125 <= x <= +.125 | |
626 | 1.125 <= x+1.25 <= 1.375 | |
627 | Peak relative error = 4.9e-36 */ | |
1f5649f8 UD |
628 | static const long double lgam1r25a = -9.82818603515625E-2L; |
629 | static const long double lgam1r25b = 1.0023929749338536146197303364159774377296E-5L; | |
a3c937ce | 630 | #define NRN1r25 9 |
1f5649f8 | 631 | static const long double RN1r25[NRN1r25 + 1] = |
a3c937ce AJ |
632 | { |
633 | -9.054787275312026472896002240379580536760E4L, | |
634 | -8.685076892989927640126560802094680794471E4L, | |
635 | 2.797898965448019916967849727279076547109E5L, | |
636 | 6.175520827134342734546868356396008898299E5L, | |
637 | 5.179626599589134831538516906517372619641E5L, | |
638 | 2.253076616239043944538380039205558242161E5L, | |
639 | 5.312653119599957228630544772499197307195E4L, | |
640 | 6.434329437514083776052669599834938898255E3L, | |
641 | 3.385414416983114598582554037612347549220E2L, | |
642 | 4.907821957946273805080625052510832015792E0L | |
643 | }; | |
644 | #define NRD1r25 8 | |
1f5649f8 | 645 | static const long double RD1r25[NRD1r25 + 1] = |
a3c937ce AJ |
646 | { |
647 | 3.980939377333448005389084785896660309000E5L, | |
648 | 1.429634893085231519692365775184490465542E6L, | |
649 | 2.145438946455476062850151428438668234336E6L, | |
650 | 1.743786661358280837020848127465970357893E6L, | |
651 | 8.316364251289743923178092656080441655273E5L, | |
652 | 2.355732939106812496699621491135458324294E5L, | |
653 | 3.822267399625696880571810137601310855419E4L, | |
654 | 3.228463206479133236028576845538387620856E3L, | |
655 | 1.152133170470059555646301189220117965514E2L | |
656 | /* 1.0E0L */ | |
657 | }; | |
658 | ||
659 | ||
660 | /* log gamma(x + 1) = x P(x)/Q(x) | |
661 | 0.0 <= x <= +0.125 | |
662 | 1.0 <= x+1 <= 1.125 | |
663 | Peak relative error 1.1e-35 */ | |
664 | #define NRN1 8 | |
1f5649f8 | 665 | static const long double RN1[NRN1 + 1] = |
a3c937ce AJ |
666 | { |
667 | -9.987560186094800756471055681088744738818E3L, | |
668 | -2.506039379419574361949680225279376329742E4L, | |
669 | -1.386770737662176516403363873617457652991E4L, | |
670 | 1.439445846078103202928677244188837130744E4L, | |
671 | 2.159612048879650471489449668295139990693E4L, | |
672 | 1.047439813638144485276023138173676047079E4L, | |
673 | 2.250316398054332592560412486630769139961E3L, | |
674 | 1.958510425467720733041971651126443864041E2L, | |
675 | 4.516830313569454663374271993200291219855E0L | |
676 | }; | |
677 | #define NRD1 7 | |
1f5649f8 | 678 | static const long double RD1[NRD1 + 1] = |
a3c937ce AJ |
679 | { |
680 | 1.730299573175751778863269333703788214547E4L, | |
681 | 6.807080914851328611903744668028014678148E4L, | |
682 | 1.090071629101496938655806063184092302439E5L, | |
683 | 9.124354356415154289343303999616003884080E4L, | |
684 | 4.262071638655772404431164427024003253954E4L, | |
685 | 1.096981664067373953673982635805821283581E4L, | |
686 | 1.431229503796575892151252708527595787588E3L, | |
687 | 7.734110684303689320830401788262295992921E1L | |
688 | /* 1.0E0 */ | |
689 | }; | |
690 | ||
691 | ||
692 | /* log gamma(x + 1) = x P(x)/Q(x) | |
693 | -0.125 <= x <= 0 | |
694 | 0.875 <= x+1 <= 1.0 | |
695 | Peak relative error 7.0e-37 */ | |
696 | #define NRNr9 8 | |
1f5649f8 | 697 | static const long double RNr9[NRNr9 + 1] = |
a3c937ce AJ |
698 | { |
699 | 4.441379198241760069548832023257571176884E5L, | |
700 | 1.273072988367176540909122090089580368732E6L, | |
701 | 9.732422305818501557502584486510048387724E5L, | |
702 | -5.040539994443998275271644292272870348684E5L, | |
703 | -1.208719055525609446357448132109723786736E6L, | |
704 | -7.434275365370936547146540554419058907156E5L, | |
705 | -2.075642969983377738209203358199008185741E5L, | |
706 | -2.565534860781128618589288075109372218042E4L, | |
707 | -1.032901669542994124131223797515913955938E3L, | |
708 | }; | |
709 | #define NRDr9 8 | |
1f5649f8 | 710 | static const long double RDr9[NRDr9 + 1] = |
a3c937ce AJ |
711 | { |
712 | -7.694488331323118759486182246005193998007E5L, | |
713 | -3.301918855321234414232308938454112213751E6L, | |
714 | -5.856830900232338906742924836032279404702E6L, | |
715 | -5.540672519616151584486240871424021377540E6L, | |
716 | -3.006530901041386626148342989181721176919E6L, | |
717 | -9.350378280513062139466966374330795935163E5L, | |
718 | -1.566179100031063346901755685375732739511E5L, | |
719 | -1.205016539620260779274902967231510804992E4L, | |
720 | -2.724583156305709733221564484006088794284E2L | |
721 | /* 1.0E0 */ | |
722 | }; | |
723 | ||
724 | ||
725 | /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ | |
726 | ||
727 | static long double | |
60a06b7c | 728 | neval (long double x, const long double *p, int n) |
a3c937ce AJ |
729 | { |
730 | long double y; | |
731 | ||
732 | p += n; | |
733 | y = *p--; | |
734 | do | |
735 | { | |
736 | y = y * x + *p--; | |
737 | } | |
738 | while (--n > 0); | |
739 | return y; | |
740 | } | |
741 | ||
742 | ||
743 | /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ | |
744 | ||
745 | static long double | |
60a06b7c | 746 | deval (long double x, const long double *p, int n) |
a3c937ce AJ |
747 | { |
748 | long double y; | |
749 | ||
750 | p += n; | |
751 | y = x + *p--; | |
752 | do | |
753 | { | |
754 | y = y * x + *p--; | |
755 | } | |
756 | while (--n > 0); | |
757 | return y; | |
758 | } | |
759 | ||
760 | ||
a3c937ce AJ |
761 | long double |
762 | __ieee754_lgammal_r (long double x, int *signgamp) | |
a3c937ce AJ |
763 | { |
764 | long double p, q, w, z, nx; | |
765 | int i, nn; | |
766 | ||
767 | *signgamp = 1; | |
768 | ||
d81f90cc | 769 | if (! isfinite (x)) |
a3c937ce AJ |
770 | return x * x; |
771 | ||
facd1d8e UD |
772 | if (x == 0.0L) |
773 | { | |
d81f90cc | 774 | if (signbit (x)) |
0ac5ae23 | 775 | *signgamp = -1; |
facd1d8e UD |
776 | } |
777 | ||
a3c937ce AJ |
778 | if (x < 0.0L) |
779 | { | |
050f29c1 JM |
780 | if (x < -2.0L && x > (LDBL_MANT_DIG == 106 ? -48.0L : -50.0L)) |
781 | return __lgamma_negl (x, signgamp); | |
a3c937ce | 782 | q = -x; |
a3c937ce AJ |
783 | p = __floorl (q); |
784 | if (p == q) | |
52e1b618 | 785 | return (one / (p - p)); |
9bb69b60 JM |
786 | long double halfp = p * 0.5L; |
787 | if (halfp == __floorl (halfp)) | |
a3c937ce AJ |
788 | *signgamp = -1; |
789 | else | |
790 | *signgamp = 1; | |
4f40e4b3 JM |
791 | if (q < 0x1p-120L) |
792 | return -__logl (q); | |
a3c937ce AJ |
793 | z = q - p; |
794 | if (z > 0.5L) | |
795 | { | |
796 | p += 1.0L; | |
797 | z = p - q; | |
798 | } | |
799 | z = q * __sinl (PIL * z); | |
52e1b618 | 800 | w = __ieee754_lgammal_r (q, &i); |
a3c937ce AJ |
801 | z = __logl (PIL / z) - w; |
802 | return (z); | |
803 | } | |
804 | ||
805 | if (x < 13.5L) | |
806 | { | |
807 | p = 0.0L; | |
808 | nx = __floorl (x + 0.5L); | |
809 | nn = nx; | |
810 | switch (nn) | |
811 | { | |
812 | case 0: | |
813 | /* log gamma (x + 1) = log(x) + log gamma(x) */ | |
eb3fc44b JM |
814 | if (x < 0x1p-120L) |
815 | return -__logl (x); | |
816 | else if (x <= 0.125) | |
a3c937ce AJ |
817 | { |
818 | p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1); | |
819 | } | |
820 | else if (x <= 0.375) | |
821 | { | |
822 | z = x - 0.25L; | |
823 | p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25); | |
824 | p += lgam1r25b; | |
825 | p += lgam1r25a; | |
826 | } | |
827 | else if (x <= 0.625) | |
828 | { | |
829 | z = x + (1.0L - x0a); | |
830 | z = z - x0b; | |
831 | p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5); | |
832 | p = p * z * z; | |
833 | p = p + y0b; | |
834 | p = p + y0a; | |
835 | } | |
836 | else if (x <= 0.875) | |
837 | { | |
838 | z = x - 0.75L; | |
839 | p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75); | |
840 | p += lgam1r75b; | |
841 | p += lgam1r75a; | |
842 | } | |
843 | else | |
844 | { | |
845 | z = x - 1.0L; | |
846 | p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2); | |
847 | } | |
848 | p = p - __logl (x); | |
849 | break; | |
850 | ||
851 | case 1: | |
852 | if (x < 0.875L) | |
853 | { | |
854 | if (x <= 0.625) | |
855 | { | |
856 | z = x + (1.0L - x0a); | |
857 | z = z - x0b; | |
858 | p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5); | |
859 | p = p * z * z; | |
860 | p = p + y0b; | |
861 | p = p + y0a; | |
862 | } | |
863 | else if (x <= 0.875) | |
864 | { | |
865 | z = x - 0.75L; | |
866 | p = z * neval (z, RN1r75, NRN1r75) | |
0ac5ae23 | 867 | / deval (z, RD1r75, NRD1r75); |
a3c937ce AJ |
868 | p += lgam1r75b; |
869 | p += lgam1r75a; | |
870 | } | |
871 | else | |
872 | { | |
873 | z = x - 1.0L; | |
874 | p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2); | |
875 | } | |
876 | p = p - __logl (x); | |
877 | } | |
878 | else if (x < 1.0L) | |
879 | { | |
880 | z = x - 1.0L; | |
881 | p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9); | |
882 | } | |
52e1b618 UD |
883 | else if (x == 1.0L) |
884 | p = 0.0L; | |
a3c937ce AJ |
885 | else if (x <= 1.125L) |
886 | { | |
887 | z = x - 1.0L; | |
888 | p = z * neval (z, RN1, NRN1) / deval (z, RD1, NRD1); | |
889 | } | |
890 | else if (x <= 1.375) | |
891 | { | |
892 | z = x - 1.25L; | |
893 | p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25); | |
894 | p += lgam1r25b; | |
895 | p += lgam1r25a; | |
896 | } | |
897 | else | |
898 | { | |
899 | /* 1.375 <= x+x0 <= 1.625 */ | |
900 | z = x - x0a; | |
901 | z = z - x0b; | |
902 | p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5); | |
903 | p = p * z * z; | |
904 | p = p + y0b; | |
905 | p = p + y0a; | |
906 | } | |
907 | break; | |
908 | ||
909 | case 2: | |
910 | if (x < 1.625L) | |
911 | { | |
912 | z = x - x0a; | |
913 | z = z - x0b; | |
914 | p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5); | |
915 | p = p * z * z; | |
916 | p = p + y0b; | |
917 | p = p + y0a; | |
918 | } | |
919 | else if (x < 1.875L) | |
920 | { | |
921 | z = x - 1.75L; | |
922 | p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75); | |
923 | p += lgam1r75b; | |
924 | p += lgam1r75a; | |
925 | } | |
52e1b618 UD |
926 | else if (x == 2.0L) |
927 | p = 0.0L; | |
a3c937ce AJ |
928 | else if (x < 2.375L) |
929 | { | |
930 | z = x - 2.0L; | |
931 | p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2); | |
932 | } | |
933 | else | |
934 | { | |
935 | z = x - 2.5L; | |
936 | p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5); | |
937 | p += lgam2r5b; | |
938 | p += lgam2r5a; | |
939 | } | |
940 | break; | |
941 | ||
942 | case 3: | |
943 | if (x < 2.75) | |
944 | { | |
945 | z = x - 2.5L; | |
946 | p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5); | |
947 | p += lgam2r5b; | |
948 | p += lgam2r5a; | |
949 | } | |
950 | else | |
951 | { | |
952 | z = x - 3.0L; | |
953 | p = z * neval (z, RN3, NRN3) / deval (z, RD3, NRD3); | |
954 | p += lgam3b; | |
955 | p += lgam3a; | |
956 | } | |
957 | break; | |
958 | ||
959 | case 4: | |
960 | z = x - 4.0L; | |
961 | p = z * neval (z, RN4, NRN4) / deval (z, RD4, NRD4); | |
962 | p += lgam4b; | |
963 | p += lgam4a; | |
964 | break; | |
965 | ||
966 | case 5: | |
967 | z = x - 5.0L; | |
968 | p = z * neval (z, RN5, NRN5) / deval (z, RD5, NRD5); | |
969 | p += lgam5b; | |
970 | p += lgam5a; | |
971 | break; | |
972 | ||
973 | case 6: | |
974 | z = x - 6.0L; | |
975 | p = z * neval (z, RN6, NRN6) / deval (z, RD6, NRD6); | |
976 | p += lgam6b; | |
977 | p += lgam6a; | |
978 | break; | |
979 | ||
980 | case 7: | |
981 | z = x - 7.0L; | |
982 | p = z * neval (z, RN7, NRN7) / deval (z, RD7, NRD7); | |
983 | p += lgam7b; | |
984 | p += lgam7a; | |
985 | break; | |
986 | ||
987 | case 8: | |
988 | z = x - 8.0L; | |
989 | p = z * neval (z, RN8, NRN8) / deval (z, RD8, NRD8); | |
990 | p += lgam8b; | |
991 | p += lgam8a; | |
992 | break; | |
993 | ||
994 | case 9: | |
995 | z = x - 9.0L; | |
996 | p = z * neval (z, RN9, NRN9) / deval (z, RD9, NRD9); | |
997 | p += lgam9b; | |
998 | p += lgam9a; | |
999 | break; | |
1000 | ||
1001 | case 10: | |
1002 | z = x - 10.0L; | |
1003 | p = z * neval (z, RN10, NRN10) / deval (z, RD10, NRD10); | |
1004 | p += lgam10b; | |
1005 | p += lgam10a; | |
1006 | break; | |
1007 | ||
1008 | case 11: | |
1009 | z = x - 11.0L; | |
1010 | p = z * neval (z, RN11, NRN11) / deval (z, RD11, NRD11); | |
1011 | p += lgam11b; | |
1012 | p += lgam11a; | |
1013 | break; | |
1014 | ||
1015 | case 12: | |
1016 | z = x - 12.0L; | |
1017 | p = z * neval (z, RN12, NRN12) / deval (z, RD12, NRD12); | |
1018 | p += lgam12b; | |
1019 | p += lgam12a; | |
1020 | break; | |
1021 | ||
1022 | case 13: | |
1023 | z = x - 13.0L; | |
1024 | p = z * neval (z, RN13, NRN13) / deval (z, RD13, NRD13); | |
1025 | p += lgam13b; | |
1026 | p += lgam13a; | |
1027 | break; | |
1028 | } | |
1029 | return p; | |
1030 | } | |
1031 | ||
1032 | if (x > MAXLGM) | |
1033 | return (*signgamp * huge * huge); | |
1034 | ||
4b84e247 JM |
1035 | if (x > 0x1p120L) |
1036 | return x * (__logl (x) - 1.0L); | |
a3c937ce AJ |
1037 | q = ls2pi - x; |
1038 | q = (x - 0.5L) * __logl (x) + q; | |
1039 | if (x > 1.0e18L) | |
1040 | return (q); | |
1041 | ||
1042 | p = 1.0L / (x * x); | |
1043 | q += neval (p, RASY, NRASY) / x; | |
1044 | return (q); | |
1045 | } | |
0ac5ae23 | 1046 | strong_alias (__ieee754_lgammal_r, __lgammal_r_finite) |