1 /* lgamma expanding around zeros.
2 Copyright (C) 2015-2023 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-narrow-eval.h>
22 #include <math_private.h>
23 #include <fenv_private.h>
25 static const double lgamma_zeros
[][2] =
27 { -0x2.74ff92c01f0d8p
+0, -0x2.abec9f315f1ap
-56 },
28 { -0x2.bf6821437b202p
+0, 0x6.866a5b4b9be14p
-56 },
29 { -0x3.24c1b793cb35ep
+0, -0xf.b8be699ad3d98p
-56 },
30 { -0x3.f48e2a8f85fcap
+0, -0x1.70d4561291237p
-56 },
31 { -0x4.0a139e1665604p
+0, 0xf.3c60f4f21e7fp
-56 },
32 { -0x4.fdd5de9bbabf4p
+0, 0xa.ef2f55bf89678p
-56 },
33 { -0x5.021a95fc2db64p
+0, -0x3.2a4c56e595394p
-56 },
34 { -0x5.ffa4bd647d034p
+0, -0x1.7dd4ed62cbd32p
-52 },
35 { -0x6.005ac9625f234p
+0, 0x4.9f83d2692e9c8p
-56 },
36 { -0x6.fff2fddae1bcp
+0, 0xc.29d949a3dc03p
-60 },
37 { -0x7.000cff7b7f87cp
+0, 0x1.20bb7d2324678p
-52 },
38 { -0x7.fffe5fe05673cp
+0, -0x3.ca9e82b522b0cp
-56 },
39 { -0x8.0001a01459fc8p
+0, -0x1.f60cb3cec1cedp
-52 },
40 { -0x8.ffffd1c425e8p
+0, -0xf.fc864e9574928p
-56 },
41 { -0x9.00002e3bb47d8p
+0, -0x6.d6d843fedc35p
-56 },
42 { -0x9.fffffb606bep
+0, 0x2.32f9d51885afap
-52 },
43 { -0xa.0000049f93bb8p
+0, -0x1.927b45d95e154p
-52 },
44 { -0xa.ffffff9466eap
+0, 0xe.4c92532d5243p
-56 },
45 { -0xb.0000006b9915p
+0, -0x3.15d965a6ffea4p
-52 },
46 { -0xb.fffffff708938p
+0, -0x7.387de41acc3d4p
-56 },
47 { -0xc.00000008f76c8p
+0, 0x8.cea983f0fdafp
-56 },
48 { -0xc.ffffffff4f6ep
+0, 0x3.09e80685a0038p
-52 },
49 { -0xd.00000000b092p
+0, -0x3.09c06683dd1bap
-52 },
50 { -0xd.fffffffff3638p
+0, 0x3.a5461e7b5c1f6p
-52 },
51 { -0xe.000000000c9c8p
+0, -0x3.a545e94e75ec6p
-52 },
52 { -0xe.ffffffffff29p
+0, 0x3.f9f399fb10cfcp
-52 },
53 { -0xf.0000000000d7p
+0, -0x3.f9f399bd0e42p
-52 },
54 { -0xf.fffffffffff28p
+0, -0xc.060c6621f513p
-56 },
55 { -0x1.000000000000dp
+4, -0x7.3f9f399da1424p
-52 },
56 { -0x1.0ffffffffffffp
+4, -0x3.569c47e7a93e2p
-52 },
57 { -0x1.1000000000001p
+4, 0x3.569c47e7a9778p
-52 },
58 { -0x1.2p
+4, 0xb.413c31dcbecdp
-56 },
59 { -0x1.2p
+4, -0xb.413c31dcbeca8p
-56 },
60 { -0x1.3p
+4, 0x9.7a4da340a0ab8p
-60 },
61 { -0x1.3p
+4, -0x9.7a4da340a0ab8p
-60 },
62 { -0x1.4p
+4, 0x7.950ae90080894p
-64 },
63 { -0x1.4p
+4, -0x7.950ae90080894p
-64 },
64 { -0x1.5p
+4, 0x5.c6e3bdb73d5c8p
-68 },
65 { -0x1.5p
+4, -0x5.c6e3bdb73d5c8p
-68 },
66 { -0x1.6p
+4, 0x4.338e5b6dfe14cp
-72 },
67 { -0x1.6p
+4, -0x4.338e5b6dfe14cp
-72 },
68 { -0x1.7p
+4, 0x2.ec368262c7034p
-76 },
69 { -0x1.7p
+4, -0x2.ec368262c7034p
-76 },
70 { -0x1.8p
+4, 0x1.f2cf01972f578p
-80 },
71 { -0x1.8p
+4, -0x1.f2cf01972f578p
-80 },
72 { -0x1.9p
+4, 0x1.3f3ccdd165fa9p
-84 },
73 { -0x1.9p
+4, -0x1.3f3ccdd165fa9p
-84 },
74 { -0x1.ap
+4, 0xc.4742fe35272dp
-92 },
75 { -0x1.ap
+4, -0xc.4742fe35272dp
-92 },
76 { -0x1.bp
+4, 0x7.46ac70b733a8cp
-96 },
77 { -0x1.bp
+4, -0x7.46ac70b733a8cp
-96 },
78 { -0x1.cp
+4, 0x4.2862898d42174p
-100 },
81 static const double e_hi
= 0x2.b7e151628aed2p
+0, e_lo
= 0xa.6abf7158809dp
-56;
83 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
84 approximation to lgamma function. */
86 static const double lgamma_coeff
[] =
89 -0xb.60b60b60b60b8p
-12,
90 0x3.4034034034034p
-12,
91 -0x2.7027027027028p
-12,
92 0x3.72a3c5631fe46p
-12,
93 -0x7.daac36664f1f4p
-12,
97 -0x1.6476701181f3ap
+0,
99 -0x9.cd9292e6660d8p
+4,
102 #define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
104 /* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
105 the integer end-point of the half-integer interval containing x and
106 x0 is the zero of lgamma in that half-integer interval. Each
107 polynomial is expressed in terms of x-xm, where xm is the midpoint
108 of the interval for which the polynomial applies. */
110 static const double poly_coeff
[] =
112 /* Interval [-2.125, -2] (polynomial degree 10). */
113 -0x1.0b71c5c54d42fp
+0,
114 -0xc.73a1dc05f3758p
-4,
115 -0x1.ec84140851911p
-4,
116 -0xe.37c9da23847e8p
-4,
117 -0x1.03cd87cdc0ac6p
-4,
118 -0xe.ae9aedce12eep
-4,
119 0x9.b11a1780cfd48p
-8,
120 -0xe.f25fc460bdebp
-4,
121 0x2.6e984c61ca912p
-4,
123 0x4.760c8c8909758p
-4,
124 /* Interval [-2.25, -2.125] (polynomial degree 11). */
125 -0xf.2930890d7d678p
-4,
126 -0xc.a5cfde054eaa8p
-4,
127 0x3.9c9e0fdebd99cp
-4,
128 -0x1.02a5ad35619d9p
+0,
130 -0x1.4d8332eba090ap
+0,
132 -0x1.c9a70d138c74ep
+0,
133 0x1.d7d9cf1d4c196p
+0,
134 -0x2.91fbf4cd6abacp
+0,
135 0x2.f6751f74b8ff8p
+0,
136 -0x3.e1bb7b09e3e76p
+0,
137 /* Interval [-2.375, -2.25] (polynomial degree 12). */
138 -0xd.7d28d505d618p
-4,
139 -0xe.69649a3040958p
-4,
141 -0x1.924b09228a86ep
+0,
142 0x1.d49b12bcf6175p
+0,
143 -0x3.0898bb530d314p
+0,
144 0x4.207a6be8fda4cp
+0,
146 0x8.e2e42acbccec8p
+0,
147 -0xd.0d91c1e596a68p
+0,
148 0x1.2e20d7099c585p
+4,
149 -0x1.c4eb6691b4ca9p
+4,
150 0x2.96a1a11fd85fep
+4,
151 /* Interval [-2.5, -2.375] (polynomial degree 13). */
152 -0xb.74ea1bcfff948p
-4,
153 -0x1.2a82bd590c376p
+0,
155 -0x3.32279f040d7aep
+0,
156 0x5.57ac8252ce868p
+0,
157 -0x9.c2aedd093125p
+0,
158 0x1.12c132716e94cp
+4,
159 -0x1.ea94dfa5c0a6dp
+4,
160 0x3.66b61abfe858cp
+4,
161 -0x6.0cfceb62a26e4p
+4,
162 0xa.beeba09403bd8p
+4,
163 -0x1.3188d9b1b288cp
+8,
164 0x2.37f774dd14c44p
+8,
165 -0x3.fdf0a64cd7136p
+8,
166 /* Interval [-2.625, -2.5] (polynomial degree 13). */
167 -0x3.d10108c27ebbp
-4,
168 0x1.cd557caff7d2fp
+0,
169 0x3.819b4856d36cep
+0,
172 0x1.50a53a3ce6c73p
+4,
173 0x2.57adffbb1ec0cp
+4,
174 0x4.2b15549cf400cp
+4,
175 0x7.698cfd82b3e18p
+4,
177 0x1.7699a624d07b9p
+8,
178 0x2.98ecf617abbfcp
+8,
179 0x4.d5244d44d60b4p
+8,
180 0x8.e962bf7395988p
+8,
181 /* Interval [-2.75, -2.625] (polynomial degree 12). */
182 -0x6.b5d252a56e8a8p
-4,
183 0x1.28d60383da3a6p
+0,
184 0x1.db6513ada89bep
+0,
185 0x2.e217118fa8c02p
+0,
186 0x4.450112c651348p
+0,
187 0x6.4af990f589b8cp
+0,
188 0x9.2db5963d7a238p
+0,
190 0x1.379f81f6416afp
+4,
192 0x2.9342d0af2ac4ep
+4,
193 0x3.d9cdf56d2b186p
+4,
194 0x5.ab9f91d5a27a4p
+4,
195 /* Interval [-2.875, -2.75] (polynomial degree 11). */
196 -0x8.a41b1e4f36ff8p
-4,
197 0xc.da87d3b69dbe8p
-4,
198 0x1.1474ad5c36709p
+0,
199 0x1.761ecb90c8c5cp
+0,
200 0x1.d279bff588826p
+0,
201 0x2.4e5d003fb36a8p
+0,
202 0x2.d575575566842p
+0,
203 0x3.85152b0d17756p
+0,
205 0x5.55da7dfcf69c4p
+0,
207 0x8.483cc21dd0668p
+0,
208 /* Interval [-3, -2.875] (polynomial degree 11). */
209 -0xa.046d667e468f8p
-4,
215 0x1.04b1365a9adfcp
+0,
216 0x1.22b54ef213798p
+0,
217 0x1.2c52c25206bf5p
+0,
218 0x1.4aa3d798aace4p
+0,
219 0x1.5c3f278b504e3p
+0,
220 0x1.7e08292cc347bp
+0,
223 static const size_t poly_deg
[] =
235 static const size_t poly_end
[] =
247 /* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
253 return __sin (M_PI
* x
);
255 return __cos (M_PI
* (0.5 - x
));
258 /* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
264 return __cos (M_PI
* x
);
266 return __sin (M_PI
* (0.5 - x
));
269 /* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
274 return lg_cospi (x
) / lg_sinpi (x
);
277 /* Compute lgamma of a negative argument -28 < X < -2, setting
278 *SIGNGAMP accordingly. */
281 __lgamma_neg (double x
, int *signgamp
)
283 /* Determine the half-integer region X lies in, handle exact
284 integers and determine the sign of the result. */
285 int i
= floor (-2 * x
);
286 if ((i
& 1) == 0 && i
== -2 * x
)
288 double xn
= ((i
& 1) == 0 ? -i
/ 2 : (-i
- 1) / 2);
290 *signgamp
= ((i
& 2) == 0 ? -1 : 1);
292 SET_RESTORE_ROUND (FE_TONEAREST
);
294 /* Expand around the zero X0 = X0_HI + X0_LO. */
295 double x0_hi
= lgamma_zeros
[i
][0], x0_lo
= lgamma_zeros
[i
][1];
296 double xdiff
= x
- x0_hi
- x0_lo
;
298 /* For arguments in the range -3 to -2, use polynomial
299 approximations to an adjusted version of the gamma function. */
302 int j
= floor (-8 * x
) - 16;
303 double xm
= (-33 - 2 * j
) * 0.0625;
304 double x_adj
= x
- xm
;
305 size_t deg
= poly_deg
[j
];
306 size_t end
= poly_end
[j
];
307 double g
= poly_coeff
[end
];
308 for (size_t j
= 1; j
<= deg
; j
++)
309 g
= g
* x_adj
+ poly_coeff
[end
- j
];
310 return __log1p (g
* xdiff
/ (x
- xn
));
313 /* The result we want is log (sinpi (X0) / sinpi (X))
314 + log (gamma (1 - X0) / gamma (1 - X)). */
315 double x_idiff
= fabs (xn
- x
), x0_idiff
= fabs (xn
- x0_hi
- x0_lo
);
316 double log_sinpi_ratio
;
317 if (x0_idiff
< x_idiff
* 0.5)
318 /* Use log not log1p to avoid inaccuracy from log1p of arguments
320 log_sinpi_ratio
= __ieee754_log (lg_sinpi (x0_idiff
)
321 / lg_sinpi (x_idiff
));
324 /* Use log1p not log to avoid inaccuracy from log of arguments
325 close to 1. X0DIFF2 has positive sign if X0 is further from
326 XN than X is from XN, negative sign otherwise. */
327 double x0diff2
= ((i
& 1) == 0 ? xdiff
: -xdiff
) * 0.5;
328 double sx0d2
= lg_sinpi (x0diff2
);
329 double cx0d2
= lg_cospi (x0diff2
);
330 log_sinpi_ratio
= __log1p (2 * sx0d2
331 * (-sx0d2
+ cx0d2
* lg_cotpi (x_idiff
)));
334 double log_gamma_ratio
;
335 double y0
= math_narrow_eval (1 - x0_hi
);
336 double y0_eps
= -x0_hi
+ (1 - y0
) - x0_lo
;
337 double y
= math_narrow_eval (1 - x
);
338 double y_eps
= -x
+ (1 - y
);
339 /* We now wish to compute LOG_GAMMA_RATIO
340 = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
341 accurately approximates the difference Y0 + Y0_EPS - Y -
342 Y_EPS. Use Stirling's approximation. First, we may need to
343 adjust into the range where Stirling's approximation is
344 sufficiently accurate. */
345 double log_gamma_adj
= 0;
348 int n_up
= (7 - i
) / 2;
349 double ny0
, ny0_eps
, ny
, ny_eps
;
350 ny0
= math_narrow_eval (y0
+ n_up
);
351 ny0_eps
= y0
- (ny0
- n_up
) + y0_eps
;
354 ny
= math_narrow_eval (y
+ n_up
);
355 ny_eps
= y
- (ny
- n_up
) + y_eps
;
358 double prodm1
= __lgamma_product (xdiff
, y
- n_up
, y_eps
, n_up
);
359 log_gamma_adj
= -__log1p (prodm1
);
361 double log_gamma_high
362 = (xdiff
* __log1p ((y0
- e_hi
- e_lo
+ y0_eps
) / e_hi
)
363 + (y
- 0.5 + y_eps
) * __log1p (xdiff
/ y
) + log_gamma_adj
);
364 /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
365 double y0r
= 1 / y0
, yr
= 1 / y
;
366 double y0r2
= y0r
* y0r
, yr2
= yr
* yr
;
367 double rdiff
= -xdiff
/ (y
* y0
);
368 double bterm
[NCOEFF
];
369 double dlast
= rdiff
, elast
= rdiff
* yr
* (yr
+ y0r
);
370 bterm
[0] = dlast
* lgamma_coeff
[0];
371 for (size_t j
= 1; j
< NCOEFF
; j
++)
373 double dnext
= dlast
* y0r2
+ elast
;
374 double enext
= elast
* yr2
;
375 bterm
[j
] = dnext
* lgamma_coeff
[j
];
379 double log_gamma_low
= 0;
380 for (size_t j
= 0; j
< NCOEFF
; j
++)
381 log_gamma_low
+= bterm
[NCOEFF
- 1 - j
];
382 log_gamma_ratio
= log_gamma_high
+ log_gamma_low
;
384 return log_sinpi_ratio
+ log_gamma_ratio
;