1 /* lgammal expanding around zeros.
2 Copyright (C) 2015-2024 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
21 #include <math_private.h>
22 #include <fenv_private.h>
24 static const long double lgamma_zeros
[][2] =
26 { -0x2.74ff92c01f0d82acp
+0L, 0x1.360cea0e5f8ed3ccp
-68L },
27 { -0x2.bf6821437b201978p
+0L, -0x1.95a4b4641eaebf4cp
-64L },
28 { -0x3.24c1b793cb35efb8p
+0L, -0xb.e699ad3d9ba6545p
-68L },
29 { -0x3.f48e2a8f85fca17p
+0L, -0xd.4561291236cc321p
-68L },
30 { -0x4.0a139e16656030cp
+0L, -0x3.9f0b0de18112ac18p
-64L },
31 { -0x4.fdd5de9bbabf351p
+0L, -0xd.0aa4076988501d8p
-68L },
32 { -0x5.021a95fc2db64328p
+0L, -0x2.4c56e595394decc8p
-64L },
33 { -0x5.ffa4bd647d0357ep
+0L, 0x2.b129d342ce12071cp
-64L },
34 { -0x6.005ac9625f233b6p
+0L, -0x7.c2d96d16385cb868p
-68L },
35 { -0x6.fff2fddae1bbff4p
+0L, 0x2.9d949a3dc02de0cp
-64L },
36 { -0x7.000cff7b7f87adf8p
+0L, 0x3.b7d23246787d54d8p
-64L },
37 { -0x7.fffe5fe05673c3c8p
+0L, -0x2.9e82b522b0ca9d3p
-64L },
38 { -0x8.0001a01459fc9f6p
+0L, -0xc.b3cec1cec857667p
-68L },
39 { -0x8.ffffd1c425e81p
+0L, 0x3.79b16a8b6da6181cp
-64L },
40 { -0x9.00002e3bb47d86dp
+0L, -0x6.d843fedc351deb78p
-64L },
41 { -0x9.fffffb606bdfdcdp
+0L, -0x6.2ae77a50547c69dp
-68L },
42 { -0xa.0000049f93bb992p
+0L, -0x7.b45d95e15441e03p
-64L },
43 { -0xa.ffffff9466e9f1bp
+0L, -0x3.6dacd2adbd18d05cp
-64L },
44 { -0xb.0000006b9915316p
+0L, 0x2.69a590015bf1b414p
-64L },
45 { -0xb.fffffff70893874p
+0L, 0x7.821be533c2c36878p
-64L },
46 { -0xc.00000008f76c773p
+0L, -0x1.567c0f0250f38792p
-64L },
47 { -0xc.ffffffff4f6dcf6p
+0L, -0x1.7f97a5ffc757d548p
-64L },
48 { -0xd.00000000b09230ap
+0L, 0x3.f997c22e46fc1c9p
-64L },
49 { -0xd.fffffffff36345bp
+0L, 0x4.61e7b5c1f62ee89p
-64L },
50 { -0xe.000000000c9cba5p
+0L, -0x4.5e94e75ec5718f78p
-64L },
51 { -0xe.ffffffffff28c06p
+0L, -0xc.6604ef30371f89dp
-68L },
52 { -0xf.0000000000d73fap
+0L, 0xc.6642f1bdf07a161p
-68L },
53 { -0xf.fffffffffff28cp
+0L, -0x6.0c6621f512e72e5p
-64L },
54 { -0x1.000000000000d74p
+4L, 0x6.0c6625ebdb406c48p
-64L },
55 { -0x1.0ffffffffffff356p
+4L, -0x9.c47e7a93e1c46a1p
-64L },
56 { -0x1.1000000000000caap
+4L, 0x9.c47e7a97778935ap
-64L },
57 { -0x1.1fffffffffffff4cp
+4L, 0x1.3c31dcbecd2f74d4p
-64L },
58 { -0x1.20000000000000b4p
+4L, -0x1.3c31dcbeca4c3b3p
-64L },
59 { -0x1.2ffffffffffffff6p
+4L, -0x8.5b25cbf5f545ceep
-64L },
60 { -0x1.300000000000000ap
+4L, 0x8.5b25cbf5f547e48p
-64L },
61 { -0x1.4p
+4L, 0x7.950ae90080894298p
-64L },
62 { -0x1.4p
+4L, -0x7.950ae9008089414p
-64L },
63 { -0x1.5p
+4L, 0x5.c6e3bdb73d5c63p
-68L },
64 { -0x1.5p
+4L, -0x5.c6e3bdb73d5c62f8p
-68L },
65 { -0x1.6p
+4L, 0x4.338e5b6dfe14a518p
-72L },
66 { -0x1.6p
+4L, -0x4.338e5b6dfe14a51p
-72L },
67 { -0x1.7p
+4L, 0x2.ec368262c7033b3p
-76L },
68 { -0x1.7p
+4L, -0x2.ec368262c7033b3p
-76L },
69 { -0x1.8p
+4L, 0x1.f2cf01972f577ccap
-80L },
70 { -0x1.8p
+4L, -0x1.f2cf01972f577ccap
-80L },
71 { -0x1.9p
+4L, 0x1.3f3ccdd165fa8d4ep
-84L },
72 { -0x1.9p
+4L, -0x1.3f3ccdd165fa8d4ep
-84L },
73 { -0x1.ap
+4L, 0xc.4742fe35272cd1cp
-92L },
74 { -0x1.ap
+4L, -0xc.4742fe35272cd1cp
-92L },
75 { -0x1.bp
+4L, 0x7.46ac70b733a8c828p
-96L },
76 { -0x1.bp
+4L, -0x7.46ac70b733a8c828p
-96L },
77 { -0x1.cp
+4L, 0x4.2862898d42174ddp
-100L },
78 { -0x1.cp
+4L, -0x4.2862898d42174ddp
-100L },
79 { -0x1.dp
+4L, 0x2.4b3f31686b15af58p
-104L },
80 { -0x1.dp
+4L, -0x2.4b3f31686b15af58p
-104L },
81 { -0x1.ep
+4L, 0x1.3932c5047d60e60cp
-108L },
82 { -0x1.ep
+4L, -0x1.3932c5047d60e60cp
-108L },
83 { -0x1.fp
+4L, 0xa.1a6973c1fade217p
-116L },
84 { -0x1.fp
+4L, -0xa.1a6973c1fade217p
-116L },
85 { -0x2p
+4L, 0x5.0d34b9e0fd6f10b8p
-120L },
86 { -0x2p
+4L, -0x5.0d34b9e0fd6f10b8p
-120L },
87 { -0x2.1p
+4L, 0x2.73024a9ba1aa36a8p
-124L },
90 static const long double e_hi
= 0x2.b7e151628aed2a6cp
+0L;
91 static const long double e_lo
= -0x1.408ea77f630b0c38p
-64L;
93 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
94 approximation to lgamma function. */
96 static const long double lgamma_coeff
[] =
98 0x1.5555555555555556p
-4L,
99 -0xb.60b60b60b60b60bp
-12L,
100 0x3.4034034034034034p
-12L,
101 -0x2.7027027027027028p
-12L,
102 0x3.72a3c5631fe46aep
-12L,
103 -0x7.daac36664f1f208p
-12L,
104 0x1.a41a41a41a41a41ap
-8L,
105 -0x7.90a1b2c3d4e5f708p
-8L,
106 0x2.dfd2c703c0cfff44p
-4L,
107 -0x1.6476701181f39edcp
+0L,
108 0xd.672219167002d3ap
+0L,
109 -0x9.cd9292e6660d55bp
+4L,
110 0x8.911a740da740da7p
+8L,
111 -0x8.d0cc570e255bf5ap
+12L,
112 0xa.8d1044d3708d1c2p
+16L,
113 -0xe.8844d8a169abbc4p
+20L,
116 #define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
118 /* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
119 the integer end-point of the half-integer interval containing x and
120 x0 is the zero of lgamma in that half-integer interval. Each
121 polynomial is expressed in terms of x-xm, where xm is the midpoint
122 of the interval for which the polynomial applies. */
124 static const long double poly_coeff
[] =
126 /* Interval [-2.125, -2] (polynomial degree 13). */
127 -0x1.0b71c5c54d42eb6cp
+0L,
128 -0xc.73a1dc05f349517p
-4L,
129 -0x1.ec841408528b6baep
-4L,
130 -0xe.37c9da26fc3b492p
-4L,
131 -0x1.03cd87c5178991ap
-4L,
132 -0xe.ae9ada65ece2f39p
-4L,
133 0x9.b1185505edac18dp
-8L,
134 -0xe.f28c130b54d3cb2p
-4L,
135 0x2.6ec1666cf44a63bp
-4L,
136 -0xf.57cb2774193bbd5p
-4L,
137 0x4.5ae64671a41b1c4p
-4L,
138 -0xf.f48ea8b5bd3a7cep
-4L,
139 0x6.7d73788a8d30ef58p
-4L,
140 -0x1.11e0e4b506bd272ep
+0L,
141 /* Interval [-2.25, -2.125] (polynomial degree 13). */
142 -0xf.2930890d7d675a8p
-4L,
143 -0xc.a5cfde054eab5cdp
-4L,
144 0x3.9c9e0fdebb0676e4p
-4L,
145 -0x1.02a5ad35605f0d8cp
+0L,
146 0x9.6e9b1185d0b92edp
-4L,
147 -0x1.4d8332f3d6a3959p
+0L,
148 0x1.1c0c8cacd0ced3eap
+0L,
149 -0x1.c9a6f592a67b1628p
+0L,
150 0x1.d7e9476f96aa4bd6p
+0L,
151 -0x2.921cedb488bb3318p
+0L,
152 0x2.e8b3fd6ca193e4c8p
+0L,
153 -0x3.cb69d9d6628e4a2p
+0L,
154 0x4.95f12c73b558638p
+0L,
155 -0x5.d392d0b97c02ab6p
+0L,
156 /* Interval [-2.375, -2.25] (polynomial degree 14). */
157 -0xd.7d28d505d618122p
-4L,
158 -0xe.69649a304098532p
-4L,
159 0xb.0d74a2827d055c5p
-4L,
160 -0x1.924b09228531c00ep
+0L,
161 0x1.d49b12bccee4f888p
+0L,
162 -0x3.0898bb7dbb21e458p
+0L,
163 0x4.207a6cad6fa10a2p
+0L,
164 -0x6.39ee630b46093ad8p
+0L,
165 0x8.e2e25211a3fb5ccp
+0L,
166 -0xd.0e85ccd8e79c08p
+0L,
167 0x1.2e45882bc17f9e16p
+4L,
168 -0x1.b8b6e841815ff314p
+4L,
169 0x2.7ff8bf7504fa04dcp
+4L,
170 -0x3.c192e9c903352974p
+4L,
171 0x5.8040b75f4ef07f98p
+4L,
172 /* Interval [-2.5, -2.375] (polynomial degree 15). */
173 -0xb.74ea1bcfff94b2cp
-4L,
174 -0x1.2a82bd590c375384p
+0L,
175 0x1.88020f828b968634p
+0L,
176 -0x3.32279f040eb80fa4p
+0L,
177 0x5.57ac825175943188p
+0L,
178 -0x9.c2aedcfe10f129ep
+0L,
179 0x1.12c132f2df02881ep
+4L,
180 -0x1.ea94e26c0b6ffa6p
+4L,
181 0x3.66b4a8bb0290013p
+4L,
182 -0x6.0cf735e01f5990bp
+4L,
183 0xa.c10a8db7ae99343p
+4L,
184 -0x1.31edb212b315feeap
+8L,
185 0x2.1f478592298b3ebp
+8L,
186 -0x3.c546da5957ace6ccp
+8L,
187 0x7.0e3d2a02579ba4bp
+8L,
188 -0xc.b1ea961c39302f8p
+8L,
189 /* Interval [-2.625, -2.5] (polynomial degree 16). */
190 -0x3.d10108c27ebafad4p
-4L,
191 0x1.cd557caff7d2b202p
+0L,
192 0x3.819b4856d3995034p
+0L,
193 0x6.8505cbad03dd3bd8p
+0L,
194 0xb.c1b2e653aa0b924p
+0L,
195 0x1.50a53a38f05f72d6p
+4L,
196 0x2.57ae00cbd06efb34p
+4L,
197 0x4.2b1563077a577e9p
+4L,
198 0x7.6989ed790138a7f8p
+4L,
199 0xd.2dd28417b4f8406p
+4L,
200 0x1.76e1b71f0710803ap
+8L,
201 0x2.9a7a096254ac032p
+8L,
202 0x4.a0e6109e2a039788p
+8L,
203 0x8.37ea17a93c877b2p
+8L,
204 0xe.9506a641143612bp
+8L,
205 0x1.b680ed4ea386d52p
+12L,
206 0x3.28a2130c8de0ae84p
+12L,
207 /* Interval [-2.75, -2.625] (polynomial degree 15). */
208 -0x6.b5d252a56e8a7548p
-4L,
209 0x1.28d60383da3ac72p
+0L,
210 0x1.db6513ada8a6703ap
+0L,
211 0x2.e217118f9d34aa7cp
+0L,
212 0x4.450112c5cbd6256p
+0L,
213 0x6.4af99151e972f92p
+0L,
214 0x9.2db598b5b183cd6p
+0L,
215 0xd.62bef9c9adcff6ap
+0L,
216 0x1.379f290d743d9774p
+4L,
217 0x1.c58271ff823caa26p
+4L,
218 0x2.93a871b87a06e73p
+4L,
219 0x3.bf9db66103d7ec98p
+4L,
220 0x5.73247c111fbf197p
+4L,
221 0x7.ec8b9973ba27d008p
+4L,
222 0xb.eca5f9619b39c03p
+4L,
223 0x1.18f2e46411c78b1cp
+8L,
224 /* Interval [-2.875, -2.75] (polynomial degree 14). */
225 -0x8.a41b1e4f36ff88ep
-4L,
226 0xc.da87d3b69dc0f34p
-4L,
227 0x1.1474ad5c36158ad2p
+0L,
228 0x1.761ecb90c5553996p
+0L,
229 0x1.d279bff9ae234f8p
+0L,
230 0x2.4e5d0055a16c5414p
+0L,
231 0x2.d57545a783902f8cp
+0L,
232 0x3.8514eec263aa9f98p
+0L,
233 0x4.5235e338245f6fe8p
+0L,
234 0x5.562b1ef200b256c8p
+0L,
235 0x6.8ec9782b93bd565p
+0L,
236 0x8.14baf4836483508p
+0L,
237 0x9.efaf35dc712ea79p
+0L,
238 0xc.8431f6a226507a9p
+0L,
239 0xf.80358289a768401p
+0L,
240 /* Interval [-3, -2.875] (polynomial degree 13). */
241 -0xa.046d667e468f3e4p
-4L,
242 0x9.70b88dcc006c216p
-4L,
243 0xa.a8a39421c86ce9p
-4L,
244 0xd.2f4d1363f321e89p
-4L,
245 0xd.ca9aa1a3ab2f438p
-4L,
246 0xf.cf09c31f05d02cbp
-4L,
247 0x1.04b133a195686a38p
+0L,
248 0x1.22b54799d0072024p
+0L,
249 0x1.2c5802b869a36ae8p
+0L,
250 0x1.4aadf23055d7105ep
+0L,
251 0x1.5794078dd45c55d6p
+0L,
252 0x1.7759069da18bcf0ap
+0L,
253 0x1.8e672cefa4623f34p
+0L,
254 0x1.b2acfa32c17145e6p
+0L,
257 static const size_t poly_deg
[] =
269 static const size_t poly_end
[] =
281 /* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
284 lg_sinpi (long double x
)
287 return __sinl (M_PIl
* x
);
289 return __cosl (M_PIl
* (0.5L - x
));
292 /* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
295 lg_cospi (long double x
)
298 return __cosl (M_PIl
* x
);
300 return __sinl (M_PIl
* (0.5L - x
));
303 /* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
306 lg_cotpi (long double x
)
308 return lg_cospi (x
) / lg_sinpi (x
);
311 /* Compute lgamma of a negative argument -33 < X < -2, setting
312 *SIGNGAMP accordingly. */
315 __lgamma_negl (long double x
, int *signgamp
)
317 /* Determine the half-integer region X lies in, handle exact
318 integers and determine the sign of the result. */
319 int i
= floorl (-2 * x
);
320 if ((i
& 1) == 0 && i
== -2 * x
)
322 long double xn
= ((i
& 1) == 0 ? -i
/ 2 : (-i
- 1) / 2);
324 *signgamp
= ((i
& 2) == 0 ? -1 : 1);
326 SET_RESTORE_ROUNDL (FE_TONEAREST
);
328 /* Expand around the zero X0 = X0_HI + X0_LO. */
329 long double x0_hi
= lgamma_zeros
[i
][0], x0_lo
= lgamma_zeros
[i
][1];
330 long double xdiff
= x
- x0_hi
- x0_lo
;
332 /* For arguments in the range -3 to -2, use polynomial
333 approximations to an adjusted version of the gamma function. */
336 int j
= floorl (-8 * x
) - 16;
337 long double xm
= (-33 - 2 * j
) * 0.0625L;
338 long double x_adj
= x
- xm
;
339 size_t deg
= poly_deg
[j
];
340 size_t end
= poly_end
[j
];
341 long double g
= poly_coeff
[end
];
342 for (size_t j
= 1; j
<= deg
; j
++)
343 g
= g
* x_adj
+ poly_coeff
[end
- j
];
344 return __log1pl (g
* xdiff
/ (x
- xn
));
347 /* The result we want is log (sinpi (X0) / sinpi (X))
348 + log (gamma (1 - X0) / gamma (1 - X)). */
349 long double x_idiff
= fabsl (xn
- x
), x0_idiff
= fabsl (xn
- x0_hi
- x0_lo
);
350 long double log_sinpi_ratio
;
351 if (x0_idiff
< x_idiff
* 0.5L)
352 /* Use log not log1p to avoid inaccuracy from log1p of arguments
354 log_sinpi_ratio
= __ieee754_logl (lg_sinpi (x0_idiff
)
355 / lg_sinpi (x_idiff
));
358 /* Use log1p not log to avoid inaccuracy from log of arguments
359 close to 1. X0DIFF2 has positive sign if X0 is further from
360 XN than X is from XN, negative sign otherwise. */
361 long double x0diff2
= ((i
& 1) == 0 ? xdiff
: -xdiff
) * 0.5L;
362 long double sx0d2
= lg_sinpi (x0diff2
);
363 long double cx0d2
= lg_cospi (x0diff2
);
364 log_sinpi_ratio
= __log1pl (2 * sx0d2
365 * (-sx0d2
+ cx0d2
* lg_cotpi (x_idiff
)));
368 long double log_gamma_ratio
;
369 long double y0
= 1 - x0_hi
;
370 long double y0_eps
= -x0_hi
+ (1 - y0
) - x0_lo
;
371 long double y
= 1 - x
;
372 long double y_eps
= -x
+ (1 - y
);
373 /* We now wish to compute LOG_GAMMA_RATIO
374 = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
375 accurately approximates the difference Y0 + Y0_EPS - Y -
376 Y_EPS. Use Stirling's approximation. First, we may need to
377 adjust into the range where Stirling's approximation is
378 sufficiently accurate. */
379 long double log_gamma_adj
= 0;
382 int n_up
= (9 - i
) / 2;
383 long double ny0
, ny0_eps
, ny
, ny_eps
;
385 ny0_eps
= y0
- (ny0
- n_up
) + y0_eps
;
389 ny_eps
= y
- (ny
- n_up
) + y_eps
;
392 long double prodm1
= __lgamma_productl (xdiff
, y
- n_up
, y_eps
, n_up
);
393 log_gamma_adj
= -__log1pl (prodm1
);
395 long double log_gamma_high
396 = (xdiff
* __log1pl ((y0
- e_hi
- e_lo
+ y0_eps
) / e_hi
)
397 + (y
- 0.5L + y_eps
) * __log1pl (xdiff
/ y
) + log_gamma_adj
);
398 /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
399 long double y0r
= 1 / y0
, yr
= 1 / y
;
400 long double y0r2
= y0r
* y0r
, yr2
= yr
* yr
;
401 long double rdiff
= -xdiff
/ (y
* y0
);
402 long double bterm
[NCOEFF
];
403 long double dlast
= rdiff
, elast
= rdiff
* yr
* (yr
+ y0r
);
404 bterm
[0] = dlast
* lgamma_coeff
[0];
405 for (size_t j
= 1; j
< NCOEFF
; j
++)
407 long double dnext
= dlast
* y0r2
+ elast
;
408 long double enext
= elast
* yr2
;
409 bterm
[j
] = dnext
* lgamma_coeff
[j
];
413 long double log_gamma_low
= 0;
414 for (size_t j
= 0; j
< NCOEFF
; j
++)
415 log_gamma_low
+= bterm
[NCOEFF
- 1 - j
];
416 log_gamma_ratio
= log_gamma_high
+ log_gamma_low
;
418 return log_sinpi_ratio
+ log_gamma_ratio
;