61 lines
1.4 KiB
C
61 lines
1.4 KiB
C
#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;
|
||
}
|