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