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

/* Compute the Euler-Mascheroni constant through the formula
 * of the reciprocal Γ function:
 *
 *   1/Γ(z) = z*e^{γz} pro_{k = 1}^{+ inf} ((1 + z/k) e^{- z/k})
 *
 * Thus:
 *
 *   γ = - 1/z * ln(z*Γ(z)*pro{k = 1}^{+ inf} (1 + z/k) e^{- z/k})
 *
 * Product stops when there is no difference between two consecutive
 * terms of the productory due to limited precision.
 *
 * Range of z given in input as min, max and step.
 *
 */

double exact = 0.577215664901532860;

int main(int argc, char** argv) {
  double min  = 8.95;
  double max  = 9.04;
  double step = 0.01;

  double prev_gamma, gamma = 0;
  double Z;
  double best = 0;

  for (double z = min; z <= max; z += step) {
    prev_gamma = gamma;

    double pro = 1;
    double prev = -1;
    for (double k = 1; pro != prev; k++) {
      prev = pro;
      pro *= (1 + z/k) * exp(- z/k);
      if (z == 9 && pro == prev) printf("k:\t%.0f\n", k);
    }
    double gamma = - 1/z * log(z*tgamma(z)*pro);
    printf("z:\t%.2f\t", z);
    printf("diff:\t%.20f\n", fabs(gamma - exact));

    if ((fabs(gamma - exact) < fabs(prev_gamma - exact)) &&
        (fabs(gamma - exact) < fabs(best - exact))){
      best = gamma;
      Z = z;
    }
  }

  printf("z:\t%.2f\n", Z);
  printf("approx:\t%.20f\n", best);
  printf("true:\t%.20f\n", exact);
  printf("diff:\t%.20f\n", best - exact);
  printf("\t   123456789 123456789 123456789\n");

  return EXIT_SUCCESS;
}