#include <stdio.h>
#include <math.h>
#include <stdlib.h>

// The Euler-Mascheroni constant is computed through the formula:
//
//    γ = A(N)/B(N) - ln(N)
//
// with:
//
//    A(N) = Σ_(k = 1)^(k_max) (N^k/k!) * H(k)
//    B(N) = Σ_(k = 0)^(k_max) (N^k/k!)
//    H(k) = Σ_(j = 1)^(k) (1/k)
//
// where N is computed from D as written below and k_max is the value
// at which there is no difference between two consecutive terms of
// the sum because of double precision.
//
// source: http://www.numberworld.org/y-cruncher/internals/formulas.html

// 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;
}

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);
}

// Takes in input the number D of desired correct decimals
// Best result obtained with D = 15, N = 10

int main(int argc, char** argv) {
  double exact =
    0.57721566490153286060651209008240243104215933593992;

  // If no argument is given is input, an error signal is displayed
  // 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]));
  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;
}