1/**
2 * This file has no copyright assigned and is placed in the Public Domain.
3 * This file is part of the mingw-w64 runtime package.
4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
5 */
6#include "cephes_mconf.h"
7
8#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__
9double lgamma(double x);
10
11long double lgammal(long double x)
12{
13 return lgamma(x);
14}
15#else
16
17#if UNK
18static uLD S[9] = {
19 { { -1.193945051381510095614E-3L } },
20 { { 7.220599478036909672331E-3L } },
21 { { -9.622023360406271645744E-3L } },
22 { { -4.219773360705915470089E-2L } },
23 { { 1.665386113720805206758E-1L } },
24 { { -4.200263503403344054473E-2L } },
25 { { -6.558780715202540684668E-1L } },
26 { { 5.772156649015328608253E-1L } },
27 { { 1.000000000000000000000E0L } }
28};
29#endif
30#if IBMPC
31static const uLD S[] = {
32 { { 0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, 0, 0, 0 } },
33 { { 0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, 0, 0, 0 } },
34 { { 0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, 0, 0, 0 } },
35 { { 0x10b0,0xec17,0x87dc,0xacd7,0xbffa, 0, 0, 0 } },
36 { { 0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, 0, 0, 0 } },
37 { { 0xf183,0x126b,0xf47d,0xac0a,0xbffa, 0, 0, 0 } },
38 { { 0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, 0, 0, 0 } },
39 { { 0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, 0, 0, 0 } },
40 { { 0x0000,0x0000,0x0000,0x8000,0x3fff, 0, 0, 0 } }
41};
42#endif
43#if MIEEE
44static uLD S[27] = {
45 { { 0xbff50000,0x9c7e25e5,0xd6d3baeb, 0 } },
46 { { 0x3ff70000,0xec9ac74e,0xceb4fe9a, 0 } },
47 { { 0xbff80000,0x9da5b0e9,0xdfef9225, 0 } },
48 { { 0xbffa0000,0xacd787dc,0xec1710b0, 0 } },
49 { { 0x3ffc0000,0xaa891905,0x75156b8d, 0 } },
50 { { 0xbffa0000,0xac0af47d,0x126bf183, 0 } },
51 { { 0xbffe0000,0xa7e7a013,0x57d17bf6, 0 } },
52 { { 0x3ffe0000,0x93c467e3,0x7db0c7a9, 0 } },
53 { { 0x3fff0000,0x80000000,0x00000000, 0 } }
54};
55#endif
56
57#if UNK
58static uLD SN[9] = {
59 { { 1.133374167243894382010E-3L } },
60 { { 7.220837261893170325704E-3L } },
61 { { 9.621911155035976733706E-3L } },
62 { { -4.219773343731191721664E-2L } },
63 { { -1.665386113944413519335E-1L } },
64 { { -4.200263503402112910504E-2L } },
65 { { 6.558780715202536547116E-1L } },
66 { { 5.772156649015328608727E-1L } },
67 { { -1.000000000000000000000E0L } }
68};
69#endif
70#if IBMPC
71static const uLD SN[] = {
72 { { 0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, 0, 0, 0 } },
73 { { 0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, 0, 0, 0 } },
74 { { 0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, 0, 0, 0 } },
75 { { 0x783f,0x41dd,0x87d1,0xacd7,0xbffa, 0, 0, 0 } },
76 { { 0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, 0, 0, 0 } },
77 { { 0x7f64,0x1234,0xf47d,0xac0a,0xbffa, 0, 0, 0 } },
78 { { 0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, 0, 0, 0 } },
79 { { 0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, 0, 0, 0 } },
80 { { 0x0000,0x0000,0x0000,0x8000,0xbfff, 0, 0, 0 } }
81};
82#endif
83#if MIEEE
84static uLD SN[] = {
85 { { 0x3ff50000,0x948db9f7,0x02de5dd1, 0 } },
86 { { 0x3ff70000,0xec9cc5f1,0xdd68989b, 0 } },
87 { { 0x3ff80000,0x9da5386f,0x18f02ca1, 0 } },
88 { { 0xbffa0000,0xacd787d1,0x41dd783f, 0 } },
89 { { 0xbffc0000,0xaa891905,0xd76d7a5b, 0 } },
90 { { 0xbffa0000,0xac0af47d,0x12347f64, 0 } },
91 { { 0x3ffe0000,0xa7e7a013,0x57d15e26, 0 } },
92 { { 0x3ffe0000,0x93c467e3,0x7db0c7aa, 0 } },
93 { { 0xbfff0000,0x80000000,0x00000000, 0 } }
94};
95#endif
96
97
98/* A[]: Stirling's formula expansion of log gamma
99 * B[], C[]: log gamma function between 2 and 3
100 */
101
102
103/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2)
104 * x >= 8
105 * Peak relative error 1.51e-21
106 * Relative spread of error peaks 5.67e-21
107 */
108#if UNK
109static uLD A[7] = {
110 { { 4.885026142432270781165E-3L } },
111 { { -1.880801938119376907179E-3L } },
112 { { 8.412723297322498080632E-4L } },
113 { { -5.952345851765688514613E-4L } },
114 { { 7.936507795855070755671E-4L } },
115 { { -2.777777777750349603440E-3L } },
116 { { 8.333333333333331447505E-2L } }
117};
118#endif
119#if IBMPC
120static const uLD A[] = {
121 { { 0xd984,0xcc08,0x91c2,0xa012,0x3ff7, 0, 0, 0 } },
122 { { 0x3d91,0x0304,0x3da1,0xf685,0xbff5, 0, 0, 0 } },
123 { { 0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, 0, 0, 0 } },
124 { { 0x8b20,0x9fce,0x844e,0x9c09,0xbff4, 0, 0, 0 } },
125 { { 0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, 0, 0, 0 } },
126 { { 0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, 0, 0, 0 } },
127 { { 0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, 0, 0, 0 } }
128};
129#endif
130#if MIEEE
131static uLD A[] = {
132 { { 0x3ff70000,0xa01291c2,0xcc08d984, 0 } },
133 { { 0xbff50000,0xf6853da1,0x03043d91, 0 } },
134 { { 0x3ff40000,0xdc88d492,0xaad13bdc, 0 } },
135 { { 0xbff40000,0x9c09844e,0x9fce8b20, 0 } },
136 { { 0x3ff40000,0xd00d0092,0x30e5f8f2, 0 } },
137 { { 0xbff60000,0xb60b60b6,0x03a84d88, 0 } },
138 { { 0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc, 0 } }
139};
140#endif
141
142/* log gamma(x+2) = x B(x)/C(x)
143 * 0 <= x <= 1
144 * Peak relative error 7.16e-22
145 * Relative spread of error peaks 4.78e-20
146 */
147#if UNK
148static uLD B[7] = {
149 { { -2.163690827643812857640E3L } },
150 { { -8.723871522843511459790E4L } },
151 { { -1.104326814691464261197E6L } },
152 { { -6.111225012005214299996E6L } },
153 { { -1.625568062543700591014E7L } },
154 { { -2.003937418103815175475E7L } },
155 { { -8.875666783650703802159E6L } }
156};
157static uLD C[7] = {
158 { { -5.139481484435370143617E2L } },
159 { { -3.403570840534304670537E4L } },
160 { { -6.227441164066219501697E5L } },
161 { { -4.814940379411882186630E6L } },
162 { { -1.785433287045078156959E7L } },
163 { { -3.138646407656182662088E7L } },
164 { { -2.099336717757895876142E7L } }
165};
166#endif
167#if IBMPC
168static const uLD B[] = {
169 { { 0x9557,0x4995,0x0da1,0x873b,0xc00a, 0, 0, 0 } },
170 { { 0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, 0, 0, 0 } },
171 { { 0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, 0, 0, 0 } },
172 { { 0x259a,0x258c,0xf206,0xba7f,0xc015, 0, 0, 0 } },
173 { { 0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, 0, 0, 0 } },
174 { { 0x168f,0x2c42,0x6717,0x98e3,0xc017, 0, 0, 0 } },
175 { { 0x2051,0x9d55,0x92c8,0x876e,0xc016, 0, 0, 0 } }
176};
177static const uLD C[] = {
178 { { 0xaa77,0xcf2f,0xae76,0x807c,0xc008, 0, 0, 0 } },
179 { { 0xb280,0x0d74,0xb55a,0x84f3,0xc00e, 0, 0, 0 } },
180 { { 0xa505,0xcd30,0x81dc,0x9809,0xc012, 0, 0, 0 } },
181 { { 0x3369,0x4246,0xb8c2,0x92f0,0xc015, 0, 0, 0 } },
182 { { 0x63cf,0x6aee,0xbe6f,0x8837,0xc017, 0, 0, 0 } },
183 { { 0x26bb,0xccc7,0xb009,0xef75,0xc017, 0, 0, 0 } },
184 { { 0x462b,0xbae8,0xab96,0xa02a,0xc017, 0, 0, 0 } }
185};
186#endif
187#if MIEEE
188static uLD B[] = {
189 { { 0xc00a0000,0x873b0da1,0x49959557, 0 } },
190 { { 0xc00f0000,0xaa635b8c,0x9af8fe44, 0 } },
191 { { 0xc0130000,0x86ce3684,0x7cf55aa8, 0 } },
192 { { 0xc0150000,0xba7ff206,0x258c259a, 0 } },
193 { { 0xc0160000,0xf80ac0a0,0x1ca3be18, 0 } },
194 { { 0xc0170000,0x98e36717,0x2c42168f, 0 } },
195 { { 0xc0160000,0x876e92c8,0x9d552051, 0 } }
196};
197static uLD C[] = {
198 { { 0xc0080000,0x807cae76,0xcf2faa77, 0 } },
199 { { 0xc00e0000,0x84f3b55a,0x0d74b280, 0 } },
200 { { 0xc0120000,0x980981dc,0xcd30a505, 0 } },
201 { { 0xc0150000,0x92f0b8c2,0x42463369, 0 } },
202 { { 0xc0170000,0x8837be6f,0x6aee63cf, 0 } },
203 { { 0xc0170000,0xef75b009,0xccc726bb, 0 } },
204 { { 0xc0170000,0xa02aab96,0xbae8462b, 0 } }
205};
206#endif
207
208/* log( sqrt( 2*pi ) ) */
209static const long double LS2PI = 0.91893853320467274178L;
210#if defined(__arm__) || defined(_ARM_) || defined(__aarch64__) || defined(_ARM64_)
211#define MAXLGM 2.035093e36
212#else
213#define MAXLGM 1.04848146839019521116e+4928L
214#endif /* defined(__arm__) || defined(_ARM_) || defined(__aarch64__) || defined(_ARM64_) */
215
216/* Logarithm of gamma function */
217/* Reentrant version */
218long double __lgammal_r(long double x, int* sgngaml);
219
220long double __lgammal_r(long double x, int* sgngaml)
221{
222 long double p, q, w, z, f, nx;
223 int i;
224
225 *sgngaml = 1;
226#ifdef NANS
227 if (isnanl(x))
228 return x;
229#endif
230#ifdef INFINITIES
231 if (!isfinitel(x))
232 return (INFINITYL);
233#endif
234 if (x < -34.0L)
235 {
236 q = -x;
237 w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */
238 p = floorl(q);
239 if (p == q)
240 {
241lgsing:
242 _SET_ERRNO(EDOM);
243 mtherr( "lgammal", SING );
244#ifdef INFINITIES
245 return (INFINITYL);
246#else
247 return (MAXNUML);
248#endif
249 }
250 i = p;
251 if ((i & 1) == 0)
252 *sgngaml = -1;
253 else
254 *sgngaml = 1;
255 z = q - p;
256 if (z > 0.5L)
257 {
258 p += 1.0L;
259 z = p - q;
260 }
261 z = q * sinl(PIL * z);
262 if (z == 0.0L)
263 goto lgsing;
264 /* z = LOGPI - logl( z ) - w; */
265 z = logl(PIL/z) - w;
266 return (z);
267 }
268
269 if (x < 13.0L)
270 {
271 z = 1.0L;
272 nx = floorl(x + 0.5L);
273 f = x - nx;
274 while (x >= 3.0L)
275 {
276 nx -= 1.0L;
277 x = nx + f;
278 z *= x;
279 }
280 while (x < 2.0L)
281 {
282 if (fabsl(x) <= 0.03125)
283 goto lsmall;
284 z /= nx + f;
285 nx += 1.0L;
286 x = nx + f;
287 }
288 if (z < 0.0L)
289 {
290 *sgngaml = -1;
291 z = -z;
292 }
293 else
294 *sgngaml = 1;
295 if (x == 2.0L)
296 return ( logl(z) );
297 x = (nx - 2.0L) + f;
298 p = x * polevll(x, B, 6) / p1evll(x, C, 7);
299 return ( logl(z) + p );
300 }
301
302 if (x > MAXLGM)
303 {
304 _SET_ERRNO(ERANGE);
305 mtherr("lgammal", OVERFLOW);
306#ifdef INFINITIES
307 return (*sgngaml * INFINITYL);
308#else
309 return (*sgngaml * MAXNUML);
310#endif
311 }
312
313 q = (x - 0.5L) * logl(x) - x + LS2PI;
314 if (x > 1.0e10L)
315 return(q);
316 p = 1.0L/(x*x);
317 q += polevll(p, A, 6) / x;
318 return (q);
319
320lsmall:
321 if (x == 0.0L)
322 goto lgsing;
323 if (x < 0.0L)
324 {
325 x = -x;
326 q = z / (x * polevll(x, SN, 8));
327 }
328 else
329 q = z / (x * polevll(x, S, 8));
330 if (q < 0.0L)
331 {
332 *sgngaml = -1;
333 q = -q;
334 }
335 else
336 *sgngaml = 1;
337 q = logl(q);
338 return (q);
339}
340
341/* This is the C99 version */
342long double lgammal(long double x)
343{
344 return (__lgammal_r (x, &signgam));
345}
346#endif