source/libs/mpfr/mpfr-src/src/const_euler.c

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