From c8f34f1822467cae60201560e2de422624612e71 Mon Sep 17 00:00:00 2001
From: rnhmjoj <rnhmjoj@inventati.org>
Date: Sun, 24 May 2020 12:49:47 +0000
Subject: [PATCH] ex-2: implement correct series estimation in fancier.c

---
 ex-2/fancier.c | 63 +++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 57 insertions(+), 6 deletions(-)

diff --git a/ex-2/fancier.c b/ex-2/fancier.c
index c2bc256..a0d3074 100644
--- a/ex-2/fancier.c
+++ b/ex-2/fancier.c
@@ -146,20 +146,42 @@ void mpq_log2(mpq_t rop, size_t digits) {
   }
   fprintf(stderr, "[mpq_log2] series terms: %ld\n", n);
 
-  /* calculate series */
-  for (size_t i=1; i<n+1; i++) {
+  /* Compute the nth partial sum */
+  for (size_t i=1; i<=n; 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);
 
+  /* Compute upper/lower bounds */
+  mpq_t upp, low;
+  mpq_inits(upp, low, NULL);
+
+  // upp = sum + term
+  mpq_add(upp, sum, term);
+
+  // low = sum + term * n/(n+2)
+  mpz_mul_ui(mpq_numref(low), mpq_numref(term), n);
+  mpz_mul_ui(mpq_denref(low), mpq_denref(term), n+2);
+  mpq_canonicalize(low);
+  mpq_add(low, sum, low);
+
+  /* Average the two bounds */
+  mpq_add(sum, upp, low);
+  mpz_mul_ui(mpq_denref(sum), mpq_denref(sum), 2);
+  mpq_canonicalize(sum);
+
+  size_t d = digits > 50 ? 50 : digits;
+  mpq_fprintf(stderr, d+2, "[mpq_log2] upper bound: %q...\n", &upp);
+  mpq_fprintf(stderr, d+2, "[mpq_log2] lower bound: %q...\n", &low);
+  mpq_fprintf(stderr, d,   "[mpq_log2] log(2)=%q...\n", &sum);
+
+  // return sum
   mpq_set(rop, sum);
 
   // free memory
-  mpq_clears(sum, term, NULL);
+  mpq_clears(sum, term, upp, low, NULL);
   mpz_clears(temp, prec, width, NULL);
 }
 
@@ -278,7 +300,7 @@ void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
   fprintf(stderr, "[mpq_log_ui] decimal digits: %ld\n", digits);
 
 
-  /* calculate the series of
+  /* Compute the nth partial sum 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);
@@ -299,6 +321,35 @@ void mpq_log_ui(mpq_t rop, size_t x, size_t digits) {
     mpq_add(sum, sum, term);
   }
 
+  /* Compute the upper and lower bounds */
+  mpq_t upp, low;
+  mpq_inits(upp, low, NULL);
+
+  // upp = sum + term * y2/(1-y2)
+  mpq_mul(upp, term, y2);
+  mpq_set_ui(temp, 1, 1);
+  mpq_sub(temp, temp, y2);
+  mpq_div(upp, upp, temp);
+  mpq_add(upp, sum, upp);
+
+  // low = sum + term / ((2*n + 3)/(2*n + 1)*y^-2 - 1)
+  mpq_set_ui(temp, 2*n + 3, 2*n + 1);
+  mpq_div(temp, temp, y2);
+  mpz_sub(mpq_numref(temp), mpq_numref(temp), mpq_denref(temp));
+  mpq_canonicalize(temp);
+  mpq_div(low, term, temp);
+  mpq_add(low, sum, low);
+
+  /* Average the two bounds */
+  mpq_add(sum, low, upp);
+  mpz_div_ui(mpq_numref(sum), mpq_numref(sum), 2);
+  mpq_canonicalize(sum);
+
+  size_t d = digits > 50 ? 50 : digits;
+  mpq_fprintf(stderr, d+2, "[mpq_log_ui] upper bound: %q...\n", upp);
+  mpq_fprintf(stderr, d+2, "[mpq_log_ui] lower bound: %q...\n", low);
+  mpq_fprintf(stderr, d,   "[mpq_log_ui] series: %q...\n", sum);
+
   // multiply by 2
   mpz_mul_ui(mpq_numref(sum), mpq_numref(sum), 2);
   mpq_canonicalize(sum);