Replaced Newton's method with binary search.
authorBen Lynn <benlynn@gmail.com>
Wed, 17 Sep 2008 09:04:23 +0000 (17 02:04 -0700)
committerBen Lynn <benlynn@gmail.com>
Wed, 17 Sep 2008 09:04:23 +0000 (17 02:04 -0700)
Makefile
README
cf.h
famous_test.c
mobius.c
mobius_test.c
pi2.c [deleted file]

index c61c560..3420387 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@
 
 CF_OBJS:=cf.o mobius.o famous.o bihom.o taylor.o test.o newton.o
 TESTS:=bihom_test cf_test famous_test mobius_test newton_test
-BINS:=pi pi2 sqrt2 epow
+BINS:=pi epow
 
 target : $(BINS)
 
diff --git a/README b/README
index 99afa0d..9711093 100644 (file)
--- a/README
+++ b/README
@@ -126,15 +126,19 @@ in this code.
   for each term held back, we spawn and watch over two threads on invocation to
   service the two demand channels, storing all backlogged terms in one place.
 
-=== Increasing floating-point precision ===
-
-A disclaimer to some of my anti-floating-point rhetoric: increasing the
-precision of floating-point operations without redoing work is in fact possible
-for many common cases. For example, Newton's method is self-correcting, meaning
-we may use low precision at first, and increase it for future iterations.
-Even pi enjoys such methods: there is a formula computing any hexadecimal
-digit of pi without computing any previous digits, though it requires about the
-same amount of work.
+=== Disclaimer ===
+
+Let me put my salesmanship on hold, and temper my anti-floating-point rhetoric.
+Increasing the precision of floating-point operations without redoing work is
+in fact possible for many common cases. For example, Newton's method is
+self-correcting, meaning we may use low precision at first, and increase it for
+future iterations.  Even pi enjoys such methods: there is a formula computing
+any hexadecimal digit of pi without computing any previous digits, though it
+requires about the same amount of work.
+
+Also, Newton's method, and other floating-point algorithms converge
+quadratically or faster, thus a good implementation will asymptotically
+outperform continued fractions as these converge only linearly.
 
 However, I'd prefer continued fractions to Newton's method when writing code
 to compute digits of, say, sqrt(2). Precision decisions are made automatically:
diff --git a/cf.h b/cf.h
index f250a9c..aab54a4 100644 (file)
--- a/cf.h
+++ b/cf.h
@@ -47,7 +47,7 @@ cf_t cf_new_nonregular_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_
 // Output: Regular continued fraction. Assumes input fraction is well-behaved.
 cf_t cf_new_nonregular_to_cf(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
 // Does both of the above at once. Seems slow.
-cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
+cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a[4]);
 
 // Well-known continued fraction expansions.
 // cf_famous.c:
@@ -78,6 +78,7 @@ cf_t cf_new_cos1();
 //     a0 xy + a1 x + a2 y + a3
 // y = ------------------------
 //     a4 xy - a0 x + a5 y - a2
-cf_t cf_new_newton(cf_t x, mpz_t a[6]);
+cf_t cf_new_newton(cf_t x, mpz_t a[6], mpz_t lower);
+cf_t cf_new_sqrt(cf_t x);
 
 #endif  // __CF_H__
index 7347660..6f9b951 100644 (file)
@@ -15,5 +15,21 @@ int main() {
   CF_NEW_EXPECT_DEC(cf_new_tan1, "1.5574077246549022305");
   CF_NEW_EXPECT_DEC(cf_new_sin1, "0.8414709848078965066");
   CF_NEW_EXPECT_DEC(cf_new_cos1, "0.5403023058681397174");
+
+  mpz_t z;
+  mpz_init(z);
+  cf_t x;
+
+  mpz_set_si(z, 3);
+  x = cf_new_tanh(z);
+  CF_EXPECT_DEC(x, "0.99505475368673045133188018525548847509781385470028");
+  cf_free(x);
+
+  mpz_set_ui(z, 2);
+  x = cf_new_epow(z);
+  CF_EXPECT_DEC(x, "7.38905609893065022723042746057500781318031557055184");
+  cf_free(x);
+
+  mpz_clear(z);
   return 0;
 }
index 0d00f12..0382587 100644 (file)
--- a/mobius.c
+++ b/mobius.c
@@ -257,7 +257,7 @@ static void *nonregular_mobius_decimal(cf_t cf) {
        mpz_add(t2, t2, pq->q);
        if (mpz_cmp(t2, pq->p) > 0) {
          // Output a decimal digit.
-         cf_put(cf, pq->p);
+         cf_put(cf, t1);
          // Subtract: remainder of p/q.
          mpz_sub(t2, t2, pq->p);
          mpz_sub(t2, pq->q, t2);
@@ -288,10 +288,11 @@ static void *nonregular_mobius_decimal(cf_t cf) {
   return NULL;
 }
 
-cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
+cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a[4]) {
   mobius_data_ptr md = malloc(sizeof(*md));
   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
-  mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
+  mpz_set(md->a, a[0]); mpz_set(md->b, a[1]);
+  mpz_set(md->c, a[2]); mpz_set(md->d, a[3]);
   md->input = x;
   return cf_new(nonregular_mobius_decimal, md);
 }
index 150d840..41f46bc 100644 (file)
@@ -12,6 +12,34 @@ static void *sqrt2(cf_t cf) {
   return NULL;
 }
 
+// Converges extremely slowly.
+static void *slow_pi(cf_t cf) {
+  mpz_t num, denom, t;
+  mpz_init(num);
+  mpz_init(denom);
+  mpz_init(t);
+
+  mpz_set_ui(denom, 3);
+  cf_put(cf, denom);
+
+  mpz_set_ui(num, 1);
+  cf_put(cf, num);
+  mpz_set_ui(t, 8);
+
+  mpz_set_ui(denom, 6);
+  while(cf_wait(cf)) {
+    cf_put(cf, denom);
+    mpz_add(num, num, t);
+    cf_put(cf, num);
+    mpz_add_ui(t, t, 8);
+  }
+
+  mpz_clear(num);
+  mpz_clear(denom);
+  mpz_clear(t);
+  return NULL;
+}
+
 int main() {
   cf_t x, conv;
   x = cf_new_const(sqrt2);
@@ -46,13 +74,37 @@ int main() {
   cf_free(conv);
   cf_free(x);
 
-  x = cf_new_const(sqrt2);
-  cf_t mob;
   mpz_t z[4];
   for (int i = 0; i < 4; i++) mpz_init(z[i]);
   // Identity Mobius transformation.
   mpz_set_si(z[0], 1);
   mpz_set_si(z[3], 1);
+  mpz_t digit;
+  mpz_init(digit);
+
+  x = cf_new_const(slow_pi);
+  conv = cf_new_nonregular_mobius_to_decimal(x, z);
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 3));
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 1));
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 4));
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 1));
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 5));
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 9));
+  cf_get(digit, conv);
+  EXPECT(!mpz_cmp_ui(digit, 2));
+  
+  mpz_clear(digit);
+  cf_free(x);
+  cf_free(conv);
+
+  x = cf_new_const(sqrt2);
+  cf_t mob;
   mob = cf_new_mobius_to_cf(x, z);
   CF_EXPECT_DEC(mob, "1.4142135623730");
   cf_free(mob);
diff --git a/pi2.c b/pi2.c
deleted file mode 100644 (file)
index 501f0eb..0000000
--- a/pi2.c
+++ /dev/null
@@ -1,56 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <gmp.h>
-#include "cf.h"
-
-// Converges extremely slowly. Not worth converting to decimal or
-// continued fraction form.
-static void *pi_six_sequence(cf_t cf) {
-  mpz_t num, denom, t;
-  mpz_init(num);
-  mpz_init(denom);
-  mpz_init(t);
-
-  mpz_set_ui(denom, 3);
-  cf_put(cf, denom);
-
-  mpz_set_ui(num, 1);
-  cf_put(cf, num);
-  mpz_set_ui(t, 8);
-
-  mpz_set_ui(denom, 6);
-  while(cf_wait(cf)) {
-    cf_put(cf, denom);
-    mpz_add(num, num, t);
-    cf_put(cf, num);
-    mpz_add_ui(t, t, 8);
-  }
-
-  mpz_clear(num);
-  mpz_clear(denom);
-  mpz_clear(t);
-  return NULL;
-}
-
-int main() {
-  mpz_t z;
-  mpz_t a, b, c, d;
-  mpz_init(z);
-  cf_t pi, conv;
-  pi = cf_new(pi_six_sequence, NULL);
-  mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
-  mpz_set_ui(a, 1); mpz_set_ui(b, 0);
-  mpz_set_ui(c, 0); mpz_set_ui(d, 1);
-
-  conv = cf_new_nonregular_mobius_to_decimal(pi, a, b, c, d);
-  for (int i = 1; i <= 5000; i++) {
-    cf_get(z, conv);
-    gmp_printf("%Zd\n",z);
-  }
-  printf("\n");
-  cf_free(conv);
-  cf_free(pi);
-  mpz_clear(z);
-  mpz_clear(d); mpz_clear(c); mpz_clear(b); mpz_clear(a);
-  return 0;
-}