How to resolve the algorithm Euler's constant 0.5772... step by step in the Processing programming language

Published on 12 May 2024 09:40 PM
#Processing

How to resolve the algorithm Euler's constant 0.5772... step by step in the Processing programming language

Table of Contents

Problem Statement

Compute the Euler constant 0.5772... Discovered by Leonhard Euler around 1730, it is the most ubiquitous mathematical constant after pi and e, but appears more arcane than these. Denoted gamma (γ), it measures the amount by which the partial sums of the harmonic series (the simplest diverging series) differ from the logarithmic function (its approximating integral): lim n → ∞ (1 + 1/2 + 1/3 + … + 1/n − log(n)). The definition of γ converges too slowly to be numerically useful, but in 1735 Euler himself applied his recently discovered summation formula to compute ‘the notable number’ accurate to 15 places. For a single-precision implementation this is still the most economic algorithm. In 1961, the young Donald Knuth used Euler's method to evaluate γ to 1271 places. Knuth found that the computation of the Bernoulli numbers required in the Euler-Maclaurin formula was the most time-consuming part of the procedure. The next year Dura Sweeney computed 3566 places, using a formula based on the expansion of the exponential integral which didn't need Bernoulli numbers. It's a bit-hungry method though: 2d digits of working precision obtain d correct places only. This was remedied in 1988 by David Bailey; meanwhile Richard Brent and Ed McMillan had published an even more efficient algorithm based on Bessel function identities and found 30100 places in 20 hours time. Nowadays the old records have far been exceeded: over 6·1011 decimal places are already known. These massive computations suggest that γ is neither rational nor algebraic, but this is yet to be proven.

[1] Gourdon and Sebah, The Euler constant γ. (for all formulas) [2] Euler's original journal article translated from the latin (p. 9)

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Euler's constant 0.5772... step by step in the Processing programming language

Source code in the processing programming language

/*********************************************
 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");
}

You may also check:How to resolve the algorithm Matrix multiplication step by step in the Octave programming language
You may also check:How to resolve the algorithm Letter frequency step by step in the Tcl programming language
You may also check:How to resolve the algorithm Ethiopian multiplication step by step in the CLU programming language
You may also check:How to resolve the algorithm Sort using a custom comparator step by step in the FutureBasic programming language
You may also check:How to resolve the algorithm Hamming numbers step by step in the ATS programming language