90 lines
2.1 KiB
C
90 lines
2.1 KiB
C
#include <stdio.h>
|
||
#include <math.h>
|
||
#include <stdlib.h>
|
||
|
||
// The Euler-Mascheroni constant is computed through the formula:
|
||
//
|
||
// γ = A(N)/B(N) - ln(N)
|
||
//
|
||
// with:
|
||
//
|
||
// A(N) = Σ_(k = 1)^(k_max) (N^k/k!) * H(k)
|
||
// B(N) = Σ_(k = 0)^(k_max) (N^k/k!)
|
||
// H(k) = Σ_(j = 1)^(k) (1/k)
|
||
//
|
||
// where N is computed from D as written below and k_max is the value
|
||
// at which there is no difference between two consecutive terms of
|
||
// the sum because of double precision.
|
||
//
|
||
// source: http://www.numberworld.org/y-cruncher/internals/formulas.html
|
||
|
||
// 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;
|
||
}
|
||
|
||
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);
|
||
}
|
||
|
||
// Takes in input the number D of desired correct decimals
|
||
// Best result obtained with D = 15, N = 10
|
||
|
||
int main(int argc, char** argv) {
|
||
double exact =
|
||
0.57721566490153286060651209008240243104215933593992;
|
||
|
||
// If no argument is given is input, an error signal is displayed
|
||
// and the program quits
|
||
|
||
if (argc != 2) {
|
||
fprintf(stderr, "usage: %s D\n", argv[0]);
|
||
fprintf(stderr, "Computes γ up to D decimal places.\n");
|
||
return EXIT_FAILURE;
|
||
}
|
||
|
||
int N = floor(2.0 + 1.0/4 * log(10) * (double)atoi(argv[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;
|
||
}
|