1 | /* lgammaf expanding around zeros. |
2 | Copyright (C) 2015-2018 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 | <http://www.gnu.org/licenses/>. */ |
18 | |
19 | #include <float.h> |
20 | #include <math.h> |
21 | #include <math_private.h> |
22 | |
23 | static const float lgamma_zeros[][2] = |
24 | { |
25 | { -0x2.74ff94p+0f, 0x1.3fe0f2p-24f }, |
26 | { -0x2.bf682p+0f, -0x1.437b2p-24f }, |
27 | { -0x3.24c1b8p+0f, 0x6.c34cap-28f }, |
28 | { -0x3.f48e2cp+0f, 0x1.707a04p-24f }, |
29 | { -0x4.0a13ap+0f, 0x1.e99aap-24f }, |
30 | { -0x4.fdd5ep+0f, 0x1.64454p-24f }, |
31 | { -0x5.021a98p+0f, 0x2.03d248p-24f }, |
32 | { -0x5.ffa4cp+0f, 0x2.9b82fcp-24f }, |
33 | { -0x6.005ac8p+0f, -0x1.625f24p-24f }, |
34 | { -0x6.fff3p+0f, 0x2.251e44p-24f }, |
35 | { -0x7.000dp+0f, 0x8.48078p-28f }, |
36 | { -0x7.fffe6p+0f, 0x1.fa98c4p-28f }, |
37 | { -0x8.0001ap+0f, -0x1.459fcap-28f }, |
38 | { -0x8.ffffdp+0f, -0x1.c425e8p-24f }, |
39 | { -0x9.00003p+0f, 0x1.c44b82p-24f }, |
40 | { -0xap+0f, 0x4.9f942p-24f }, |
41 | { -0xap+0f, -0x4.9f93b8p-24f }, |
42 | { -0xbp+0f, 0x6.b9916p-28f }, |
43 | { -0xbp+0f, -0x6.b9915p-28f }, |
44 | { -0xcp+0f, 0x8.f76c8p-32f }, |
45 | { -0xcp+0f, -0x8.f76c7p-32f }, |
46 | { -0xdp+0f, 0xb.09231p-36f }, |
47 | { -0xdp+0f, -0xb.09231p-36f }, |
48 | { -0xep+0f, 0xc.9cba5p-40f }, |
49 | { -0xep+0f, -0xc.9cba5p-40f }, |
50 | { -0xfp+0f, 0xd.73f9fp-44f }, |
51 | }; |
52 | |
53 | static const float e_hi = 0x2.b7e15p+0f, e_lo = 0x1.628aeep-24f; |
54 | |
55 | /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's |
56 | approximation to lgamma function. */ |
57 | |
58 | static const float lgamma_coeff[] = |
59 | { |
60 | 0x1.555556p-4f, |
61 | -0xb.60b61p-12f, |
62 | 0x3.403404p-12f, |
63 | }; |
64 | |
65 | #define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0])) |
66 | |
67 | /* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is |
68 | the integer end-point of the half-integer interval containing x and |
69 | x0 is the zero of lgamma in that half-integer interval. Each |
70 | polynomial is expressed in terms of x-xm, where xm is the midpoint |
71 | of the interval for which the polynomial applies. */ |
72 | |
73 | static const float poly_coeff[] = |
74 | { |
75 | /* Interval [-2.125, -2] (polynomial degree 5). */ |
76 | -0x1.0b71c6p+0f, |
77 | -0xc.73a1ep-4f, |
78 | -0x1.ec8462p-4f, |
79 | -0xe.37b93p-4f, |
80 | -0x1.02ed36p-4f, |
81 | -0xe.cbe26p-4f, |
82 | /* Interval [-2.25, -2.125] (polynomial degree 5). */ |
83 | -0xf.29309p-4f, |
84 | -0xc.a5cfep-4f, |
85 | 0x3.9c93fcp-4f, |
86 | -0x1.02a2fp+0f, |
87 | 0x9.896bep-4f, |
88 | -0x1.519704p+0f, |
89 | /* Interval [-2.375, -2.25] (polynomial degree 5). */ |
90 | -0xd.7d28dp-4f, |
91 | -0xe.6964cp-4f, |
92 | 0xb.0d4f1p-4f, |
93 | -0x1.9240aep+0f, |
94 | 0x1.dadabap+0f, |
95 | -0x3.1778c4p+0f, |
96 | /* Interval [-2.5, -2.375] (polynomial degree 6). */ |
97 | -0xb.74ea2p-4f, |
98 | -0x1.2a82cp+0f, |
99 | 0x1.880234p+0f, |
100 | -0x3.320c4p+0f, |
101 | 0x5.572a38p+0f, |
102 | -0x9.f92bap+0f, |
103 | 0x1.1c347ep+4f, |
104 | /* Interval [-2.625, -2.5] (polynomial degree 6). */ |
105 | -0x3.d10108p-4f, |
106 | 0x1.cd5584p+0f, |
107 | 0x3.819c24p+0f, |
108 | 0x6.84cbb8p+0f, |
109 | 0xb.bf269p+0f, |
110 | 0x1.57fb12p+4f, |
111 | 0x2.7b9854p+4f, |
112 | /* Interval [-2.75, -2.625] (polynomial degree 6). */ |
113 | -0x6.b5d25p-4f, |
114 | 0x1.28d604p+0f, |
115 | 0x1.db6526p+0f, |
116 | 0x2.e20b38p+0f, |
117 | 0x4.44c378p+0f, |
118 | 0x6.62a08p+0f, |
119 | 0x9.6db3ap+0f, |
120 | /* Interval [-2.875, -2.75] (polynomial degree 5). */ |
121 | -0x8.a41b2p-4f, |
122 | 0xc.da87fp-4f, |
123 | 0x1.147312p+0f, |
124 | 0x1.7617dap+0f, |
125 | 0x1.d6c13p+0f, |
126 | 0x2.57a358p+0f, |
127 | /* Interval [-3, -2.875] (polynomial degree 5). */ |
128 | -0xa.046d6p-4f, |
129 | 0x9.70b89p-4f, |
130 | 0xa.a89a6p-4f, |
131 | 0xd.2f2d8p-4f, |
132 | 0xd.e32b4p-4f, |
133 | 0xf.fb741p-4f, |
134 | }; |
135 | |
136 | static const size_t poly_deg[] = |
137 | { |
138 | 5, |
139 | 5, |
140 | 5, |
141 | 6, |
142 | 6, |
143 | 6, |
144 | 5, |
145 | 5, |
146 | }; |
147 | |
148 | static const size_t poly_end[] = |
149 | { |
150 | 5, |
151 | 11, |
152 | 17, |
153 | 24, |
154 | 31, |
155 | 38, |
156 | 44, |
157 | 50, |
158 | }; |
159 | |
160 | /* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ |
161 | |
162 | static float |
163 | lg_sinpi (float x) |
164 | { |
165 | if (x <= 0.25f) |
166 | return __sinf ((float) M_PI * x); |
167 | else |
168 | return __cosf ((float) M_PI * (0.5f - x)); |
169 | } |
170 | |
171 | /* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ |
172 | |
173 | static float |
174 | lg_cospi (float x) |
175 | { |
176 | if (x <= 0.25f) |
177 | return __cosf ((float) M_PI * x); |
178 | else |
179 | return __sinf ((float) M_PI * (0.5f - x)); |
180 | } |
181 | |
182 | /* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ |
183 | |
184 | static float |
185 | lg_cotpi (float x) |
186 | { |
187 | return lg_cospi (x) / lg_sinpi (x); |
188 | } |
189 | |
190 | /* Compute lgamma of a negative argument -15 < X < -2, setting |
191 | *SIGNGAMP accordingly. */ |
192 | |
193 | float |
194 | __lgamma_negf (float x, int *signgamp) |
195 | { |
196 | /* Determine the half-integer region X lies in, handle exact |
197 | integers and determine the sign of the result. */ |
198 | int i = __floorf (-2 * x); |
199 | if ((i & 1) == 0 && i == -2 * x) |
200 | return 1.0f / 0.0f; |
201 | float xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); |
202 | i -= 4; |
203 | *signgamp = ((i & 2) == 0 ? -1 : 1); |
204 | |
205 | SET_RESTORE_ROUNDF (FE_TONEAREST); |
206 | |
207 | /* Expand around the zero X0 = X0_HI + X0_LO. */ |
208 | float x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; |
209 | float xdiff = x - x0_hi - x0_lo; |
210 | |
211 | /* For arguments in the range -3 to -2, use polynomial |
212 | approximations to an adjusted version of the gamma function. */ |
213 | if (i < 2) |
214 | { |
215 | int j = __floorf (-8 * x) - 16; |
216 | float xm = (-33 - 2 * j) * 0.0625f; |
217 | float x_adj = x - xm; |
218 | size_t deg = poly_deg[j]; |
219 | size_t end = poly_end[j]; |
220 | float g = poly_coeff[end]; |
221 | for (size_t j = 1; j <= deg; j++) |
222 | g = g * x_adj + poly_coeff[end - j]; |
223 | return __log1pf (g * xdiff / (x - xn)); |
224 | } |
225 | |
226 | /* The result we want is log (sinpi (X0) / sinpi (X)) |
227 | + log (gamma (1 - X0) / gamma (1 - X)). */ |
228 | float x_idiff = fabsf (xn - x), x0_idiff = fabsf (xn - x0_hi - x0_lo); |
229 | float log_sinpi_ratio; |
230 | if (x0_idiff < x_idiff * 0.5f) |
231 | /* Use log not log1p to avoid inaccuracy from log1p of arguments |
232 | close to -1. */ |
233 | log_sinpi_ratio = __ieee754_logf (lg_sinpi (x0_idiff) |
234 | / lg_sinpi (x_idiff)); |
235 | else |
236 | { |
237 | /* Use log1p not log to avoid inaccuracy from log of arguments |
238 | close to 1. X0DIFF2 has positive sign if X0 is further from |
239 | XN than X is from XN, negative sign otherwise. */ |
240 | float x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5f; |
241 | float sx0d2 = lg_sinpi (x0diff2); |
242 | float cx0d2 = lg_cospi (x0diff2); |
243 | log_sinpi_ratio = __log1pf (2 * sx0d2 |
244 | * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); |
245 | } |
246 | |
247 | float log_gamma_ratio; |
248 | float y0 = math_narrow_eval (1 - x0_hi); |
249 | float y0_eps = -x0_hi + (1 - y0) - x0_lo; |
250 | float y = math_narrow_eval (1 - x); |
251 | float y_eps = -x + (1 - y); |
252 | /* We now wish to compute LOG_GAMMA_RATIO |
253 | = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF |
254 | accurately approximates the difference Y0 + Y0_EPS - Y - |
255 | Y_EPS. Use Stirling's approximation. */ |
256 | float log_gamma_high |
257 | = (xdiff * __log1pf ((y0 - e_hi - e_lo + y0_eps) / e_hi) |
258 | + (y - 0.5f + y_eps) * __log1pf (xdiff / y)); |
259 | /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ |
260 | float y0r = 1 / y0, yr = 1 / y; |
261 | float y0r2 = y0r * y0r, yr2 = yr * yr; |
262 | float rdiff = -xdiff / (y * y0); |
263 | float bterm[NCOEFF]; |
264 | float dlast = rdiff, elast = rdiff * yr * (yr + y0r); |
265 | bterm[0] = dlast * lgamma_coeff[0]; |
266 | for (size_t j = 1; j < NCOEFF; j++) |
267 | { |
268 | float dnext = dlast * y0r2 + elast; |
269 | float enext = elast * yr2; |
270 | bterm[j] = dnext * lgamma_coeff[j]; |
271 | dlast = dnext; |
272 | elast = enext; |
273 | } |
274 | float log_gamma_low = 0; |
275 | for (size_t j = 0; j < NCOEFF; j++) |
276 | log_gamma_low += bterm[NCOEFF - 1 - j]; |
277 | log_gamma_ratio = log_gamma_high + log_gamma_low; |
278 | |
279 | return log_sinpi_ratio + log_gamma_ratio; |
280 | } |
281 | |