ex-2/fast: implement proper exit condition for the algorithm

ex-2/fast.c

@@ -17,7 +17,12 @@
  *
  * If N = 2^p (to simplify the log) to get the result to D
  * decimal digits we must have 2^(p+1) > D / log(π), which is
- * satisfied by p = floor(log(D)/log(2)).
+ * satisfied by p = floor(log2(N log(10) + log(π))) - 1.
+ *
+ * The number of series terms needed to reduce the error, due to
+ * to less than 10^-D is α⋅2^p where α is the solution of:
+ *
+ *   α log(α) = 3 + α  ⇒  α = exp(W(3e) - 1) ≈ 4.97
  *
  * source: http://www.numberworld.org/y-cruncher/internals/formulas.html
  *
@@ -104,7 +109,7 @@ void mpf_log2(mpf_t rop, size_t digits) {
   /* Note: printf() rounds floating points to the number of digits
    * in an implementation-dependent manner. To avoid printing
    * (wrongly) rounded digits we increase the count by 1 and hide
-   * the last by overwiting with a space.
+   * the last by overwriting with a space.
    */
   size_t d = digits>50? 50 : digits;
   gmp_fprintf(stderr, "[mpf_log2] upper bound: %.*Ff\b \n", d+3, upp);
@@ -120,20 +125,37 @@ void mpf_log2(mpf_t rop, size_t digits) {
 }
 
 
-void mpf_gamma(mpf_t rop, size_t digits, size_t prec) {
+/* `mpf_gamma(rop, digits)` computes γ to `digits` decimal places
-  unsigned int p = log2(digits);
+ * using the above described algorithm and set results to `rop`.
+ */
+void mpf_gamma(mpf_t rop, size_t digits) {
+  /* The number N=2^p is fixed by requiring that the error
+   * is smaller than the precision/digits required:
+   *
+   *   π exp(-4N) < 10^(-D)
+   *
+   * The first power of 2 solving the constraint is
+   *
+   *   p = floor(log2(N log(10) + log(π))) - 1
+   */
+  unsigned int p = log2(digits * log(10.0) + log(M_PI)) - 1;
   fprintf(stderr, "[mpf_gamma] p: %d\n", p);
+  fprintf(stderr, "[mpf_gamma] log2(digits): %d\n", (int)log2(digits));
+
+  /* The number of iterations needed to reduce the error, due to
+   * ignoring the C/B² term, to less than 10^-D is α⋅2^p where α
+   * is the solution of:
+   *
+   *   α log(α) = 3 + α  ⇒  α = exp(W(3e) - 1) ≈ 4.97
+   *
+   */
+  size_t kmax = 4.97 * (1 << p);
+  fprintf(stderr, "[mpf_gamma] iterations: %ld\n", kmax);
 
   // initialise variables
   mpf_t A, B, a, b, stop, temp;
   mpf_inits(A, B, a, b, stop, temp, NULL);
 
-  /* Smallest representable number */
-  mpf_set_ui(stop, 10);
-  mpf_pow_ui(stop, stop, digits);
-  mpf_ui_div(stop, 1, stop);
-  gmp_fprintf(stderr, "[mpf_gamma] stop=%Fe\n", stop);
-
   // A = a0 = -log(2^p) = -plog(2)
   mpf_log2(a, digits);
   mpf_mul_ui(a, a, p);
@@ -144,8 +166,7 @@ void mpf_gamma(mpf_t rop, size_t digits, size_t prec) {
   mpf_set_si(b, 1);
   mpf_set_si(B, 1);
 
-  size_t k;
+  for(size_t k = 1; k <= kmax; k++) {
-  for(k = 1;; k++) {
     // bk = bk-1 (n/k)^2
     mpf_mul_2exp(b, b, 2*p);
     mpf_div_ui(b, b, k);
@@ -163,23 +184,13 @@ void mpf_gamma(mpf_t rop, size_t digits, size_t prec) {
     // A += ak, B += bk
     mpf_add(A, A, a);
     mpf_add(B, B, b);
-
-    /* Exit when both ak and bk are smaller
-     * than the asked precision (10^-D).
-     */
-    if (mpf_cmp(a, b) >= 0) {
-      if (mpf_cmp(a, stop) < 0) break;
-    } else {
-      if (mpf_cmp(b, stop) < 0) break;
   }
-  }
-  fprintf(stderr, "[mpf_gamma] iterations: %ld\n", k);
 
   // return A/B
   mpf_div(rop, A, B);
 
   // free memory
-  mpf_clears(A, B, a, b, temp, stop, NULL);
+  mpf_clears(A, B, a, b, temp, NULL);
 }
 
 int main(int argc, char** argv) {
@@ -205,7 +216,7 @@ int main(int argc, char** argv) {
   mpf_set_default_prec(prec + 64);
 
   mpf_t y; mpf_init(y);
-  mpf_gamma(y, digits, prec);
+  mpf_gamma(y, digits);
 
   // again, trick to avoid rounding
   gmp_printf("%.*Ff\b \n", digits+1, y);