/*********************************************
 Subject: Comparing five methods for
 computing Euler's constant 0.5772...
 // https://rosettacode.org/wiki/Euler%27s_constant_0.5772...
 --------------------------------------------*/
double a, b, h, n2, r, u, v;
float floatA, floatB, floatN2;
int k, k2, m, n;
double eps = 1e-6;

void setup() {
  size(100, 100);
  noLoop();
}

void draw() {
  println("From the definition, err. 3e-10\n");

  n = 400;

  h = 1;

  for (int k = 2; k <= n; k++) {
    h += 1.0 / k;
  }
  //faster convergence: Negoi, 1997
  a = log(n +.5 + 1.0 / (24*n));

  println("Hn    ", h);
  println("gamma ", h - a);
  println("k = ", n);
  println("");


  println("Sweeney, 1963, err. idem");
  n = 21;

  double s[] = {0, n};
  r = n;
  k = 1;
  while (r > eps) {
    k ++;
    r *= (double) n / k;
    s[k & 1] = s[k & 1] + r / k;
  }

  // println("gamma %.16f\nk = %d\n\n", s[1] - s[0] - log(n), k);

  println("Hn    ", h);
  println("gamma ", s[1] - s[0] - log(n));
  println("k = ", k);
  println("");

  println("Bailey, 1988");
  n = 5;
  floatA = 1;
  h = 1;
  floatN2 = pow(2, n);
  r = 1;
  k = 1;
  while (abs(floatB - floatA) > eps) {
    k += 1;
    r *= floatN2 / k;
    h += 1.0 / k;
    floatB = floatA;
    floatA += r * h;
  }
  floatA *= floatN2 / exp(floatN2);

  println("gamma ", floatA - n * log(2));
  println("k = ", k);
  println("");

  println("Brent-McMillan, 1980");

  n = 13;

  floatA = -log(n);
  floatB = 1;
  u = a;
  v = b;
  n2 = n * n;
  k2 = 0;
  k = 0;


  while (abs(floatA) > eps) {
    k2 += 2*k + 1;
    k += 1;
    floatA *= n2 / k;
    floatB *= n2 / k2;
    floatA = (floatA + floatB) / k;
    u += floatA;
    v += floatB;
  }
  println("gamma  ", u / v);
  println("k      ", k);



  println("How Euler did it in 1735\n");
  //Bernoulli numbers with even indices

  double[] B2 = new double[11];

  B2[1] = 1.0;
  B2[2] = 1.0/6;
  B2[3] = -1.0/30;
  B2[4] = 1.0/42;
  B2[5] = -1.0/30;
  B2[6] = 5.0/66;
  B2[7] = -691.0/2730;
  B2[8] = 7.0/6;
  B2[9] = -3617.0/510;
  B2[10]= 43867.0/798;

  m = 7;
  n = 10;

  //n-th harmonic number
  h = 1;
  for (k = 2; k <= n; k++) {
    h += 1.0 / k;
  }
  println("Hn    ", h);

  h -= log(n);
  println("  -ln ", h);

  //expansion C = -digamma(1)
  a = -1.0 / (2*n);
  n2 = n * n;
  r = 1;
  for (k = 1; k <= m; k++) {
    r *= n2;
    a += B2[k] / (2*k * r);
  }

  println("");
  println("err  ", a);
  println("gamma ", h + a );
  println("k = ", n + m);
  println("");
  println("C  =  0.57721566490153286...\n");
}