1 /* lgammaf expanding around zeros.
2 Copyright (C) 2015-2023 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4
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.
9
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.
14
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/>. */
18
19 #include <float.h>
20 #include <math.h>
21 #include <math-narrow-eval.h>
22 #include <math_private.h>
23 #include <fenv_private.h>
24
25 static const float lgamma_zeros[][2] =
26 {
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 },
53 };
54
55 static const float e_hi = 0x2.b7e15p+0f, e_lo = 0x1.628aeep-24f;
56
57 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
58 approximation to lgamma function. */
59
60 static const float lgamma_coeff[] =
61 {
62 0x1.555556p-4f,
63 -0xb.60b61p-12f,
64 0x3.403404p-12f,
65 };
66
67 #define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
68
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. */
74
75 static const float poly_coeff[] =
76 {
77 /* Interval [-2.125, -2] (polynomial degree 5). */
78 -0x1.0b71c6p+0f,
79 -0xc.73a1ep-4f,
80 -0x1.ec8462p-4f,
81 -0xe.37b93p-4f,
82 -0x1.02ed36p-4f,
83 -0xe.cbe26p-4f,
84 /* Interval [-2.25, -2.125] (polynomial degree 5). */
85 -0xf.29309p-4f,
86 -0xc.a5cfep-4f,
87 0x3.9c93fcp-4f,
88 -0x1.02a2fp+0f,
89 0x9.896bep-4f,
90 -0x1.519704p+0f,
91 /* Interval [-2.375, -2.25] (polynomial degree 5). */
92 -0xd.7d28dp-4f,
93 -0xe.6964cp-4f,
94 0xb.0d4f1p-4f,
95 -0x1.9240aep+0f,
96 0x1.dadabap+0f,
97 -0x3.1778c4p+0f,
98 /* Interval [-2.5, -2.375] (polynomial degree 6). */
99 -0xb.74ea2p-4f,
100 -0x1.2a82cp+0f,
101 0x1.880234p+0f,
102 -0x3.320c4p+0f,
103 0x5.572a38p+0f,
104 -0x9.f92bap+0f,
105 0x1.1c347ep+4f,
106 /* Interval [-2.625, -2.5] (polynomial degree 6). */
107 -0x3.d10108p-4f,
108 0x1.cd5584p+0f,
109 0x3.819c24p+0f,
110 0x6.84cbb8p+0f,
111 0xb.bf269p+0f,
112 0x1.57fb12p+4f,
113 0x2.7b9854p+4f,
114 /* Interval [-2.75, -2.625] (polynomial degree 6). */
115 -0x6.b5d25p-4f,
116 0x1.28d604p+0f,
117 0x1.db6526p+0f,
118 0x2.e20b38p+0f,
119 0x4.44c378p+0f,
120 0x6.62a08p+0f,
121 0x9.6db3ap+0f,
122 /* Interval [-2.875, -2.75] (polynomial degree 5). */
123 -0x8.a41b2p-4f,
124 0xc.da87fp-4f,
125 0x1.147312p+0f,
126 0x1.7617dap+0f,
127 0x1.d6c13p+0f,
128 0x2.57a358p+0f,
129 /* Interval [-3, -2.875] (polynomial degree 5). */
130 -0xa.046d6p-4f,
131 0x9.70b89p-4f,
132 0xa.a89a6p-4f,
133 0xd.2f2d8p-4f,
134 0xd.e32b4p-4f,
135 0xf.fb741p-4f,
136 };
137
138 static const size_t poly_deg[] =
139 {
140 5,
141 5,
142 5,
143 6,
144 6,
145 6,
146 5,
147 5,
148 };
149
150 static const size_t poly_end[] =
151 {
152 5,
153 11,
154 17,
155 24,
156 31,
157 38,
158 44,
159 50,
160 };
161
162 /* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
163
164 static float
165 lg_sinpi (float x)
166 {
167 if (x <= 0.25f)
168 return __sinf (M_PIf * x);
169 else
170 return __cosf (M_PIf * (0.5f - x));
171 }
172
173 /* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
174
175 static float
176 lg_cospi (float x)
177 {
178 if (x <= 0.25f)
179 return __cosf (M_PIf * x);
180 else
181 return __sinf (M_PIf * (0.5f - x));
182 }
183
184 /* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
185
186 static float
187 lg_cotpi (float x)
188 {
189 return lg_cospi (x) / lg_sinpi (x);
190 }
191
192 /* Compute lgamma of a negative argument -15 < X < -2, setting
193 *SIGNGAMP accordingly. */
194
195 float
196 __lgamma_negf (float x, int *signgamp)
197 {
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)
202 return 1.0f / 0.0f;
203 float xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
204 i -= 4;
205 *signgamp = ((i & 2) == 0 ? -1 : 1);
206
207 SET_RESTORE_ROUNDF (FE_TONEAREST);
208
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;
212
213 /* For arguments in the range -3 to -2, use polynomial
214 approximations to an adjusted version of the gamma function. */
215 if (i < 2)
216 {
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));
226 }
227
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
234 close to -1. */
235 log_sinpi_ratio = __ieee754_logf (lg_sinpi (x0_idiff)
236 / lg_sinpi (x_idiff));
237 else
238 {
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)));
247 }
248
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. */
258 float log_gamma_high
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);
265 float bterm[NCOEFF];
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++)
269 {
270 float dnext = dlast * y0r2 + elast;
271 float enext = elast * yr2;
272 bterm[j] = dnext * lgamma_coeff[j];
273 dlast = dnext;
274 elast = enext;
275 }
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;
280
281 return log_sinpi_ratio + log_gamma_ratio;
282 }