1 /* lgammaf expanding around zeros.
2 Copyright (C) 2015-2016 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 <http://www.gnu.org/licenses/>. */
21 #include <math_private.h>
23 static const float lgamma_zeros
[][2] =
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
},
53 static const float e_hi
= 0x2.b7e15p
+0f
, e_lo
= 0x1.628aeep
-24f
;
55 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
56 approximation to lgamma function. */
58 static const float lgamma_coeff
[] =
65 #define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
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. */
73 static const float poly_coeff
[] =
75 /* Interval [-2.125, -2] (polynomial degree 5). */
82 /* Interval [-2.25, -2.125] (polynomial degree 5). */
89 /* Interval [-2.375, -2.25] (polynomial degree 5). */
96 /* Interval [-2.5, -2.375] (polynomial degree 6). */
104 /* Interval [-2.625, -2.5] (polynomial degree 6). */
112 /* Interval [-2.75, -2.625] (polynomial degree 6). */
120 /* Interval [-2.875, -2.75] (polynomial degree 5). */
127 /* Interval [-3, -2.875] (polynomial degree 5). */
136 static const size_t poly_deg
[] =
148 static const size_t poly_end
[] =
160 /* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
166 return __sinf ((float) M_PI
* x
);
168 return __cosf ((float) M_PI
* (0.5f
- x
));
171 /* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
177 return __cosf ((float) M_PI
* x
);
179 return __sinf ((float) M_PI
* (0.5f
- x
));
182 /* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
187 return lg_cospi (x
) / lg_sinpi (x
);
190 /* Compute lgamma of a negative argument -15 < X < -2, setting
191 *SIGNGAMP accordingly. */
194 __lgamma_negf (float x
, int *signgamp
)
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
)
201 float xn
= ((i
& 1) == 0 ? -i
/ 2 : (-i
- 1) / 2);
203 *signgamp
= ((i
& 2) == 0 ? -1 : 1);
205 SET_RESTORE_ROUNDF (FE_TONEAREST
);
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
;
211 /* For arguments in the range -3 to -2, use polynomial
212 approximations to an adjusted version of the gamma function. */
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
));
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
233 log_sinpi_ratio
= __ieee754_logf (lg_sinpi (x0_idiff
)
234 / lg_sinpi (x_idiff
));
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
)));
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. */
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
);
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
++)
268 float dnext
= dlast
* y0r2
+ elast
;
269 float enext
= elast
* yr2
;
270 bterm
[j
] = dnext
* lgamma_coeff
[j
];
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
;
279 return log_sinpi_ratio
+ log_gamma_ratio
;