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;
|
|||
|
|
}
|