512 lines
13 KiB
C
512 lines
13 KiB
C
#include <stdio.h>
|
||
#include <stdlib.h>
|
||
#include <stdarg.h>
|
||
#include <math.h>
|
||
#include <gmp.h>
|
||
|
||
/* The Euler-Mascheroni constant is here computed to
|
||
* arbitrary precision using GMP rational numbers
|
||
* with the formula:
|
||
*
|
||
* γ = A(N)/B(N) - C(N)/B(N)² - log(N)
|
||
*
|
||
* with:
|
||
*
|
||
* 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
|
||
*
|
||
* where N is computed from D as written below and k_max is
|
||
* currently 5*N, chosen in a completely arbitrary way.
|
||
*
|
||
* source: http://www.numberworld.org/y-cruncher/internals/formulas.html
|
||
*
|
||
*/
|
||
|
||
/* `mpq_dec_str(z, n)`
|
||
* returns the decimal representation of the
|
||
* rational z as a string, up to n digits.
|
||
*/
|
||
char* mpq_dec_str(mpq_t z, size_t mantissa_size) {
|
||
mpz_t num, den, quot, rem;
|
||
mpz_inits(num, den, quot, rem, NULL);
|
||
|
||
// calculate integral part
|
||
mpq_get_num(num, z);
|
||
mpq_get_den(den, z);
|
||
mpz_tdiv_qr(quot, rem, num, den);
|
||
|
||
// store the integral part
|
||
size_t integral_size = mpz_sizeinbase(quot, 10);
|
||
char* repr = calloc(integral_size + 1 + mantissa_size + 1, sizeof(char));
|
||
mpz_get_str(repr, 10, quot);
|
||
|
||
// add decimal point
|
||
repr[integral_size] = '.';
|
||
|
||
// calculate the mantissa
|
||
mpz_t cur;
|
||
mpz_init(cur);
|
||
size_t digit;
|
||
for (size_t i = 0; i < mantissa_size; i++) {
|
||
// calculate i-th digit as:
|
||
// digit = (remainder * 10) / denominator
|
||
// new_ramainder = (remainder * 10) % denominator
|
||
mpz_set(cur, rem);
|
||
mpz_mul_ui(cur, rem, 10);
|
||
mpz_tdiv_qr(quot, rem, cur, den);
|
||
|
||
// store digit in mantissa
|
||
digit = mpz_get_ui(quot);
|
||
sprintf(repr + integral_size + 1 + i, "%ld", digit);
|
||
}
|
||
|
||
// free memory
|
||
mpz_clears(num, den, quot, rem, cur, NULL);
|
||
|
||
return repr;
|
||
}
|
||
|
||
|
||
void mpq_fprintf(FILE *fp, size_t digits, char* fmt, ...) {
|
||
va_list args;
|
||
va_start(args, fmt);
|
||
|
||
char* str;
|
||
char fmtc[2] = "%x";
|
||
|
||
for (char* p = fmt; *p; p++) {
|
||
if (*p != '%') {
|
||
fprintf(fp, "%c", *p);
|
||
continue;
|
||
}
|
||
switch (*++p) {
|
||
case 's':
|
||
str = va_arg(args, char*);
|
||
fputs(str, fp);
|
||
break;
|
||
case 'l':
|
||
if (*++p == 'd')
|
||
fprintf(fp, "%ld", va_arg(args, size_t));
|
||
break;
|
||
case 'q':
|
||
str = mpq_dec_str(*va_arg(args, mpq_t*), digits);
|
||
fputs(str, fp);
|
||
free(str);
|
||
break;
|
||
default:
|
||
fmtc[1] = *p;
|
||
gmp_fprintf(fp, fmtc, *va_arg(args, mpq_t*));
|
||
}
|
||
}
|
||
va_end(args);
|
||
}
|
||
|
||
|
||
/* `mpq_log2(z, D)`
|
||
* calculate log(2) to D decimal places
|
||
* using the series
|
||
*
|
||
* log(2) = Σ_{k=1} 1/(k 2^k)
|
||
*
|
||
* The error of the n-th partial sum is
|
||
* estimated as
|
||
*
|
||
* Δn = 1/(n(n+2) 2^(n-1)).
|
||
*
|
||
* So to fix the first D decimal places
|
||
* we evaluate up to n terms such that
|
||
*
|
||
* Δn < 1/10^D
|
||
*/
|
||
void mpq_log2(mpq_t rop, size_t digits) {
|
||
mpq_t sum, term;
|
||
mpq_inits(sum, term, NULL);
|
||
|
||
fprintf(stderr, "[mpq_log2] decimal digits: %ld\n", digits);
|
||
|
||
/* calculate n, the number of required terms */
|
||
size_t n;
|
||
mpz_t temp, prec, width; mpz_inits(temp, prec, width, NULL);
|
||
// prec = 10^digits
|
||
mpz_ui_pow_ui(prec, 10, digits + 2);
|
||
|
||
for (n = 1;; n++) {
|
||
// width = n*(n + 2)*2^(n - 1)
|
||
mpz_set_ui(width, n);
|
||
mpz_set_ui(temp, n);
|
||
mpz_add_ui(temp, temp, 2);
|
||
mpz_mul(width, width, temp);
|
||
mpz_sub_ui(temp, temp, 3);
|
||
mpz_mul_2exp(width, width, mpz_get_ui(temp));
|
||
|
||
if (mpz_cmp(width, prec) > 0) break;
|
||
}
|
||
fprintf(stderr, "[mpq_log2] series terms: %ld\n", n);
|
||
|
||
/* calculate series */
|
||
for (size_t i=1; i<n+1; i++) {
|
||
mpq_set_ui(term, 1, 1);
|
||
mpz_ui_pow_ui(mpq_denref(term), 2, i);
|
||
mpz_mul_ui(mpq_denref(term), mpq_denref(term), i);
|
||
mpq_add(sum, sum, term);
|
||
}
|
||
mpq_fprintf(stderr, digits > 50 ? 50 : digits,
|
||
"[mpq_log2] log(2)=%q...\n", &sum);
|
||
|
||
mpq_set(rop, sum);
|
||
|
||
// free memory
|
||
mpq_clears(sum, term, NULL);
|
||
mpz_clears(temp, prec, width, NULL);
|
||
}
|
||
|
||
|
||
/* `mpq_log_ui(z, x, D)`
|
||
* calculates the logarithm of x up
|
||
* to D decimal places using the formula
|
||
*
|
||
* log(x) = log(a(x) * 2^(n-1))
|
||
* = log(a) + (n-1) log(2)
|
||
*
|
||
* where 1<a<2 and n is the number of binary
|
||
* digits of x.
|
||
*
|
||
* log(a) is the computed from
|
||
*
|
||
* log((1+y)/(1-y)) = 2Σ_{k=0} y^(2k+1)/(2k+1)
|
||
*
|
||
* where y = (a-1)/(a+1).
|
||
*
|
||
* The error of the n-th partial sum is
|
||
* estimated as:
|
||
*
|
||
* Δ(n, y) = 4y^(2n + 3) /
|
||
* (4n^2y^4 - 8n^2y^2 + 4n^2 + 4ny^4 - 12ny^2 + 8n + y^4 - 4y^2 + 3)
|
||
*
|
||
* So to fix the first D decimal places
|
||
* we evaluate up to n terms such that
|
||
*
|
||
* Δ(n, y) < 1/10^D
|
||
*/
|
||
void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
|
||
/* calculate `a` such that x = a⋅2^(m-1) */
|
||
mpq_t a; mpq_init(a);
|
||
mpz_t xz, power; mpz_inits(xz, power, NULL);
|
||
mpz_set_ui(xz, x);
|
||
|
||
// m = digits in base 2 of x
|
||
size_t m = mpz_sizeinbase(xz, 2);
|
||
mpz_ui_pow_ui(power, 2, m - 1);
|
||
mpq_set_num(a, xz);
|
||
mpq_set_den(a, power);
|
||
|
||
fprintf(stderr, "[mpq_log_ui] binary digits: %ld\n", m);
|
||
fprintf(stderr, "[mpq_log_ui] a=%g\n", mpq_get_d(a));
|
||
|
||
/* calculate y = (a-1)/(a+1) */
|
||
mpq_t y, temp; mpq_inits(y, temp, NULL);
|
||
|
||
mpq_set(temp, a);
|
||
// decrement: a/b -= b/b
|
||
mpz_submul_ui(mpq_numref(temp),
|
||
mpq_denref(temp), 1);
|
||
mpq_set(y, temp);
|
||
|
||
mpq_set(temp, a);
|
||
// increment: a/b += b/b
|
||
mpz_addmul_ui(mpq_numref(temp),
|
||
mpq_denref(temp), 1);
|
||
mpq_div(y, y, temp);
|
||
|
||
fprintf(stderr, "[mpq_log_ui] y=%g\n", mpq_get_d(y));
|
||
|
||
|
||
/* calculate number of required terms:
|
||
* find smallest n such that Δ(n, y) < 1/10^digits.
|
||
*/
|
||
mpq_t y4, y2, y2n, t1, t2;
|
||
mpq_inits(y4, y2, y2n, t1, t2, NULL);
|
||
|
||
// y4 = y^4
|
||
mpz_pow_ui(mpq_numref(y4), mpq_numref(y), 4);
|
||
mpz_pow_ui(mpq_denref(y4), mpq_denref(y), 4);
|
||
|
||
// y2 = y^2
|
||
mpz_pow_ui(mpq_numref(y2), mpq_numref(y), 2);
|
||
mpz_pow_ui(mpq_denref(y2), mpq_denref(y), 2);
|
||
|
||
size_t n;
|
||
size_t p1, p2, p3;
|
||
|
||
for(n = 0;; n++) {
|
||
// polynomials in n
|
||
p1 = 4*pow(n, 2) + 4*n + 1;
|
||
p2 = 8*pow(n, 2) + 12*n - 4;
|
||
p3 = 4*pow(n, 2) + 3;
|
||
|
||
// t1 = y4 * p1
|
||
mpz_mul_ui(mpq_numref(t1), mpq_numref(y4), p1);
|
||
mpq_canonicalize(t1);
|
||
// t2 = y2 * p2
|
||
mpz_mul_ui(mpq_numref(t2), mpq_numref(y2), p2);
|
||
mpq_canonicalize(t2);
|
||
|
||
// t1 -= t2
|
||
mpq_sub(t1, t1, t2);
|
||
// t1 += p3
|
||
mpz_addmul_ui(mpq_numref(t1), mpq_denref(t1), p3);
|
||
|
||
// y2n = 4*y^(2n + 3)
|
||
mpq_set(y2n, y);
|
||
mpz_pow_ui(mpq_numref(y2n), mpq_numref(y2n), 2*n+3);
|
||
mpz_pow_ui(mpq_denref(y2n), mpq_denref(y2n), 2*n+3);
|
||
mpz_mul_ui(mpq_numref(y2n), mpq_numref(y2n), 4);
|
||
mpq_canonicalize(y2n);
|
||
|
||
// t2 = 10^digits * y2n
|
||
mpq_set_ui(t2, 1, 1);
|
||
mpz_ui_pow_ui(mpq_numref(t2), 10, digits + 1);
|
||
mpq_mul(t2, t2, y2n);
|
||
|
||
if (mpq_cmp(t1, t2) > 0) break;
|
||
}
|
||
|
||
fprintf(stderr, "[mpq_log_ui] series terms: %ld\n", n);
|
||
fprintf(stderr, "[mpq_log_ui] decimal digits: %ld\n", digits);
|
||
|
||
|
||
/* calculate the series of
|
||
* log(A) = log((1+y)/(1-y)) = 2 Σ_{k=0} y^(2k+1)/(2k+1)
|
||
*/
|
||
mpq_t sum, term; mpq_inits(sum, term, NULL);
|
||
|
||
for (size_t i = 0; i <= n; i++) {
|
||
// numerator: y^(2i+1)
|
||
mpq_set(temp, y);
|
||
mpz_pow_ui(mpq_numref(temp), mpq_numref(temp), 2*i + 1);
|
||
mpz_pow_ui(mpq_denref(temp), mpq_denref(temp), 2*i + 1);
|
||
mpq_canonicalize(temp);
|
||
mpq_set(term, temp);
|
||
|
||
// denominator: 2i+1
|
||
mpq_set_ui(temp, 2*i + 1, 1);
|
||
mpq_div(term, term, temp);
|
||
|
||
// sum the term
|
||
mpq_add(sum, sum, term);
|
||
}
|
||
|
||
// multiply by 2
|
||
mpz_mul_ui(mpq_numref(sum), mpq_numref(sum), 2);
|
||
mpq_canonicalize(sum);
|
||
|
||
// calculate log(2)*(m-1)
|
||
mpq_t log2; mpq_init(log2);
|
||
mpq_log2(log2, digits);
|
||
|
||
mpz_mul_ui(mpq_numref(log2), mpq_numref(log2), m-1);
|
||
mpq_canonicalize(log2);
|
||
|
||
// series += log(2)*(m-1)
|
||
mpq_add(sum, sum, log2);
|
||
mpq_set(rop, sum);
|
||
|
||
mpq_fprintf(stderr, digits > 50 ? 50 : digits,
|
||
"[mpq_log_ui] log(%ld)=%q...\n", x, sum);
|
||
|
||
// free memory
|
||
mpq_clears(a, y, temp, sum, term, log2,
|
||
y4, y2, y2n, t1, t2, NULL);
|
||
mpz_clears(xz, power, NULL);
|
||
}
|
||
|
||
|
||
/* A series */
|
||
void a_series(mpq_t rop, size_t N) {
|
||
mpq_t sum, term, factq, harm, recip;
|
||
mpz_t fact, power, stop;
|
||
|
||
mpq_inits(sum, term, factq, harm, recip, NULL);
|
||
mpz_inits(fact, power, stop, NULL);
|
||
mpz_set_ui(stop, N);
|
||
|
||
for (size_t k = 0; k < 5*N; k++) {
|
||
mpz_fac_ui(fact, k);
|
||
mpz_pow_ui(power, stop, k);
|
||
|
||
mpq_set_z(term, power);
|
||
mpq_set_z(factq, fact);
|
||
mpq_div(term, term, factq);
|
||
mpq_mul(term, term, term);
|
||
|
||
if (k > 0)
|
||
mpq_set_ui(recip, 1, k);
|
||
mpq_add(harm, harm, recip);
|
||
mpq_mul(term, term, harm);
|
||
mpq_add(sum, sum, term);
|
||
}
|
||
|
||
mpq_set(rop, sum);
|
||
|
||
// free memory
|
||
mpq_clears(sum, term, factq, harm, recip, NULL);
|
||
mpz_clears(fact, power, stop, NULL);
|
||
}
|
||
|
||
|
||
/* B series */
|
||
void b_series(mpq_t rop, size_t N){
|
||
mpq_t sum, term, factq;
|
||
mpz_t fact, power, stop;
|
||
|
||
mpq_inits(sum, term, factq, NULL);
|
||
mpz_inits(fact, power, stop, NULL);
|
||
mpz_set_ui(stop, N);
|
||
|
||
for (size_t k = 0; k < 5*N; k++) {
|
||
mpz_fac_ui(fact, k);
|
||
mpz_pow_ui(power, stop, k);
|
||
|
||
mpq_set_z(term, power);
|
||
mpq_set_z(factq, fact);
|
||
mpq_div(term, term, factq);
|
||
|
||
mpq_mul(term, term, term);
|
||
mpq_add(sum, sum, term);
|
||
}
|
||
|
||
mpq_set(rop, sum);
|
||
|
||
// free memory
|
||
mpq_clears(sum, term, factq, NULL);
|
||
mpz_clears(fact, power, stop, NULL);
|
||
}
|
||
|
||
|
||
/* C series */
|
||
void c_series(mpq_t rop, size_t N) {
|
||
mpq_t sum, term, stop2, fact2q;
|
||
mpz_t fact1, fact2, power, stop1;
|
||
|
||
mpq_inits(sum, term, stop2, fact2q, NULL);
|
||
mpz_inits(fact1, fact2, power, stop1, NULL);
|
||
mpq_set_ui(stop2, 4*N, 1);
|
||
|
||
for (size_t k = 0; k <= 2*N; k++) {
|
||
mpz_fac_ui(fact1, 2*k);
|
||
mpz_pow_ui(fact1, fact1, 3);
|
||
|
||
mpz_fac_ui(fact2, k);
|
||
mpz_pow_ui(fact2, fact2, 4);
|
||
|
||
mpz_set_ui(stop1, 16*N);
|
||
mpz_pow_ui(stop1, stop1, 2*k);
|
||
mpz_mul(fact2, fact2, stop1);
|
||
|
||
mpq_set_z(term, fact1);
|
||
mpq_set_z(fact2q, fact2);
|
||
|
||
mpq_div(term, term, fact2q);
|
||
mpq_add(sum, sum, term);
|
||
}
|
||
mpq_div(sum, sum, stop2);
|
||
|
||
mpq_set(rop, sum);
|
||
|
||
// free memory
|
||
mpq_clears(sum, term, stop2, fact2q, NULL);
|
||
mpz_clears(fact1, fact2, power, stop1, NULL);
|
||
}
|
||
|
||
|
||
int main(int argc, char** argv) {
|
||
/* 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;
|
||
}
|
||
|
||
mpq_t exact, indicator;
|
||
mpq_inits(exact, indicator, NULL);
|
||
|
||
// value of γ up to 500 digits
|
||
mpq_set_str(
|
||
exact,
|
||
"47342467929426395252720948664321583815681737620205"
|
||
"78046578043524355342995304421105234825234769529727"
|
||
"60710541774258570118818437511973945213031180449905"
|
||
"77067110892216494927227098092412641670668896104172"
|
||
"59380467495735805599918127376547804186055507379270"
|
||
"09659478275476110597673948659342592173551684412254"
|
||
"47940426607883070678658320829056222506449877887289"
|
||
"02099638646695565284675455912794915138438516540912"
|
||
"22061900908485016328626399040675802962645141413722"
|
||
"567840383761005003563012969560659130635723002086605/"
|
||
"82018681765164048732047980836753451023877954010252"
|
||
"60062364748361667340168652059998708337602423525120"
|
||
"45225158774173869894826877890589130978987229877889"
|
||
"33367849273189687823618289122425446493605087108634"
|
||
"04387981302669131224273324182166778131513056804533"
|
||
"58955006355665628938266331979307689540884269372365"
|
||
"76288367811322713649805442241450184023209087215891"
|
||
"55369788474437679223152173114447113970483314961392"
|
||
"48250188991402851129033493732164230227458717486395"
|
||
"514436574417275149404197774547389507462779807727616",
|
||
10);
|
||
|
||
// number with decimal expansion equal to
|
||
// 0123456789 periodic: useful to count digits
|
||
mpq_set_str(
|
||
indicator,
|
||
"137174210/1111111111",
|
||
10);
|
||
|
||
// number of decimal places
|
||
size_t D = atol(argv[1]);
|
||
// number of terms
|
||
size_t N = floor(1.0/4 * log(10) * (double)D);
|
||
|
||
/* calculate series */
|
||
mpq_t A, B, C;
|
||
mpq_inits(A, B, C, NULL);
|
||
a_series(A, N);
|
||
b_series(B, N);
|
||
c_series(C, N);
|
||
|
||
/* calculate γ... */
|
||
mpq_t gamma, corr, logn;
|
||
mpq_inits(gamma, corr, logn, NULL);
|
||
|
||
mpq_log_ui(logn, N, D);
|
||
mpq_div(gamma, A, B);
|
||
mpq_div(corr, C, B);
|
||
mpq_div(corr, corr, B);
|
||
mpq_sub(gamma, gamma, corr);
|
||
mpq_sub(gamma, gamma, logn);
|
||
|
||
/* calculate difference */
|
||
mpq_t diff;
|
||
mpq_init(diff);
|
||
mpq_sub(diff, gamma, exact);
|
||
|
||
/* print results */
|
||
size_t d = D+1 > 50 ? 50 : D+1;
|
||
|
||
fprintf(stderr, "[main] N: %ld\n", N);
|
||
mpq_fprintf(stderr, d, "[main] approx:\t%q...\n", gamma);
|
||
mpq_fprintf(stderr, d, "[main] exact:\t%q...\n", exact);
|
||
mpq_fprintf(stderr, d, "[main] diff:\t%q...\n", diff);
|
||
mpq_fprintf(stderr, d, "[main] digits:\t%q\n", indicator);
|
||
|
||
mpq_fprintf(stdout, D+1, "%q\n", gamma);
|
||
|
||
// free memory
|
||
mpq_clears(A, B, C, gamma, corr, logn,
|
||
exact, indicator, diff, NULL);
|
||
|
||
return EXIT_SUCCESS;
|
||
}
|