contrib/mpfr/lngamma.c

   1 /* mpfr_lngamma -- lngamma function
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   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.
  12
  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.
  17
  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. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* assuming b[0]...b[2(n-1)] are computed, computes and stores B[2n]*(2n+1)!
  27
  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)
  32    which gives:
  33    B[n] = -sum(binomial(n+1,k)*B[k], k=0..n-1)/(n+1).
  34
  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.
  38 */
  39 static mpz_t*
  40 bernoulli (mpz_t *b, unsigned long n)
  41 {
  42   if (n == 0)
  43     {
  44       b = (mpz_t *) (*__gmp_allocate_func) (sizeof (mpz_t));
  45       mpz_init_set_ui (b[0], 1);
  46     }
  47   else
  48     {
  49       mpz_t t;
  50       unsigned long k;
  51
  52       b = (mpz_t *) (*__gmp_reallocate_func)
  53         (b, n * sizeof (mpz_t), (n + 1) * sizeof (mpz_t));
  54       mpz_init (b[n]);
  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);
  59       mpz_mul_ui (t, t, n);
  60       mpz_div_ui (t, t, 3); /* exact: t=binomial(2*n+1,2*k)*(2*n)!/(2*k+1)!
  61                                for k=n-1 */
  62       mpz_mul (b[n], t, b[n-1]);
  63       for (k = n - 1; k-- > 0;)
  64         {
  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]);
  72         }
  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);
  77       mpz_neg (b[n], b[n]);
  78       mpz_clear (t);
  79     }
  80   return b;
  81 }
  82
  83 /* given a precision p, return alpha, such that the argument reduction
  84    will use k = alpha*p*log(2).
  85
  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.
  89 */
  90 static void
  91 mpfr_gamma_alpha (mpfr_t s, mp_prec_t p)
  92 {
  93   if (p <= 100)
  94     mpfr_set_ui_2exp (s, 614, -10, GMP_RNDN); /* about 0.6 */
  95   else if (p <= 500)
  96     mpfr_set_ui_2exp (s, 819, -10, GMP_RNDN); /* about 0.8 */
  97   else if (p <= 1000)
  98     mpfr_set_ui_2exp (s, 1331, -10, GMP_RNDN); /* about 1.3 */
  99   else if (p <= 2000)
 100     mpfr_set_ui_2exp (s, 1741, -10, GMP_RNDN); /* about 1.7 */
 101   else if (p <= 5000)
 102     mpfr_set_ui_2exp (s, 2253, -10, GMP_RNDN); /* about 2.2 */
 103   else if (p <= 10000)
 104     mpfr_set_ui_2exp (s, 3482, -10, GMP_RNDN); /* about 3.4 */
 105   else
 106     mpfr_set_ui_2exp (s, 9, -1, GMP_RNDN); /* 4.5 */
 107 }
 108
 109 #ifndef IS_GAMMA
 110 static int
 111 unit_bit (mpfr_srcptr (x))
 112 {
 113   mp_exp_t expo;
 114   mp_prec_t prec;
 115   mp_limb_t x0;
 116
 117   expo = MPFR_GET_EXP (x);
 118   if (expo <= 0)
 119     return 0;  /* |x| < 1 */
 120
 121   prec = MPFR_PREC (x);
 122   if (expo > prec)
 123     return 0;  /* y is a multiple of 2^(expo-prec), thus an even integer */
 124
 125   /* Now, the unit bit is represented. */
 126
 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) */
 129
 130   x0 = *(MPFR_MANT (x) + prec / BITS_PER_MP_LIMB);
 131   /* limb containing the unit bit */
 132
 133   return (x0 >> (prec % BITS_PER_MP_LIMB)) & 1;
 134 }
 135 #endif
 136
 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.
 145  */
 146 #ifdef IS_GAMMA
 147 #define GAMMA_FUNC mpfr_gamma_aux
 148 #else
 149 #define GAMMA_FUNC mpfr_lngamma_aux
 150 #endif
 151
 152 static int
 153 GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd)
 154 {
 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[] */
 162   unsigned long oldBm;
 163   double d;
 164   MPFR_SAVE_EXPO_DECL (expo);
 165
 166   compared = mpfr_cmp_ui (z0, 1);
 167
 168   MPFR_SAVE_EXPO_MARK (expo);
 169
 170 #ifndef IS_GAMMA /* lngamma or lgamma */
 171   if (compared == 0 || (compared > 0 && mpfr_cmp_ui (z0, 2) == 0))
 172     {
 173       MPFR_SAVE_EXPO_FREE (expo);
 174       return mpfr_set_ui (y, 0, GMP_RNDN);  /* lngamma(1 or 2) = +0 */
 175     }
 176
 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))
 180     {
 181       mpfr_t l, h, g;
 182       int ok, inex2;
 183       mp_prec_t prec = MPFR_PREC(y) + 14;
 184       MPFR_ZIV_DECL (loop);
 185
 186       MPFR_ZIV_INIT (loop, prec);
 187       do
 188         {
 189           mpfr_init2 (l, prec);
 190           if (MPFR_IS_POS(z0))
 191             {
 192               mpfr_log (l, z0, GMP_RNDU); /* upper bound for log(z0) */
 193               mpfr_init2 (h, MPFR_PREC(l));
 194             }
 195           else
 196             {
 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));
 201             }
 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
 205                                  to mpfr_log */
 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;
 225           if (ok)
 226             mpfr_set (y, h, rnd); /* exact */
 227           mpfr_clear (l);
 228           mpfr_clear (h);
 229           mpfr_clear (g);
 230           if (ok)
 231             {
 232               MPFR_SAVE_EXPO_FREE (expo);
 233               return mpfr_check_range (y, inexact, rnd);
 234             }
 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);
 242         }
 243       while (prec <= -MPFR_EXP(z0));
 244       MPFR_ZIV_FREE (loop);
 245     }
 246 #endif
 247
 248   precy = MPFR_PREC(y);
 249
 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);
 255
 256   if (compared < 0)
 257     {
 258       mp_exp_t err_u;
 259
 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) */
 263
 264       w = precy + MPFR_INT_CEIL_LOG2 (precy);
 265       while (1)
 266         {
 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
 275              values. */
 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 */
 289
 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
 303              true value.
 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.
 313           */
 314           err_s += 2; /* exponent of relative error */
 315
 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). */
 327           if (err_s + 2 > w)
 328             {
 329               w += err_s + 2;
 330             }
 331           else
 332             {
 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)))
 346                 goto end;
 347             }
 348         }
 349     }
 350
 351   /* now z0 > 1 */
 352
 353   MPFR_ASSERTD (compared > 0);
 354
 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);
 358
 359   do
 360     {
 361       w += MPFR_INT_CEIL_LOG2 (w) + 13;
 362       MPFR_ASSERTD (w >= 3);
 363
 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)
 376         {
 377           mpfr_sub (s, s, z0, GMP_RNDU);
 378           k = mpfr_get_ui (s, GMP_RNDU);
 379           if (k < 3)
 380             k = 3;
 381         }
 382       else
 383         k = 3;
 384
 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);
 390
 391       mpfr_add_ui (z, z0, k, GMP_RNDN);
 392       /* z = (z0+k)*(1+t1) with |t1| <= 2^(-w) */
 393
 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);
 397
 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 */
 411
 412       mpfr_ui_div (u, 1, z, GMP_RNDN); /* 1/z * (1+u), u <= 1/4 since z >= 4 */
 413
 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 */
 418
 419       mpfr_mul (u, u, u, GMP_RNDN); /* 1/z^2 * (1+u)^3 */
 420
 421       if (Bm == 0)
 422         {
 423           B = bernoulli ((mpz_t *) 0, 0);
 424           B = bernoulli (B, 1);
 425           Bm = 2;
 426         }
 427
 428       /* m <= maxm ensures that 2*m*(2*m+1) <= ULONG_MAX */
 429       maxm = 1UL << (BITS_PER_MP_LIMB / 2 - 1);
 430
 431       /* s:(1+u)^15, t:(1+u)^2, t <= 3/128 */
 432
 433       for (m = 2; MPFR_GET_EXP(v) + (mp_exp_t) w >= MPFR_GET_EXP(s); m++)
 434         {
 435           mpfr_mul (t, t, u, GMP_RNDN); /* (1+u)^(10m-14) */
 436           if (m <= maxm)
 437             {
 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);
 441             }
 442           else
 443             {
 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);
 450             }
 451           /* (1+u)^(10m-8) */
 452           /* invariant: t=1/(2m)/(2m-1)/z^(2m-1)/(2m+1)! */
 453           if (Bm <= m)
 454             {
 455               B = bernoulli (B, m); /* B[2m]*(2m+1)!, exact */
 456               Bm ++;
 457             }
 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);
 461         }
 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));
 464
 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)
 474          <= (2m+47)*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).
 477       */
 478
 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) */
 482       if (k)
 483         {
 484           unsigned long l;
 485           mpfr_set (t, z0, GMP_RNDN); /* t = z0*(1+u) */
 486           for (l = 1; l < k; l++)
 487             {
 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) */
 490             }
 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) */
 497         }
 498 #ifdef IS_GAMMA
 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;
 517 #else
 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);
 539 #endif
 540     }
 541   while (MPFR_UNLIKELY (!MPFR_CAN_ROUND (s, w - err_s, precy, rnd)));
 542
 543   oldBm = Bm;
 544   while (Bm--)
 545     mpz_clear (B[Bm]);
 546   (*__gmp_free_func) (B, oldBm * sizeof (mpz_t));
 547
 548  end:
 549   inexact = mpfr_set (y, s, rnd);
 550
 551   mpfr_clear (s);
 552   mpfr_clear (t);
 553   mpfr_clear (u);
 554   mpfr_clear (v);
 555   mpfr_clear (z);
 556
 557   MPFR_SAVE_EXPO_FREE (expo);
 558   return mpfr_check_range (y, inexact, rnd);
 559 }
 560
 561 #ifndef IS_GAMMA
 562
 563 int
 564 mpfr_lngamma (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
 565 {
 566   int inex;
 567
 568   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
 569                  ("lngamma[%#R]=%R inexact=%d", y, y, inex));
 570
 571   /* special cases */
 572   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 573     {
 574       if (MPFR_IS_NAN (x) || MPFR_IS_NEG (x))
 575         {
 576           MPFR_SET_NAN (y);
 577           MPFR_RET_NAN;
 578         }
 579       else /* lngamma(+Inf) = lngamma(+0) = +Inf */
 580         {
 581           MPFR_SET_INF (y);
 582           MPFR_SET_POS (y);
 583           MPFR_RET (0);  /* exact */
 584         }
 585     }
 586
 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)))
 589     {
 590       MPFR_SET_NAN (y);
 591       MPFR_RET_NAN;
 592     }
 593
 594   inex = mpfr_lngamma_aux (y, x, rnd);
 595   return inex;
 596 }
 597
 598 int
 599 mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mp_rnd_t rnd)
 600 {
 601   int inex;
 602
 603   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
 604                  ("lgamma[%#R]=%R inexact=%d", y, y, inex));
 605
 606   *signp = 1;  /* most common case */
 607
 608   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 609     {
 610       if (MPFR_IS_NAN (x))
 611         {
 612           MPFR_SET_NAN (y);
 613           MPFR_RET_NAN;
 614         }
 615       else
 616         {
 617           *signp = MPFR_INT_SIGN (x);
 618           MPFR_SET_INF (y);
 619           MPFR_SET_POS (y);
 620           MPFR_RET (0);
 621         }
 622     }
 623
 624   if (MPFR_IS_NEG (x))
 625     {
 626       if (mpfr_integer_p (x))
 627         {
 628           MPFR_SET_INF (y);
 629           MPFR_SET_POS (y);
 630           MPFR_RET (0);
 631         }
 632
 633       if (unit_bit (x) == 0)
 634         *signp = -1;
 635
 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 */
 649         {
 650           mpfr_t l, h;
 651           int ok, inex2;
 652           mp_prec_t w = MPFR_PREC (y) + 14;
 653
 654           while (1)
 655             {
 656               mpfr_init2 (l, w);
 657               mpfr_init2 (h, w);
 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);
 675               if (ok)
 676                 mpfr_set (y, h, rnd); /* exact */
 677               mpfr_clear (l);
 678               mpfr_clear (h);
 679               if (ok)
 680                 return inex;
 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)
 684                 break;
 685               w += MPFR_INT_CEIL_LOG2(w) + 3;
 686             }
 687         }
 688     }
 689
 690   inex = mpfr_lngamma_aux (y, x, rnd);
 691   return inex;
 692 }
 693
 694 #endif