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