112 lines
2.9 KiB
C
112 lines
2.9 KiB
C
#include <stdio.h>
|
||
#include <stdlib.h>
|
||
#include <math.h>
|
||
#include <gsl/gsl_sf_lambert.h>
|
||
|
||
/* The Euler-Mascheroni constant is computed to D decimal
|
||
* places using the refined Brent-McMillan formula:
|
||
*
|
||
* γ = A(N)/B(N) - C(N)/B(N)² - log(N)
|
||
*
|
||
* where:
|
||
*
|
||
* A(N) = Σ_(k = 1)^(k_max) (N^k/k!)² * H(k)
|
||
* B(N) = Σ_(k = 0)^(k_max) (N^k/k!)²
|
||
* C(N) = 1/4N Σ_(k = 0)^(2N) ((2k)!)^3/((k!)^4*(16N))^2k
|
||
* H(k) = Σ_(j = 1)^(k) (1/k), the k-th harmonic number
|
||
*
|
||
* The error of the estimation is given by
|
||
*
|
||
* Δ(N) = 5/12 √(2π) exp(-8N)/√x
|
||
*
|
||
* so, given a number D of decimal digits to compute we ask
|
||
* that the error be smaller than 10^-D. The smallest integer
|
||
* to solve the inequality can be proven to be
|
||
*
|
||
* N = floor(W(5π/9 10^(2D + 1))/16) + 1
|
||
*
|
||
* where W is the principal value of the Lambert W function.
|
||
*
|
||
* The series are truncated at k_max, which is when the difference
|
||
* between two consecutive terms of the sum is zero, at double precision.
|
||
*
|
||
* source: Precise error estimate of the Brent-McMillan algorithm for
|
||
* the computation of Euler's constant, Jean-Pierre Demailly.
|
||
*/
|
||
|
||
// 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;
|
||
}
|
||
|
||
// C series
|
||
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);
|
||
}
|
||
|
||
/* The Best result is obtained with D=15, which accurately
|
||
* computes 15 digits and gives an error of 3.3e-16.
|
||
*
|
||
* Increasing to D=21 decreases the error to a minimum of
|
||
* 1.1e-16 but the program can't achieve the requested D.
|
||
*/
|
||
int main(int argc, char** argv) {
|
||
double exact =
|
||
0.57721566490153286060651209008240243104215933593992;
|
||
|
||
/* if no argument is given, show the usage */
|
||
if (argc != 2) {
|
||
fprintf(stderr, "usage: %s D\n", argv[0]);
|
||
fprintf(stderr, "Computes γ up to D decimal places.\n");
|
||
return EXIT_FAILURE;
|
||
}
|
||
|
||
// requested decimal digits
|
||
int D = atoi(argv[1]);
|
||
// Brent-McMillan number N
|
||
int N = gsl_sf_lambert_W0(5*M_PI/9*pow(10, 2*D + 1))/16 + 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;
|
||
}
|