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

/* Compute the Euler_Mascheroni constant through the formula:
 *
 * γ = lim_{M → +inf} sum_{k  = 1}^{M} bin{M, k} (-1)^k / k * ln(Γ(k+1))
 *
 * Best result is obtained for M = 41.
 * Values from 1 to 100 were checked.
 *
 */

double exact = 0.577215664901532860;

int main() {
  double best = 0, num = 0;
  double sum = 0, prev;

  for (double M = 1; M < 60; M++) {
    prev = sum;
    sum = 0;
    for (double k = 1; k <= M; k++) {
      sum += (tgamma(M + 1))/(tgamma(k + 1)*tgamma(M - k + 1)) * pow(-1, k)/k
        * log(tgamma (k + 1));
    }
    printf("M: \t%d\t", (int)M);
    printf("diff:\t%.15f\n", fabs(sum - exact));
    if ((fabs(sum - exact) < fabs(prev - exact)) &&
        (fabs(sum - exact) < fabs(best - exact))){
      best = sum;
      num = M;
    }
  }

  printf("M: %d\n", (int)num);
  printf("approx:\t%.15f\n", best);
  printf("true:\t%.15f\n", exact);
  printf("diff:\t%.15f\n", fabs(best - exact));
  printf("\t  123456789_12345\n");

  return EXIT_SUCCESS;
}