1 /* mpfr_lngamma -- lngamma function
3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* assuming b[0]...b[2(n-1)] are computed, computes and stores B[2n]*(2n+1)!
28 t/(exp(t)-1) = sum(B[j]*t^j/j!, j=0..infinity)
29 thus t = (exp(t)-1) * sum(B[j]*t^j/j!, n=0..infinity).
30 Taking the coefficient of degree n+1 > 1, we get:
31 0 = sum(1/(n+1-k)!*B[k]/k!, k=0..n)
33 B[n] = -sum(binomial(n+1,k)*B[k], k=0..n-1)/(n+1).
35 Let C[n] = B[n]*(n+1)!.
36 Then C[n] = -sum(binomial(n+1,k)*C[k]*n!/(k+1)!, k=0..n-1),
37 which proves that the C[n] are integers.
40 bernoulli (mpz_t *b, unsigned long n)
44 b = (mpz_t *) (*__gmp_allocate_func) (sizeof (mpz_t));
45 mpz_init_set_ui (b[0], 1);
52 b = (mpz_t *) (*__gmp_reallocate_func)
53 (b, n * sizeof (mpz_t), (n + 1) * sizeof (mpz_t));
55 /* b[n] = -sum(binomial(2n+1,2k)*C[k]*(2n)!/(2k+1)!, k=0..n-1) */
56 mpz_init_set_ui (t, 2 * n + 1);
57 mpz_mul_ui (t, t, 2 * n - 1);
58 mpz_mul_ui (t, t, 2 * n);
60 mpz_div_ui (t, t, 3); /* exact: t=binomial(2*n+1,2*k)*(2*n)!/(2*k+1)!
62 mpz_mul (b[n], t, b[n-1]);
63 for (k = n - 1; k-- > 0;)
65 mpz_mul_ui (t, t, 2 * k + 1);
66 mpz_mul_ui (t, t, 2 * k + 2);
67 mpz_mul_ui (t, t, 2 * k + 2);
68 mpz_mul_ui (t, t, 2 * k + 3);
69 mpz_div_ui (t, t, 2 * (n - k) + 1);
70 mpz_div_ui (t, t, 2 * (n - k));
71 mpz_addmul (b[n], t, b[k]);
73 /* take into account C[1] */
74 mpz_mul_ui (t, t, 2 * n + 1);
75 mpz_div_2exp (t, t, 1);
76 mpz_sub (b[n], b[n], t);
83 /* given a precision p, return alpha, such that the argument reduction
84 will use k = alpha*p*log(2).
86 Warning: we should always have alpha >= log(2)/(2Pi) ~ 0.11,
87 and the smallest value of alpha multiplied by the smallest working
88 precision should be >= 4.
91 mpfr_gamma_alpha (mpfr_t s, mp_prec_t p)
94 mpfr_set_ui_2exp (s, 614, -10, GMP_RNDN); /* about 0.6 */
96 mpfr_set_ui_2exp (s, 819, -10, GMP_RNDN); /* about 0.8 */
98 mpfr_set_ui_2exp (s, 1331, -10, GMP_RNDN); /* about 1.3 */
100 mpfr_set_ui_2exp (s, 1741, -10, GMP_RNDN); /* about 1.7 */
102 mpfr_set_ui_2exp (s, 2253, -10, GMP_RNDN); /* about 2.2 */
104 mpfr_set_ui_2exp (s, 3482, -10, GMP_RNDN); /* about 3.4 */
106 mpfr_set_ui_2exp (s, 9, -1, GMP_RNDN); /* 4.5 */
111 unit_bit (mpfr_srcptr (x))
117 expo = MPFR_GET_EXP (x);
119 return 0; /* |x| < 1 */
121 prec = MPFR_PREC (x);
123 return 0; /* y is a multiple of 2^(expo-prec), thus an even integer */
125 /* Now, the unit bit is represented. */
127 prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
128 /* number of represented fractional bits (including the trailing 0's) */
130 x0 = *(MPFR_MANT (x) + prec / BITS_PER_MP_LIMB);
131 /* limb containing the unit bit */
133 return (x0 >> (prec % BITS_PER_MP_LIMB)) & 1;
137 /* lngamma(x) = log(gamma(x)).
138 We use formula [6.1.40] from Abramowitz&Stegun:
139 lngamma(z) = (z-1/2)*log(z) - z + 1/2*log(2*Pi)
140 + sum (Bernoulli[2m]/(2m)/(2m-1)/z^(2m-1),m=1..infinity)
141 According to [6.1.42], if the sum is truncated after m=n, the error
142 R_n(z) is bounded by |B[2n+2]|*K(z)/(2n+1)/(2n+2)/|z|^(2n+1)
143 where K(z) = max (z^2/(u^2+z^2)) for u >= 0.
144 For z real, |K(z)| <= 1 thus R_n(z) is bounded by the first neglected term.
147 #define GAMMA_FUNC mpfr_gamma_aux
149 #define GAMMA_FUNC mpfr_lngamma_aux
153 GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd)
155 mp_prec_t precy, w; /* working precision */
156 mpfr_t s, t, u, v, z;
157 unsigned long m, k, maxm;
158 mpz_t *INITIALIZED(B); /* variable B declared as initialized */
159 int inexact, compared;
160 mp_exp_t err_s, err_t;
161 unsigned long Bm = 0; /* number of allocated B[] */
164 MPFR_SAVE_EXPO_DECL (expo);
166 compared = mpfr_cmp_ui (z0, 1);
168 MPFR_SAVE_EXPO_MARK (expo);
170 #ifndef IS_GAMMA /* lngamma or lgamma */
171 if (compared == 0 || (compared > 0 && mpfr_cmp_ui (z0, 2) == 0))
173 MPFR_SAVE_EXPO_FREE (expo);
174 return mpfr_set_ui (y, 0, GMP_RNDN); /* lngamma(1 or 2) = +0 */
177 /* Deal here with tiny inputs. We have for -0.3 <= x <= 0.3:
178 - log|x| - gamma*x <= log|gamma(x)| <= - log|x| - gamma*x + x^2 */
179 if (MPFR_EXP(z0) <= - (mp_exp_t) MPFR_PREC(y))
183 mp_prec_t prec = MPFR_PREC(y) + 14;
184 MPFR_ZIV_DECL (loop);
186 MPFR_ZIV_INIT (loop, prec);
189 mpfr_init2 (l, prec);
192 mpfr_log (l, z0, GMP_RNDU); /* upper bound for log(z0) */
193 mpfr_init2 (h, MPFR_PREC(l));
197 mpfr_init2 (h, MPFR_PREC(z0));
198 mpfr_neg (h, z0, GMP_RNDN); /* exact */
199 mpfr_log (l, h, GMP_RNDU); /* upper bound for log(-z0) */
200 mpfr_set_prec (h, MPFR_PREC(l));
202 mpfr_neg (l, l, GMP_RNDD); /* lower bound for -log(|z0|) */
203 mpfr_set (h, l, GMP_RNDD); /* exact */
204 mpfr_nextabove (h); /* upper bound for -log(|z0|), avoids two calls
206 mpfr_init2 (g, MPFR_PREC(l));
207 /* if z0>0, we need an upper approximation of Euler's constant
208 for the left bound */
209 mpfr_const_euler (g, MPFR_IS_POS(z0) ? GMP_RNDU : GMP_RNDD);
210 mpfr_mul (g, g, z0, GMP_RNDD);
211 mpfr_sub (l, l, g, GMP_RNDD);
212 mpfr_const_euler (g, MPFR_IS_POS(z0) ? GMP_RNDD : GMP_RNDU); /* cached */
213 mpfr_mul (g, g, z0, GMP_RNDU);
214 mpfr_sub (h, h, g, GMP_RNDD);
215 mpfr_mul (g, z0, z0, GMP_RNDU);
216 mpfr_add (h, h, g, GMP_RNDU);
217 inexact = mpfr_prec_round (l, MPFR_PREC(y), rnd);
218 inex2 = mpfr_prec_round (h, MPFR_PREC(y), rnd);
219 /* Caution: we not only need l = h, but both inexact flags should
220 agree. Indeed, one of the inexact flags might be zero. In that
221 case if we assume lngamma(z0) cannot be exact, the other flag
222 should be correct. We are conservative here and request that both
223 inexact flags agree. */
224 ok = SAME_SIGN (inexact, inex2) && mpfr_cmp (l, h) == 0;
226 mpfr_set (y, h, rnd); /* exact */
232 MPFR_SAVE_EXPO_FREE (expo);
233 return mpfr_check_range (y, inexact, rnd);
235 /* since we have log|gamma(x)| = - log|x| - gamma*x + O(x^2),
236 if x ~ 2^(-n), then we have a n-bit approximation, thus
237 we can try again with a working precision of n bits,
238 especially when n >> PREC(y).
239 Otherwise we would use the reflection formula evaluating x-1,
240 which would need precision n. */
241 MPFR_ZIV_NEXT (loop, prec);
243 while (prec <= -MPFR_EXP(z0));
244 MPFR_ZIV_FREE (loop);
248 precy = MPFR_PREC(y);
250 mpfr_init2 (s, MPFR_PREC_MIN);
251 mpfr_init2 (t, MPFR_PREC_MIN);
252 mpfr_init2 (u, MPFR_PREC_MIN);
253 mpfr_init2 (v, MPFR_PREC_MIN);
254 mpfr_init2 (z, MPFR_PREC_MIN);
260 /* use reflection formula:
261 gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x)
262 thus lngamma(x) = log(Pi*(x-1)/sin(Pi*(2-x))) - lngamma(2-x) */
264 w = precy + MPFR_INT_CEIL_LOG2 (precy);
267 w += MPFR_INT_CEIL_LOG2 (w) + 14;
268 MPFR_ASSERTD(w >= 3);
269 mpfr_set_prec (s, w);
270 mpfr_set_prec (t, w);
271 mpfr_set_prec (u, w);
272 mpfr_set_prec (v, w);
273 /* In the following, we write r for a real of absolute value
274 at most 2^(-w). Different instances of r may represent different
276 mpfr_ui_sub (s, 2, z0, GMP_RNDD); /* s = (2-z0) * (1+2r) >= 1 */
277 mpfr_const_pi (t, GMP_RNDN); /* t = Pi * (1+r) */
278 mpfr_lngamma (u, s, GMP_RNDN); /* lngamma(2-x) */
279 /* Let s = (2-z0) + h. By construction, -(2-z0)*2^(1-w) <= h <= 0.
280 We have lngamma(s) = lngamma(2-z0) + h*Psi(z), z in [2-z0+h,2-z0].
281 Since 2-z0+h = s >= 1 and |Psi(x)| <= max(1,log(x)) for x >= 1,
282 the error on u is bounded by
283 ulp(u)/2 + (2-z0)*max(1,log(2-z0))*2^(1-w)
284 = (1/2 + (2-z0)*max(1,log(2-z0))*2^(1-E(u))) ulp(u) */
285 d = (double) MPFR_GET_EXP(s) * 0.694; /* upper bound for log(2-z0) */
286 err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u);
287 err_u = (err_u >= 0) ? err_u + 1 : 0;
288 /* now the error on u is bounded by 2^err_u ulps */
290 mpfr_mul (s, s, t, GMP_RNDN); /* Pi*(2-x) * (1+r)^4 */
291 err_s = MPFR_GET_EXP(s); /* 2-x <= 2^err_s */
292 mpfr_sin (s, s, GMP_RNDN); /* sin(Pi*(2-x)) */
293 /* the error on s is bounded by 1/2*ulp(s) + [(1+2^(-w))^4-1]*(2-x)
294 <= 1/2*ulp(s) + 5*2^(-w)*(2-x) for w >= 3
295 <= (1/2 + 5 * 2^(-E(s)) * (2-x)) ulp(s) */
296 err_s += 3 - MPFR_GET_EXP(s);
297 err_s = (err_s >= 0) ? err_s + 1 : 0;
298 /* the error on s is bounded by 2^err_s ulp(s), thus by
299 2^(err_s+1)*2^(-w)*|s| since ulp(s) <= 2^(1-w)*|s|.
300 Now n*2^(-w) can always be written |(1+r)^n-1| for some
301 |r|<=2^(-w), thus taking n=2^(err_s+1) we see that
302 |S - s| <= |(1+r)^(2^(err_s+1))-1| * |s|, where S is the
304 In fact if ulp(s) <= ulp(S) the same inequality holds for
305 |S| instead of |s| in the right hand side, i.e., we can
306 write s = (1+r)^(2^(err_s+1)) * S.
307 But if ulp(S) < ulp(s), we need to add one ``bit'' to the error,
308 to get s = (1+r)^(2^(err_s+2)) * S. This is true since with
309 E = n*2^(-w) we have |s - S| <= E * |s|, thus
310 |s - S| <= E/(1-E) * |S|.
311 Now E/(1-E) is bounded by 2E as long as E<=1/2,
312 and 2E can be written (1+r)^(2n)-1 as above.
314 err_s += 2; /* exponent of relative error */
316 mpfr_sub_ui (v, z0, 1, GMP_RNDN); /* v = (x-1) * (1+r) */
317 mpfr_mul (v, v, t, GMP_RNDN); /* v = Pi*(x-1) * (1+r)^3 */
318 mpfr_div (v, v, s, GMP_RNDN); /* Pi*(x-1)/sin(Pi*(2-x)) */
319 mpfr_abs (v, v, GMP_RNDN);
320 /* (1+r)^(3+2^err_s+1) */
321 err_s = (err_s <= 1) ? 3 : err_s + 1;
322 /* now (1+r)^M with M <= 2^err_s */
323 mpfr_log (v, v, GMP_RNDN);
324 /* log(v*(1+e)) = log(v)+log(1+e) where |e| <= 2^(err_s-w).
325 Since |log(1+e)| <= 2*e for |e| <= 1/4, the error on v is
326 bounded by ulp(v)/2 + 2^(err_s+1-w). */
333 err_s += 1 - MPFR_GET_EXP(v);
334 err_s = (err_s >= 0) ? err_s + 1 : 0;
335 /* the error on v is bounded by 2^err_s ulps */
336 err_u += MPFR_GET_EXP(u); /* absolute error on u */
337 err_s += MPFR_GET_EXP(v); /* absolute error on v */
338 mpfr_sub (s, v, u, GMP_RNDN);
339 /* the total error on s is bounded by ulp(s)/2 + 2^(err_u-w)
340 + 2^(err_s-w) <= ulp(s)/2 + 2^(max(err_u,err_s)+1-w) */
341 err_s = (err_s >= err_u) ? err_s : err_u;
342 err_s += 1 - MPFR_GET_EXP(s); /* error is 2^err_s ulp(s) */
343 err_s = (err_s >= 0) ? err_s + 1 : 0;
344 if (mpfr_can_round (s, w - err_s, GMP_RNDN, GMP_RNDZ, precy
345 + (rnd == GMP_RNDN)))
353 MPFR_ASSERTD (compared > 0);
355 /* since k is O(w), the value of log(z0*...*(z0+k-1)) is about w*log(w),
356 so there is a cancellation of ~log(w) in the argument reconstruction */
357 w = precy + MPFR_INT_CEIL_LOG2 (precy);
361 w += MPFR_INT_CEIL_LOG2 (w) + 13;
362 MPFR_ASSERTD (w >= 3);
364 /* argument reduction: we compute gamma(z0 + k), where the series
365 has error term B_{2n}/(z0+k)^(2n) ~ (n/(Pi*e*(z0+k)))^(2n)
366 and we need k steps of argument reconstruction. Assuming k is large
367 with respect to z0, and k = n, we get 1/(Pi*e)^(2n) ~ 2^(-w), i.e.,
368 k ~ w*log(2)/2/log(Pi*e) ~ 0.1616 * w.
369 However, since the series is more expensive to compute, the optimal
370 value seems to be k ~ 4.5 * w experimentally. */
371 mpfr_set_prec (s, 53);
372 mpfr_gamma_alpha (s, w);
373 mpfr_set_ui_2exp (s, 9, -1, GMP_RNDU);
374 mpfr_mul_ui (s, s, w, GMP_RNDU);
375 if (mpfr_cmp (z0, s) < 0)
377 mpfr_sub (s, s, z0, GMP_RNDU);
378 k = mpfr_get_ui (s, GMP_RNDU);
385 mpfr_set_prec (s, w);
386 mpfr_set_prec (t, w);
387 mpfr_set_prec (u, w);
388 mpfr_set_prec (v, w);
389 mpfr_set_prec (z, w);
391 mpfr_add_ui (z, z0, k, GMP_RNDN);
392 /* z = (z0+k)*(1+t1) with |t1| <= 2^(-w) */
394 /* z >= 4 ensures the relative error on log(z) is small,
395 and also (z-1/2)*log(z)-z >= 0 */
396 MPFR_ASSERTD (mpfr_cmp_ui (z, 4) >= 0);
398 mpfr_log (s, z, GMP_RNDN); /* log(z) */
399 /* we have s = log((z0+k)*(1+t1))*(1+t2) with |t1|, |t2| <= 2^(-w).
400 Since w >= 2 and z0+k >= 4, we can write log((z0+k)*(1+t1))
401 = log(z0+k) * (1+t3) with |t3| <= 2^(-w), thus we have
402 s = log(z0+k) * (1+t4)^2 with |t4| <= 2^(-w) */
403 mpfr_mul_2ui (t, z, 1, GMP_RNDN); /* t = 2z * (1+t5) */
404 mpfr_sub_ui (t, t, 1, GMP_RNDN); /* t = 2z-1 * (1+t6)^3 */
405 /* since we can write 2z*(1+t5) = (2z-1)*(1+t5') with
406 t5' = 2z/(2z-1) * t5, thus |t5'| <= 8/7 * t5 */
407 mpfr_mul (s, s, t, GMP_RNDN); /* (2z-1)*log(z) * (1+t7)^6 */
408 mpfr_div_2ui (s, s, 1, GMP_RNDN); /* (z-1/2)*log(z) * (1+t7)^6 */
409 mpfr_sub (s, s, z, GMP_RNDN); /* (z-1/2)*log(z)-z */
410 /* s = [(z-1/2)*log(z)-z]*(1+u)^14, s >= 1/2 */
412 mpfr_ui_div (u, 1, z, GMP_RNDN); /* 1/z * (1+u), u <= 1/4 since z >= 4 */
414 /* the first term is B[2]/2/z = 1/12/z: t=1/12/z, C[2]=1 */
415 mpfr_div_ui (t, u, 12, GMP_RNDN); /* 1/(12z) * (1+u)^2, t <= 3/128 */
416 mpfr_set (v, t, GMP_RNDN); /* (1+u)^2, v < 2^(-5) */
417 mpfr_add (s, s, v, GMP_RNDN); /* (1+u)^15 */
419 mpfr_mul (u, u, u, GMP_RNDN); /* 1/z^2 * (1+u)^3 */
423 B = bernoulli ((mpz_t *) 0, 0);
424 B = bernoulli (B, 1);
428 /* m <= maxm ensures that 2*m*(2*m+1) <= ULONG_MAX */
429 maxm = 1UL << (BITS_PER_MP_LIMB / 2 - 1);
431 /* s:(1+u)^15, t:(1+u)^2, t <= 3/128 */
433 for (m = 2; MPFR_GET_EXP(v) + (mp_exp_t) w >= MPFR_GET_EXP(s); m++)
435 mpfr_mul (t, t, u, GMP_RNDN); /* (1+u)^(10m-14) */
438 mpfr_mul_ui (t, t, 2*(m-1)*(2*m-3), GMP_RNDN);
439 mpfr_div_ui (t, t, 2*m*(2*m-1), GMP_RNDN);
440 mpfr_div_ui (t, t, 2*m*(2*m+1), GMP_RNDN);
444 mpfr_mul_ui (t, t, 2*(m-1), GMP_RNDN);
445 mpfr_mul_ui (t, t, 2*m-3, GMP_RNDN);
446 mpfr_div_ui (t, t, 2*m, GMP_RNDN);
447 mpfr_div_ui (t, t, 2*m-1, GMP_RNDN);
448 mpfr_div_ui (t, t, 2*m, GMP_RNDN);
449 mpfr_div_ui (t, t, 2*m+1, GMP_RNDN);
452 /* invariant: t=1/(2m)/(2m-1)/z^(2m-1)/(2m+1)! */
455 B = bernoulli (B, m); /* B[2m]*(2m+1)!, exact */
458 mpfr_mul_z (v, t, B[m], GMP_RNDN); /* (1+u)^(10m-7) */
459 MPFR_ASSERTD(MPFR_GET_EXP(v) <= - (2 * m + 3));
460 mpfr_add (s, s, v, GMP_RNDN);
462 /* m <= 1/2*Pi*e*z ensures that |v[m]| < 1/2^(2m+3) */
463 MPFR_ASSERTD ((double) m <= 4.26 * mpfr_get_d (z, GMP_RNDZ));
465 /* We have sum([(1+u)^(10m-7)-1]*1/2^(2m+3), m=2..infinity)
466 <= 1.46*u for u <= 2^(-3).
467 We have 0 < lngamma(z) - [(z - 1/2) ln(z) - z + 1/2 ln(2 Pi)] < 0.021
468 for z >= 4, thus since the initial s >= 0.85, the different values of
469 s differ by at most one binade, and the total rounding error on s
470 in the for-loop is bounded by 2*(m-1)*ulp(final_s).
471 The error coming from the v's is bounded by
472 1.46*2^(-w) <= 2*ulp(final_s).
473 Thus the total error so far is bounded by [(1+u)^15-1]*s+2m*ulp(s)
475 Taking into account the truncation error (which is bounded by the last
476 term v[] according to 6.1.42 in A&S), the bound is (2m+48)*ulp(s).
479 /* add 1/2*log(2*Pi) and subtract log(z0*(z0+1)*...*(z0+k-1)) */
480 mpfr_const_pi (v, GMP_RNDN); /* v = Pi*(1+u) */
481 mpfr_mul_2ui (v, v, 1, GMP_RNDN); /* v = 2*Pi * (1+u) */
485 mpfr_set (t, z0, GMP_RNDN); /* t = z0*(1+u) */
486 for (l = 1; l < k; l++)
488 mpfr_add_ui (u, z0, l, GMP_RNDN); /* u = (z0+l)*(1+u) */
489 mpfr_mul (t, t, u, GMP_RNDN); /* (1+u)^(2l+1) */
491 /* now t: (1+u)^(2k-1) */
492 /* instead of computing log(sqrt(2*Pi)/t), we compute
493 1/2*log(2*Pi/t^2), which trades a square root for a square */
494 mpfr_mul (t, t, t, GMP_RNDN); /* (z0*...*(z0+k-1))^2, (1+u)^(4k-1) */
495 mpfr_div (v, v, t, GMP_RNDN);
496 /* 2*Pi/(z0*...*(z0+k-1))^2 (1+u)^(4k+1) */
499 err_s = MPFR_GET_EXP(s);
500 mpfr_exp (s, s, GMP_RNDN);
501 /* before the exponential, we have s = s0 + h where
502 |h| <= (2m+48)*ulp(s), thus exp(s0) = exp(s) * exp(-h).
503 For |h| <= 1/4, we have |exp(h)-1| <= 1.2*|h| thus
504 |exp(s) - exp(s0)| <= 1.2 * exp(s) * (2m+48)* 2^(EXP(s)-w). */
505 d = 1.2 * (2.0 * (double) m + 48.0);
506 /* the error on s is bounded by d*2^err_s * 2^(-w) */
507 mpfr_sqrt (t, v, GMP_RNDN);
508 /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1),
509 thus t = sqrt(v0)*(1+u)^(2k+3/2). */
510 mpfr_mul (s, s, t, GMP_RNDN);
511 /* the error on input s is bounded by (1+u)^(d*2^err_s),
512 and that on t is (1+u)^(2k+3/2), thus the
513 total error is (1+u)^(d*2^err_s+2k+5/2) */
514 err_s += __gmpfr_ceil_log2 (d);
515 err_t = __gmpfr_ceil_log2 (2.0 * (double) k + 2.5);
516 err_s = (err_s >= err_t) ? err_s + 1 : err_t + 1;
518 mpfr_log (t, v, GMP_RNDN);
519 /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1),
520 thus log(v) = log(v0) + (4k+1)*log(1+u). Since |log(1+u)/u| <= 1.07
521 for |u| <= 2^(-3), the absolute error on log(v) is bounded by
522 1.07*(4k+1)*u, and the rounding error by ulp(t). */
523 mpfr_div_2ui (t, t, 1, GMP_RNDN);
524 /* the error on t is now bounded by ulp(t) + 0.54*(4k+1)*2^(-w).
525 We have sqrt(2*Pi)/(z0*(z0+1)*...*(z0+k-1)) <= sqrt(2*Pi)/k! <= 0.5
526 since k>=3, thus t <= -0.5 and ulp(t) >= 2^(-w).
527 Thus the error on t is bounded by (2.16*k+1.54)*ulp(t). */
528 err_t = MPFR_GET_EXP(t) + (mp_exp_t)
529 __gmpfr_ceil_log2 (2.2 * (double) k + 1.6);
530 err_s = MPFR_GET_EXP(s) + (mp_exp_t)
531 __gmpfr_ceil_log2 (2.0 * (double) m + 48.0);
532 mpfr_add (s, s, t, GMP_RNDN); /* this is a subtraction in fact */
533 /* the final error in ulp(s) is
534 <= 1 + 2^(err_t-EXP(s)) + 2^(err_s-EXP(s))
535 <= 2^(1+max(err_t,err_s)-EXP(s)) if err_t <> err_s
536 <= 2^(2+max(err_t,err_s)-EXP(s)) if err_t = err_s */
537 err_s = (err_t == err_s) ? 1 + err_s : ((err_t > err_s) ? err_t : err_s);
538 err_s += 1 - MPFR_GET_EXP(s);
541 while (MPFR_UNLIKELY (!MPFR_CAN_ROUND (s, w - err_s, precy, rnd)));
546 (*__gmp_free_func) (B, oldBm * sizeof (mpz_t));
549 inexact = mpfr_set (y, s, rnd);
557 MPFR_SAVE_EXPO_FREE (expo);
558 return mpfr_check_range (y, inexact, rnd);
564 mpfr_lngamma (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
568 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
569 ("lngamma[%#R]=%R inexact=%d", y, y, inex));
572 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
574 if (MPFR_IS_NAN (x) || MPFR_IS_NEG (x))
579 else /* lngamma(+Inf) = lngamma(+0) = +Inf */
583 MPFR_RET (0); /* exact */
587 /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */
588 if (MPFR_IS_NEG (x) && (unit_bit (x) == 0 || mpfr_integer_p (x)))
594 inex = mpfr_lngamma_aux (y, x, rnd);
599 mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mp_rnd_t rnd)
603 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
604 ("lgamma[%#R]=%R inexact=%d", y, y, inex));
606 *signp = 1; /* most common case */
608 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
617 *signp = MPFR_INT_SIGN (x);
626 if (mpfr_integer_p (x))
633 if (unit_bit (x) == 0)
636 /* For tiny negative x, we have gamma(x) = 1/x - euler + O(x),
637 thus |gamma(x)| = -1/x + euler + O(x), and
638 log |gamma(x)| = -log(-x) - euler*x + O(x^2).
639 More precisely we have for -0.4 <= x < 0:
640 -log(-x) <= log |gamma(x)| <= -log(-x) - x.
641 Since log(x) is not representable, we may have an instance of the
642 Table Maker Dilemma. The only way to ensure correct rounding is to
643 compute an interval [l,h] such that l <= -log(-x) and
644 -log(-x) - x <= h, and check whether l and h round to the same number
645 for the target precision and rounding modes. */
646 if (MPFR_EXP(x) + 1 <= - (mp_exp_t) MPFR_PREC(y))
647 /* since PREC(y) >= 1, this ensures EXP(x) <= -2,
648 thus |x| <= 0.25 < 0.4 */
652 mp_prec_t w = MPFR_PREC (y) + 14;
658 /* we want a lower bound on -log(-x), thus an upper bound
659 on log(-x), thus an upper bound on -x. */
660 mpfr_neg (l, x, GMP_RNDU); /* upper bound on -x */
661 mpfr_log (l, l, GMP_RNDU); /* upper bound for log(-x) */
662 mpfr_neg (l, l, GMP_RNDD); /* lower bound for -log(-x) */
663 mpfr_neg (h, x, GMP_RNDD); /* lower bound on -x */
664 mpfr_log (h, h, GMP_RNDD); /* lower bound on log(-x) */
665 mpfr_neg (h, h, GMP_RNDU); /* upper bound for -log(-x) */
666 mpfr_sub (h, h, x, GMP_RNDU); /* upper bound for -log(-x) - x */
667 inex = mpfr_prec_round (l, MPFR_PREC (y), rnd);
668 inex2 = mpfr_prec_round (h, MPFR_PREC (y), rnd);
669 /* Caution: we not only need l = h, but both inexact flags
670 should agree. Indeed, one of the inexact flags might be
671 zero. In that case if we assume ln|gamma(x)| cannot be
672 exact, the other flag should be correct. We are conservative
673 here and request that both inexact flags agree. */
674 ok = SAME_SIGN (inex, inex2) && mpfr_equal_p (l, h);
676 mpfr_set (y, h, rnd); /* exact */
681 /* if ulp(log(-x)) <= |x| there is no reason to loop,
682 since the width of [l, h] will be at least |x| */
683 if (MPFR_EXP(l) < MPFR_EXP(x) + (mp_exp_t) w)
685 w += MPFR_INT_CEIL_LOG2(w) + 3;
690 inex = mpfr_lngamma_aux (y, x, rnd);