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/* A[]: Stirling's formula expansion of log gamma
9 * B[], C[]: log gamma function between 2 and 3
10 */
11#ifdef UNK
12static uD A[] = {
13 { { 8.11614167470508450300E-4 } },
14 { { -5.95061904284301438324E-4 } },
15 { { 7.93650340457716943945E-4 } },
16 { { -2.77777777730099687205E-3 } },
17 { { 8.33333333333331927722E-2 } }
18};
19static uD B[] = {
20 { { -1.37825152569120859100E3 } },
21 { { -3.88016315134637840924E4 } },
22 { { -3.31612992738871184744E5 } },
23 { { -1.16237097492762307383E6 } },
24 { { -1.72173700820839662146E6 } },
25 { { -8.53555664245765465627E5 } }
26};
27static uD C[] = {
28 { { -3.51815701436523470549E2 } },
29 { { -1.70642106651881159223E4 } },
30 { { -2.20528590553854454839E5 } },
31 { { -1.13933444367982507207E6 } },
32 { { -2.53252307177582951285E6 } },
33 { { -2.01889141433532773231E6 } }
34};
35/* log( sqrt( 2*pi ) ) */
36static double LS2PI = 0.91893853320467274178;
37#define MAXLGM 2.556348e305
38static double LOGPI = 1.14472988584940017414;
39#endif
40
41#ifdef DEC
42static const uD A[] = {
43 { { 0035524,0141201,0034633,0031405 } },
44 { { 0135433,0176755,0126007,0045030 } },
45 { { 0035520,0006371,0003342,0172730 } },
46 { { 0136066,0005540,0132605,0026407 } },
47 { { 0037252,0125252,0125252,0125132 } }
48};
49static const uD B[] = {
50 { { 0142654,0044014,0077633,0035410 } },
51 { { 0144027,0110641,0125335,0144760 } },
52 { { 0144641,0165637,0142204,0047447 } },
53 { { 0145215,0162027,0146246,0155211 } },
54 { { 0145322,0026110,0010317,0110130 } },
55 { { 0145120,0061472,0120300,0025363 } }
56};
57static const uD C[] = {
58 { { 0142257,0164150,0163630,0112622 } },
59 { { 0143605,0050153,0156116,0135272 } },
60 { { 0144527,0056045,0145642,0062332 } },
61 { { 0145213,0012063,0106250,0001025 } },
62 { { 0145432,0111254,0044577,0115142 } },
63 { { 0145366,0071133,0050217,0005122 } }
64};
65/* log( sqrt( 2*pi ) ) */
66static const uD LS2P[] = { {040153,037616,041445,0172645,} };
67#define LS2PI LS2P[0].d
68#define MAXLGM 2.035093e36
69static const uD LPI[] = { { 0040222,0103202,0043475,0006750, } };
70#define LOGPI LPI[0].d
71
72#endif
73
74#ifdef IBMPC
75static const uD A[] = {
76 { { 0x6661,0x2733,0x9850,0x3f4a } },
77 { { 0xe943,0xb580,0x7fbd,0xbf43 } },
78 { { 0x5ebb,0x20dc,0x019f,0x3f4a } },
79 { { 0xa5a1,0x16b0,0xc16c,0xbf66 } },
80 { { 0x554b,0x5555,0x5555,0x3fb5 } }
81};
82static const uD B[] = {
83 { { 0x6761,0x8ff3,0x8901,0xc095 } },
84 { { 0xb93e,0x355b,0xf234,0xc0e2 } },
85 { { 0x89e5,0xf890,0x3d73,0xc114 } },
86 { { 0xdb51,0xf994,0xbc82,0xc131 } },
87 { { 0xf20b,0x0219,0x4589,0xc13a } },
88 { { 0x055e,0x5418,0x0c67,0xc12a } }
89};
90static const uD C[] = {
91 { { 0x12b2,0x1cf3,0xfd0d,0xc075 } },
92 { { 0xd757,0x7b89,0xaa0d,0xc0d0 } },
93 { { 0x4c9b,0xb974,0xeb84,0xc10a } },
94 { { 0x0043,0x7195,0x6286,0xc131 } },
95 { { 0xf34c,0x892f,0x5255,0xc143 } },
96 { { 0xe14a,0x6a11,0xce4b,0xc13e } }
97};
98/* log( sqrt( 2*pi ) ) */
99static const union
100{
101 unsigned short s[4];
102 double d;
103} ls2p = {{0xbeb5,0xc864,0x67f1,0x3fed}};
104#define LS2PI (ls2p.d)
105#define MAXLGM 2.556348e305
106/* log (pi) */
107static const union
108{
109 unsigned short s[4];
110 double d;
111} lpi = {{0xa1bd,0x48e7,0x50d0,0x3ff2}};
112#define LOGPI (lpi.d)
113#endif
114
115#ifdef MIEEE
116static const uD A[] = {
117 { { 0x3f4a,0x9850,0x2733,0x6661 } },
118 { { 0xbf43,0x7fbd,0xb580,0xe943 } },
119 { { 0x3f4a,0x019f,0x20dc,0x5ebb } },
120 { { 0xbf66,0xc16c,0x16b0,0xa5a1 } },
121 { { 0x3fb5,0x5555,0x5555,0x554b } }
122};
123static const uD B[] = {
124 { { 0xc095,0x8901,0x8ff3,0x6761 } },
125 { { 0xc0e2,0xf234,0x355b,0xb93e } },
126 { { 0xc114,0x3d73,0xf890,0x89e5 } },
127 { { 0xc131,0xbc82,0xf994,0xdb51 } },
128 { { 0xc13a,0x4589,0x0219,0xf20b } },
129 { { 0xc12a,0x0c67,0x5418,0x055e } }
130};
131static const uD C[] = {
132 { { 0xc075,0xfd0d,0x1cf3,0x12b2 } },
133 { { 0xc0d0,0xaa0d,0x7b89,0xd757 } },
134 { { 0xc10a,0xeb84,0xb974,0x4c9b } },
135 { { 0xc131,0x6286,0x7195,0x0043 } },
136 { { 0xc143,0x5255,0x892f,0xf34c } },
137 { { 0xc13e,0xce4b,0x6a11,0xe14a } }
138};
139/* log( sqrt( 2*pi ) ) */
140static const union
141{
142 unsigned short s[4];
143 double d;
144} ls2p = {{0x3fed,0x67f1,0xc864,0xbeb5}};
145#define LS2PI ls2p.d
146#define MAXLGM 2.556348e305
147/* log (pi) */
148static const union
149{
150 unsigned short s[4];
151 double d;
152} lpi = {{0x3ff2, 0x50d0, 0x48e7, 0xa1bd}};
153#define LOGPI (lpi.d)
154#endif
155
156
157/* Logarithm of gamma function */
158/* Reentrant version */
159double __lgamma_r(double x, int* sgngam);
160
161double __lgamma_r(double x, int* sgngam)
162{
163 double p, q, u, w, z;
164 int i;
165
166 *sgngam = 1;
167#ifdef NANS
168 if (isnan(x))
169 return (x);
170#endif
171
172#ifdef INFINITIES
173 if (!isfinite(x))
174 return (INFINITY);
175#endif
176
177 if (x < -34.0)
178 {
179 q = -x;
180 w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */
181 p = floor(q);
182 if (p == q)
183 {
184lgsing:
185 _SET_ERRNO(EDOM);
186 mtherr( "lgam", SING );
187#ifdef INFINITIES
188 return (INFINITY);
189#else
190 return (MAXNUM);
191#endif
192 }
193 i = p;
194 if ((i & 1) == 0)
195 *sgngam = -1;
196 else
197 *sgngam = 1;
198 z = q - p;
199 if (z > 0.5)
200 {
201 p += 1.0;
202 z = p - q;
203 }
204 z = q * sin( PI * z );
205 if (z == 0.0)
206 goto lgsing;
207 /* z = log(PI) - log( z ) - w;*/
208 z = LOGPI - log( z ) - w;
209 return (z);
210 }
211
212 if (x < 13.0)
213 {
214 z = 1.0;
215 p = 0.0;
216 u = x;
217 while (u >= 3.0)
218 {
219 p -= 1.0;
220 u = x + p;
221 z *= u;
222 }
223 while (u < 2.0)
224 {
225 if (u == 0.0)
226 goto lgsing;
227 z /= u;
228 p += 1.0;
229 u = x + p;
230 }
231 if (z < 0.0)
232 {
233 *sgngam = -1;
234 z = -z;
235 }
236 else
237 *sgngam = 1;
238 if (u == 2.0)
239 return ( log(z) );
240 p -= 2.0;
241 x = x + p;
242 p = x * polevl(x, B, 5) / p1evl(x, C, 6);
243 return ( log(z) + p );
244 }
245
246 if (x > MAXLGM)
247 {
248 _SET_ERRNO(ERANGE);
249 mtherr("lgamma", OVERFLOW);
250#ifdef INFINITIES
251 return (*sgngam * INFINITY);
252#else
253 return (*sgngam * MAXNUM);
254#endif
255 }
256
257 q = (x - 0.5) * log(x) - x + LS2PI;
258 if (x > 1.0e8)
259 return (q);
260
261 p = 1.0/(x*x);
262 if (x >= 1000.0)
263 q += (( 7.9365079365079365079365e-4 * p
264 - 2.7777777777777777777778e-3) *p
265 + 0.0833333333333333333333) / x;
266 else
267 q += polevl( p, A, 4 ) / x;
268 return (q);
269}
270
271/* This is the C99 version */
272double lgamma(double x)
273{
274 return (__lgamma_r(x, &signgam));
275}
276