1 /* lgammaf 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 float lgamma_zeros
[][2] =
27 { -0x2.74ff94p
+0f
, 0x1.3fe0f2p
-24f
},
28 { -0x2.bf682p
+0f
, -0x1.437b2p
-24f
},
29 { -0x3.24c1b8p
+0f
, 0x6.c34cap
-28f
},
30 { -0x3.f48e2cp
+0f
, 0x1.707a04p
-24f
},
31 { -0x4.0a13ap
+0f
, 0x1.e99aap
-24f
},
32 { -0x4.fdd5ep
+0f
, 0x1.64454p
-24f
},
33 { -0x5.021a98p
+0f
, 0x2.03d248p
-24f
},
34 { -0x5.ffa4cp
+0f
, 0x2.9b82fcp
-24f
},
35 { -0x6.005ac8p
+0f
, -0x1.625f24p
-24f
},
36 { -0x6.fff3p
+0f
, 0x2.251e44p
-24f
},
37 { -0x7.000dp
+0f
, 0x8.48078p
-28f
},
38 { -0x7.fffe6p
+0f
, 0x1.fa98c4p
-28f
},
39 { -0x8.0001ap
+0f
, -0x1.459fcap
-28f
},
40 { -0x8.ffffdp
+0f
, -0x1.c425e8p
-24f
},
41 { -0x9.00003p
+0f
, 0x1.c44b82p
-24f
},
42 { -0xap
+0f
, 0x4.9f942p
-24f
},
43 { -0xap
+0f
, -0x4.9f93b8p
-24f
},
44 { -0xbp
+0f
, 0x6.b9916p
-28f
},
45 { -0xbp
+0f
, -0x6.b9915p
-28f
},
46 { -0xcp
+0f
, 0x8.f76c8p
-32f
},
47 { -0xcp
+0f
, -0x8.f76c7p
-32f
},
48 { -0xdp
+0f
, 0xb.09231p
-36f
},
49 { -0xdp
+0f
, -0xb.09231p
-36f
},
50 { -0xep
+0f
, 0xc.9cba5p
-40f
},
51 { -0xep
+0f
, -0xc.9cba5p
-40f
},
52 { -0xfp
+0f
, 0xd.73f9fp
-44f
},
55 static const float e_hi
= 0x2.b7e15p
+0f
, e_lo
= 0x1.628aeep
-24f
;
57 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
58 approximation to lgamma function. */
60 static const float lgamma_coeff
[] =
67 #define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
69 /* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
70 the integer end-point of the half-integer interval containing x and
71 x0 is the zero of lgamma in that half-integer interval. Each
72 polynomial is expressed in terms of x-xm, where xm is the midpoint
73 of the interval for which the polynomial applies. */
75 static const float poly_coeff
[] =
77 /* Interval [-2.125, -2] (polynomial degree 5). */
84 /* Interval [-2.25, -2.125] (polynomial degree 5). */
91 /* Interval [-2.375, -2.25] (polynomial degree 5). */
98 /* Interval [-2.5, -2.375] (polynomial degree 6). */
106 /* Interval [-2.625, -2.5] (polynomial degree 6). */
114 /* Interval [-2.75, -2.625] (polynomial degree 6). */
122 /* Interval [-2.875, -2.75] (polynomial degree 5). */
129 /* Interval [-3, -2.875] (polynomial degree 5). */
138 static const size_t poly_deg
[] =
150 static const size_t poly_end
[] =
162 /* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
168 return __sinf (M_PIf
* x
);
170 return __cosf (M_PIf
* (0.5f
- x
));
173 /* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
179 return __cosf (M_PIf
* x
);
181 return __sinf (M_PIf
* (0.5f
- x
));
184 /* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
189 return lg_cospi (x
) / lg_sinpi (x
);
192 /* Compute lgamma of a negative argument -15 < X < -2, setting
193 *SIGNGAMP accordingly. */
196 __lgamma_negf (float x
, int *signgamp
)
198 /* Determine the half-integer region X lies in, handle exact
199 integers and determine the sign of the result. */
200 int i
= floorf (-2 * x
);
201 if ((i
& 1) == 0 && i
== -2 * x
)
203 float xn
= ((i
& 1) == 0 ? -i
/ 2 : (-i
- 1) / 2);
205 *signgamp
= ((i
& 2) == 0 ? -1 : 1);
207 SET_RESTORE_ROUNDF (FE_TONEAREST
);
209 /* Expand around the zero X0 = X0_HI + X0_LO. */
210 float x0_hi
= lgamma_zeros
[i
][0], x0_lo
= lgamma_zeros
[i
][1];
211 float xdiff
= x
- x0_hi
- x0_lo
;
213 /* For arguments in the range -3 to -2, use polynomial
214 approximations to an adjusted version of the gamma function. */
217 int j
= floorf (-8 * x
) - 16;
218 float xm
= (-33 - 2 * j
) * 0.0625f
;
219 float x_adj
= x
- xm
;
220 size_t deg
= poly_deg
[j
];
221 size_t end
= poly_end
[j
];
222 float g
= poly_coeff
[end
];
223 for (size_t j
= 1; j
<= deg
; j
++)
224 g
= g
* x_adj
+ poly_coeff
[end
- j
];
225 return __log1pf (g
* xdiff
/ (x
- xn
));
228 /* The result we want is log (sinpi (X0) / sinpi (X))
229 + log (gamma (1 - X0) / gamma (1 - X)). */
230 float x_idiff
= fabsf (xn
- x
), x0_idiff
= fabsf (xn
- x0_hi
- x0_lo
);
231 float log_sinpi_ratio
;
232 if (x0_idiff
< x_idiff
* 0.5f
)
233 /* Use log not log1p to avoid inaccuracy from log1p of arguments
235 log_sinpi_ratio
= __ieee754_logf (lg_sinpi (x0_idiff
)
236 / lg_sinpi (x_idiff
));
239 /* Use log1p not log to avoid inaccuracy from log of arguments
240 close to 1. X0DIFF2 has positive sign if X0 is further from
241 XN than X is from XN, negative sign otherwise. */
242 float x0diff2
= ((i
& 1) == 0 ? xdiff
: -xdiff
) * 0.5f
;
243 float sx0d2
= lg_sinpi (x0diff2
);
244 float cx0d2
= lg_cospi (x0diff2
);
245 log_sinpi_ratio
= __log1pf (2 * sx0d2
246 * (-sx0d2
+ cx0d2
* lg_cotpi (x_idiff
)));
249 float log_gamma_ratio
;
250 float y0
= math_narrow_eval (1 - x0_hi
);
251 float y0_eps
= -x0_hi
+ (1 - y0
) - x0_lo
;
252 float y
= math_narrow_eval (1 - x
);
253 float y_eps
= -x
+ (1 - y
);
254 /* We now wish to compute LOG_GAMMA_RATIO
255 = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
256 accurately approximates the difference Y0 + Y0_EPS - Y -
257 Y_EPS. Use Stirling's approximation. */
259 = (xdiff
* __log1pf ((y0
- e_hi
- e_lo
+ y0_eps
) / e_hi
)
260 + (y
- 0.5f
+ y_eps
) * __log1pf (xdiff
/ y
));
261 /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
262 float y0r
= 1 / y0
, yr
= 1 / y
;
263 float y0r2
= y0r
* y0r
, yr2
= yr
* yr
;
264 float rdiff
= -xdiff
/ (y
* y0
);
266 float dlast
= rdiff
, elast
= rdiff
* yr
* (yr
+ y0r
);
267 bterm
[0] = dlast
* lgamma_coeff
[0];
268 for (size_t j
= 1; j
< NCOEFF
; j
++)
270 float dnext
= dlast
* y0r2
+ elast
;
271 float enext
= elast
* yr2
;
272 bterm
[j
] = dnext
* lgamma_coeff
[j
];
276 float log_gamma_low
= 0;
277 for (size_t j
= 0; j
< NCOEFF
; j
++)
278 log_gamma_low
+= bterm
[NCOEFF
- 1 - j
];
279 log_gamma_ratio
= log_gamma_high
+ log_gamma_low
;
281 return log_sinpi_ratio
+ log_gamma_ratio
;