#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_sf_lambert.h>

/* The Euler-Mascheroni constant is computed to D decimal
 * places using the refined Brent-McMillan formula:
 *
 *    γ = A(N)/B(N) - C(N)/B(N)² - log(N)
 *
 * 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
 *
 * The error of the estimation is given by
 *
 *   Δ(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
double harmonic_sum(double n) {
  double sum = 0;
  for (double k = 1; k < n+1; k++) {
    sum += 1/k;
  }
  return sum;
}

// A series 
double a_series(int N) {
  double sum = 0;
  double prev = -1;
  for (double k = 1; sum != prev; k++) {
    prev = sum;
    sum += pow(((pow(N, k))/(tgamma(k+1))), 2) * harmonic_sum(k);
  }
  return sum;
}

// B series
double b_series(int N){
  double sum = 0;
  double prev = -1;
  for (double k = 0; sum != prev; k++) {
    prev = sum;
    sum += pow(((pow(N, k))/(tgamma(k+1))), 2);
  }
  return sum;
}

// C series
double c_series(int N) {
  double sum = 0;
  for (double k = 0; k < N; k++) {
    sum += pow(tgamma(2*k + 1), 3)/(pow(tgamma(k + 1), 4) * pow(16.0*N, (int)2*k));
  }
  return sum/(4.0*N);
}

/* The Best result is obtained with D=15, which accurately
 * 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, show the usage */
  if (argc != 2) {
    fprintf(stderr, "usage: %s D\n", argv[0]);
    fprintf(stderr, "Computes γ up to D decimal places.\n");
    return EXIT_FAILURE;
  }

  // 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);
  double gamma = A/B - C/(B*B) - log(N);

  printf("N: %d\n", N);
  printf("approx:\t%.30f\n", gamma);
  printf("true:\t%.30f\n", exact);
  printf("diff:\t%.30f\n", fabs(gamma - exact));
  printf("\t  123456789 123456789 123456789\n");

  return EXIT_SUCCESS;
}