contrib/mpfr/src/const_euler.c

   1 /* mpfr_const_euler -- Euler's constant
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel 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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* Declare the cache */
  27 MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_euler, mpfr_const_euler_internal);
  28
  29 /* Set User Interface */
  30 #undef mpfr_const_euler
  31 int
  32 mpfr_const_euler (mpfr_ptr x, mpfr_rnd_t rnd_mode) {
  33   return mpfr_cache (x, __gmpfr_cache_const_euler, rnd_mode);
  34 }
  35
  36
  37 static void mpfr_const_euler_S2 (mpfr_ptr, unsigned long);
  38 static void mpfr_const_euler_R (mpfr_ptr, unsigned long);
  39
  40 int
  41 mpfr_const_euler_internal (mpfr_t x, mpfr_rnd_t rnd)
  42 {
  43   mpfr_prec_t prec = MPFR_PREC(x), m, log2m;
  44   mpfr_t y, z;
  45   unsigned long n;
  46   int inexact;
  47   MPFR_ZIV_DECL (loop);
  48
  49   log2m = MPFR_INT_CEIL_LOG2 (prec);
  50   m = prec + 2 * log2m + 23;
  51
  52   mpfr_init2 (y, m);
  53   mpfr_init2 (z, m);
  54
  55   MPFR_ZIV_INIT (loop, m);
  56   for (;;)
  57     {
  58       mpfr_exp_t exp_S, err;
  59       /* since prec >= 1, we have m >= 24 here, which ensures n >= 9 below */
  60       n = 1 + (unsigned long) ((double) m * LOG2 / 2.0);
  61       MPFR_ASSERTD (n >= 9);
  62       mpfr_const_euler_S2 (y, n); /* error <= 3 ulps */
  63       exp_S = MPFR_EXP(y);
  64       mpfr_set_ui (z, n, MPFR_RNDN);
  65       mpfr_log (z, z, MPFR_RNDD); /* error <= 1 ulp */
  66       mpfr_sub (y, y, z, MPFR_RNDN); /* S'(n) - log(n) */
  67       /* the error is less than 1/2 + 3*2^(exp_S-EXP(y)) + 2^(EXP(z)-EXP(y))
  68          <= 1/2 + 2^(exp_S+2-EXP(y)) + 2^(EXP(z)-EXP(y))
  69          <= 1/2 + 2^(1+MAX(exp_S+2,EXP(z))-EXP(y)) */
  70       err = 1 + MAX(exp_S + 2, MPFR_EXP(z)) - MPFR_EXP(y);
  71       err = (err >= -1) ? err + 1 : 0; /* error <= 2^err ulp(y) */
  72       exp_S = MPFR_EXP(y);
  73       mpfr_const_euler_R (z, n); /* err <= ulp(1/2) = 2^(-m) */
  74       mpfr_sub (y, y, z, MPFR_RNDN);
  75       /* err <= 1/2 ulp(y) + 2^(-m) + 2^(err + exp_S - EXP(y)) ulp(y).
  76          Since the result is between 0.5 and 1, ulp(y) = 2^(-m).
  77          So we get 3/2*ulp(y) + 2^(err + exp_S - EXP(y)) ulp(y).
  78          3/2 + 2^e <= 2^(e+1) for e>=1, and <= 2^2 otherwise */
  79       err = err + exp_S - MPFR_EXP(y);
  80       err = (err >= 1) ? err + 1 : 2;
  81       if (MPFR_LIKELY (MPFR_CAN_ROUND (y, m - err, prec, rnd)))
  82         break;
  83       MPFR_ZIV_NEXT (loop, m);
  84       mpfr_set_prec (y, m);
  85       mpfr_set_prec (z, m);
  86     }
  87   MPFR_ZIV_FREE (loop);
  88
  89   inexact = mpfr_set (x, y, rnd);
  90
  91   mpfr_clear (y);
  92   mpfr_clear (z);
  93
  94   return inexact; /* always inexact */
  95 }
  96
  97 static void
  98 mpfr_const_euler_S2_aux (mpz_t P, mpz_t Q, mpz_t T, unsigned long n,
  99                          unsigned long a, unsigned long b, int need_P)
 100 {
 101   if (a + 1 == b)
 102     {
 103       mpz_set_ui (P, n);
 104       if (a > 1)
 105         mpz_mul_si (P, P, 1 - (long) a);
 106       mpz_set (T, P);
 107       mpz_set_ui (Q, a);
 108       mpz_mul_ui (Q, Q, a);
 109     }
 110   else
 111     {
 112       unsigned long c = (a + b) / 2;
 113       mpz_t P2, Q2, T2;
 114       mpfr_const_euler_S2_aux (P, Q, T, n, a, c, 1);
 115       mpz_init (P2);
 116       mpz_init (Q2);
 117       mpz_init (T2);
 118       mpfr_const_euler_S2_aux (P2, Q2, T2, n, c, b, 1);
 119       mpz_mul (T, T, Q2);
 120       mpz_mul (T2, T2, P);
 121       mpz_add (T, T, T2);
 122       if (need_P)
 123         mpz_mul (P, P, P2);
 124       mpz_mul (Q, Q, Q2);
 125       mpz_clear (P2);
 126       mpz_clear (Q2);
 127       mpz_clear (T2);
 128       /* divide by 2 if possible */
 129       {
 130         unsigned long v2;
 131         v2 = mpz_scan1 (P, 0);
 132         c = mpz_scan1 (Q, 0);
 133         if (c < v2)
 134           v2 = c;
 135         c = mpz_scan1 (T, 0);
 136         if (c < v2)
 137           v2 = c;
 138         if (v2)
 139           {
 140             mpz_tdiv_q_2exp (P, P, v2);
 141             mpz_tdiv_q_2exp (Q, Q, v2);
 142             mpz_tdiv_q_2exp (T, T, v2);
 143           }
 144       }
 145     }
 146 }
 147
 148 /* computes S(n) = sum(n^k*(-1)^(k-1)/k!/k, k=1..ceil(4.319136566 * n))
 149    using binary splitting.
 150    We have S(n) = sum(f(k), k=1..N) with N=ceil(4.319136566 * n)
 151    and f(k) = n^k*(-1)*(k-1)/k!/k,
 152    thus f(k)/f(k-1) = -n*(k-1)/k^2
 153 */
 154 static void
 155 mpfr_const_euler_S2 (mpfr_t x, unsigned long n)
 156 {
 157   mpz_t P, Q, T;
 158   unsigned long N = (unsigned long) (ALPHA * (double) n + 1.0);
 159   mpz_init (P);
 160   mpz_init (Q);
 161   mpz_init (T);
 162   mpfr_const_euler_S2_aux (P, Q, T, n, 1, N + 1, 0);
 163   mpfr_set_z (x, T, MPFR_RNDN);
 164   mpfr_div_z (x, x, Q, MPFR_RNDN);
 165   mpz_clear (P);
 166   mpz_clear (Q);
 167   mpz_clear (T);
 168 }
 169
 170 /* computes R(n) = exp(-n)/n * sum(k!/(-n)^k, k=0..n-2)
 171    with error at most 4*ulp(x). Assumes n>=2.
 172    Since x <= exp(-n)/n <= 1/8, then 4*ulp(x) <= ulp(1).
 173 */
 174 static void
 175 mpfr_const_euler_R (mpfr_t x, unsigned long n)
 176 {
 177   unsigned long k, m;
 178   mpz_t a, s;
 179   mpfr_t y;
 180
 181   MPFR_ASSERTN (n >= 2); /* ensures sum(k!/(-n)^k, k=0..n-2) >= 2/3 */
 182
 183   /* as we multiply the sum by exp(-n), we need only PREC(x) - n/LOG2 bits */
 184   m = MPFR_PREC(x) - (unsigned long) ((double) n / LOG2);
 185
 186   mpz_init_set_ui (a, 1);
 187   mpz_mul_2exp (a, a, m);
 188   mpz_init_set (s, a);
 189
 190   for (k = 1; k <= n; k++)
 191     {
 192       mpz_mul_ui (a, a, k);
 193       mpz_fdiv_q_ui (a, a, n);
 194       /* the error e(k) on a is e(k) <= 1 + k/n*e(k-1) with e(0)=0,
 195          i.e. e(k) <= k */
 196       if (k % 2)
 197         mpz_sub (s, s, a);
 198       else
 199         mpz_add (s, s, a);
 200     }
 201   /* the error on s is at most 1+2+...+n = n*(n+1)/2 */
 202   mpz_fdiv_q_ui (s, s, n); /* err <= 1 + (n+1)/2 */
 203   MPFR_ASSERTN (MPFR_PREC(x) >= mpz_sizeinbase(s, 2));
 204   mpfr_set_z (x, s, MPFR_RNDD); /* exact */
 205   mpfr_div_2ui (x, x, m, MPFR_RNDD);
 206   /* now x = 1/n * sum(k!/(-n)^k, k=0..n-2) <= 1/n */
 207   /* err(x) <= (n+1)/2^m <= (n+1)*exp(n)/2^PREC(x) */
 208
 209   mpfr_init2 (y, m);
 210   mpfr_set_si (y, -(long)n, MPFR_RNDD); /* assumed exact */
 211   mpfr_exp (y, y, MPFR_RNDD); /* err <= ulp(y) <= exp(-n)*2^(1-m) */
 212   mpfr_mul (x, x, y, MPFR_RNDD);
 213   /* err <= ulp(x) + (n + 1 + 2/n) / 2^prec(x)
 214      <= ulp(x) + (n + 1 + 2/n) ulp(x)/x since x*2^(-prec(x)) < ulp(x)
 215      <= ulp(x) + (n + 1 + 2/n) 3/(2n) ulp(x) since x >= 2/3*n for n >= 2
 216      <= 4 * ulp(x) for n >= 2 */
 217   mpfr_clear (y);
 218
 219   mpz_clear (a);
 220   mpz_clear (s);
 221 }