ex-2/fanc{y,ier}: implement proper computation of N

ex-2/fancier.c

@@ -2,27 +2,42 @@
 #include <stdlib.h>
 #include <stdarg.h>
 #include <math.h>
+#include <gsl/gsl_sf_lambert.h>
 #include <gmp.h>
 #include <time.h>
 
 /* The Euler-Mascheroni constant is here computed to
  * arbitrary precision using GMP rational numbers
- * with the formula:
+ * and the standard Brent-McMillan algorithm.
  *
  *    γ = A(N)/B(N) - C(N)/B(N)² - log(N)
  *
- * with:
+ * where:
  *
  *    A(N) = Σ_(k = 1)^(k_max) (N^k/k!)² * H(k)
  *    B(N) = Σ_(k = 0)^(k_max) (N^k/k!)²
  *    C(N) = 1/4N Σ_(k = 0)^(2N) ((2k)!)^3/((k!)^4*(16N))^2k
  *    H(k) = Σ_(j = 1)^(k) (1/k), the k-th harmonic number
  *
- * where N is computed from D as written below and k_max is
+ * The error of the estimation is given by
- * currently 5*N, chosen in a completely arbitrary way.
  *
- * source: http://www.numberworld.org/y-cruncher/internals/formulas.html
+ *   Δ(N) = 5/12 √(2π) exp(-8N)/√x
  *
+ * so, given a number D of decimal digits to compute we ask
+ * that the error be smaller than 10^-D. The smallest integer
+ * to solve the inequality can be proven to be
+ *
+ *   N = floor(W(5π/9 10^(2D + 1))/16) + 1
+ *
+ *  where W is the principal value of the Lambert W function.
+ *
+ * The number of series terms needed to reduce the error, due to
+ * to less than 10^-D is α⋅2^p where α is the solution of:
+ *
+ *   α log(α) = 3 + α  ⇒  α = exp(W(3e) - 1) ≈ 4.97
+ *
+ * source: Precise error estimate of the Brent-McMillan algorithm for
+ * the computation of Euler's constant, Jean-Pierre Demailly.
  */
 
 /* `mpq_dec_str(z, n)`
@@ -490,6 +505,11 @@ int main(int argc, char** argv) {
     return EXIT_FAILURE;
   }
 
+  // requested decimal digits
+  size_t D = atol(argv[1]);
+  // Brent-McMillan number N
+  size_t N = gsl_sf_lambert_W0(5*M_PI/9*pow(10, 2*D + 1))/16 + 1;
+
   mpq_t exact, indicator;
   mpq_inits(exact, indicator, NULL);
 
@@ -525,11 +545,6 @@ int main(int argc, char** argv) {
       "137174210/1111111111",
       10);
 
-  // number of decimal places
-  size_t D = atol(argv[1]);
-  // number of terms
-  size_t N = 10 + floor(1.0/8 * log(10) * (double)D);
-
   fprintf(stderr, "[main] N: %ld\n", N);
 
   /* calculate series */

ex-2/fancy.c

@@ -1,24 +1,40 @@
 #include <stdio.h>
-#include <math.h>
 #include <stdlib.h>
+#include <math.h>
+#include <gsl/gsl_sf_lambert.h>
 
-// The Euler-Mascheroni constant is computed through the formula:
+/* The Euler-Mascheroni constant is computed to D decimal
-//
+ * places using the refined Brent-McMillan formula:
-//    γ = A(N)/B(N) - ln(N)
+ *
-//
+ *    γ = A(N)/B(N) - C(N)/B(N)² - log(N)
-// with:
+ *
-//
+ * where:
-//    A(N) = Σ_(k = 1)^(k_max) (N^k/k!) * H(k)
+ *
-//    B(N) = Σ_(k = 0)^(k_max) (N^k/k!)
+ *    A(N) = Σ_(k = 1)^(k_max) (N^k/k!)² * H(k)
-//    H(k) = Σ_(j = 1)^(k) (1/k)
+ *    B(N) = Σ_(k = 0)^(k_max) (N^k/k!)²
-//
+ *    C(N) = 1/4N Σ_(k = 0)^(2N) ((2k)!)^3/((k!)^4*(16N))^2k
-// where N is computed from D as written below and k_max is the value
+ *    H(k) = Σ_(j = 1)^(k) (1/k), the k-th harmonic number
-// at which there is no difference between two consecutive terms of
+ *
-// the sum because of double precision.
+ * The error of the estimation is given by
-//
+ *
-// source: http://www.numberworld.org/y-cruncher/internals/formulas.html
+ *   Δ(N) = 5/12 √(2π) exp(-8N)/√x
+ *
+ * so, given a number D of decimal digits to compute we ask
+ * that the error be smaller than 10^-D. The smallest integer
+ * to solve the inequality can be proven to be
+ *
+ *   N = floor(W(5π/9 10^(2D + 1))/16) + 1
+ *
+ *  where W is the principal value of the Lambert W function.
+ *
+ * The series are truncated at k_max, which is when the difference
+ * between two consecutive terms of the sum is zero, at double precision.
+ *
+ * source: Precise error estimate of the Brent-McMillan algorithm for
+ * the computation of Euler's constant, Jean-Pierre Demailly.
+ */
 
-// Partial harmonic sum h
+// Partial harmonic sum H
 double harmonic_sum(double n) {
   double sum = 0;
   for (double k = 1; k < n+1; k++) {
@@ -49,6 +65,7 @@ double b_series(int N){
   return sum;
 }
 
+// C series
 double c_series(int N) {
   double sum = 0;
   for (double k = 0; k < N; k++) {
@@ -57,23 +74,28 @@ double c_series(int N) {
   return sum/(4.0*N);
 }
 
-// Takes in input the number D of desired correct decimals
+/* The Best result is obtained with D=15, which accurately
-// Best result obtained with D = 15, N = 10
+ * computes 15 digits and gives an error of 3.3e-16.
-
+ *
+ * Increasing to D=21 decreases the error to a minimum of
+ * 1.1e-16 but the program can't achieve the requested D.
+ */
 int main(int argc, char** argv) {
   double exact =
     0.57721566490153286060651209008240243104215933593992;
 
-  // If no argument is given is input, an error signal is displayed
+  /* if no argument is given, show the usage */
-  // and the program quits
-
   if (argc != 2) {
     fprintf(stderr, "usage: %s D\n", argv[0]);
     fprintf(stderr, "Computes γ up to D decimal places.\n");
     return EXIT_FAILURE;
   }
 
-  int N = floor(2.0 + 1.0/4 * log(10) * (double)atoi(argv[1]));
+  // requested decimal digits
+  int D = atoi(argv[1]);
+  // Brent-McMillan number N
+  int N = gsl_sf_lambert_W0(5*M_PI/9*pow(10, 2*D + 1))/16 + 1;
+
   double A = a_series(N);
   double B = b_series(N);
   double C = c_series(N);