Improved tests. Appears to be a bug in bihom.c.
authorBen Lynn <benlynn@gmail.com>
Mon, 15 Sep 2008 11:34:01 +0000 (15 04:34 -0700)
committerBen Lynn <benlynn@gmail.com>
Mon, 15 Sep 2008 11:34:01 +0000 (15 04:34 -0700)
Makefile
README
cf.h
famous_test.c
taylor.c
test.c [new file with mode: 0644]
test.h

index 16338a1..5f61db1 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,13 @@
-CF_FILES:=cf.c mobius.c famous.c bihom.c
+CF_FILES:=cf.c mobius.c famous.c bihom.c taylor.c
+TEST_FILES:=$(CF_FILES) test.c
 
-target : pi pi2 sqrt2 epow taylor
+target : pi pi2 sqrt2 epow
 
 % : %.c $(CF_FILES)
        #gcc -O3 -std=c99 -Wall -o $@ $^ -lgmp -lpthread
        gcc -g -std=c99 -Wall -o $@ $^ -lgmp -lpthread
 
 test :
-       gcc -g -std=c99 -Wall -o cf_test cf_test.c $(CF_FILES) -lgmp -lpthread
-       gcc -g -std=c99 -Wall -o famous_test famous_test.c $(CF_FILES) -lgmp -lpthread
+       gcc -g -std=c99 -Wall -o cf_test cf_test.c $(TEST_FILES) -lgmp -lpthread
+       gcc -g -std=c99 -Wall -o famous_test famous_test.c $(TEST_FILES) -lgmp -lpthread
+       gcc -g -std=c99 -Wall -o bihom_test bihom_test.c $(TEST_FILES) -lgmp -lpthread
diff --git a/README b/README
index b9d6f93..e0db81d 100644 (file)
--- a/README
+++ b/README
@@ -30,7 +30,9 @@ If only we had used continued fractions, for then we could arbitrarily
 increase the precision of our answer without redoing any work. Furthermore,
 the error analysis becomes trivial.
 
-However, we must pay a price for the advantages of continued fractions.  We
+== Disadvantages ==
+
+We must pay a price for the advantages of continued fractions.  We
 have a representation from which we can instantly infer error bounds. A
 representation for which we may arbitrarily increase precision without redoing
 any work. A representation which in some sense is the closest we can hope to
@@ -49,12 +51,14 @@ approximations are coarser but the show goes on.
 
 I use the GMP library for multi-precision arithmetic. Type
  
- $ make
+  $ make
+
+to compile. Try "pi 2000" to generate 2000 digits of pi using a continued
+fraction. It takes longer than I had hoped, but it's faster than
 
-to compile. The "pi" binary generates 5000 digits of pi using a continued
-fraction. It takes longer than I had hoped.
+  $ echo "scale = 2000; 4*a(1)" | bc -l
 
-There's no API documentation yet; see the source.
+There's no API documentation; see the source.
 
 == Demand channels ==
 
@@ -118,4 +122,4 @@ in this code.
 - I'd like a split function as described by McIlroy, though for continued
   fractions a "tee"-like function is more apt: rather than spawning a thread
   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 thread.
+  service the two demand channels, storing all backlogged terms in one place.
diff --git a/cf.h b/cf.h
index 07f223a..13b379c 100644 (file)
--- a/cf.h
+++ b/cf.h
@@ -59,4 +59,8 @@ cf_t cf_new_tan(mpz_t z);
 // Gosper's method for computing bihomographic functions of continued fractions.
 cf_t cf_new_bihom(cf_t x, cf_t y, mpz_t a[8]);
 
+// From taylor.c:
+cf_t cf_new_sin1();
+cf_t cf_new_cos1();
+
 #endif  // __CF_H__
dissimilarity index 84%
index 5e7adb9..ec2ac79 100644 (file)
@@ -1,59 +1,19 @@
-// Test famous continued fraction expansions are correct.
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <gmp.h>
-#include "cf.h"
-#include "test.h"
-
-void check_match(cf_t (*cf_new_fn)(), char *result) {
-  mpz_t z;
-  mpz_init(z);
-  cf_t x, conv;
-  x = cf_new_fn();
-  conv = cf_new_cf_to_decimal(x);
-  int len = strlen(result);
-  char s[len + 1];
-  for (int i = 0; i <= len; i++) {
-    cf_get(z, conv);
-    s[i] = mpz_get_ui(z) + '0';
-  }
-  s[len] = 0;
-  EXPECT(!strcmp(s, result));
-
-  cf_free(conv);
-  cf_free(x);
-  mpz_clear(z);
-}
-
-int main() {
-  check_match(cf_new_e, "2718281828");
-  check_match(cf_new_pi, "31415926535897932384");
-  check_match(cf_new_tan1, "15574077246549022305");
-
-  mpz_t z;
-  mpz_init(z);
-  cf_t e = cf_new_e();
-  cf_t pi = cf_new_pi();
-  mpz_t addarray[8];
-  for (int i = 0; i < 8; i++) {
-    mpz_init(addarray[i]);
-    mpz_set_ui(addarray[i], 0);
-  }
-  mpz_set_ui(addarray[1], 1);
-  mpz_set_ui(addarray[2], 1);
-  mpz_set_ui(addarray[7], 1);
-  cf_t b = cf_new_bihom(e, pi, addarray);
-  cf_t conv = cf_new_cf_to_decimal(b);
-  cf_get(z, conv);
-  gmp_printf("e + pi = %Zd.", z);
-  for (int i = 0; i <= 20; i++) {
-    cf_get(z, conv);
-    gmp_printf("%Zd", z);
-  }
-  printf("\n");
-
-  mpz_clear(z);
-  return 0;
-}
+// Test famous continued fraction expansions are correct.
+// as well as expansions based on Taylor series.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+#include "test.h"
+
+int main() {
+  CF_NEW_EXPECT_DEC(cf_new_e, "2718281828");
+  CF_NEW_EXPECT_DEC(cf_new_pi, "31415926535897932384");
+  CF_NEW_EXPECT_DEC(cf_new_tan1, "15574077246549022305");
+
+  CF_NEW_EXPECT_DEC(cf_new_sin1, "08414709848078965066");
+  CF_NEW_EXPECT_DEC(cf_new_cos1, "05403023058681397174");
+  return 0;
+}
index a834a48..838b793 100644 (file)
--- a/taylor.c
+++ b/taylor.c
@@ -33,7 +33,7 @@ static void *sin1_expansion(cf_t cf) {
              // numbered from 0, unlike continued fraction terms.
   int i = 2; // Number of convergent we're computing.
 
-  void increase_limit() {
+  void more_taylor() {
     k++;
     mpz_mul_ui(mpq_numref(r1), mpq_numref(r0), (2 * k) * (2 * k + 1));
     mpz_mul_ui(mpq_denref(r1), mpq_denref(r0), (2 * k) * (2 * k + 1));
@@ -45,7 +45,6 @@ static void *sin1_expansion(cf_t cf) {
   }
   while (cf_wait(cf)) {
     for(;;) {
-      printf("find x\n");
       // Find largest x such that (p1 x + p0)/(q1 x + q0)
       // <= r0 for odd i, and >= r1 for even i
       //  =>  x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
@@ -62,7 +61,7 @@ static void *sin1_expansion(cf_t cf) {
       // Then x is the next convergent provided substituting x + 1 overshoots
       // the other bound. (If not, we need better bounds to decide.)
       if (!mpz_sgn(t0)) {
-       increase_limit();
+       more_taylor();
       } else {
        mpz_add_ui(t1, t0, 1);
        mpz_set(mpq_numref(tq0), p0);
@@ -76,7 +75,7 @@ static void *sin1_expansion(cf_t cf) {
          if (mpq_cmp(tq0, r0) < 0) break;
          // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
        }
-       increase_limit();
+       more_taylor();
       }
     };
     cf_put(cf, t0);
@@ -94,14 +93,100 @@ static void *sin1_expansion(cf_t cf) {
   return NULL;
 }
 
-int main(int argc, char *argv[]) {
-  cf_t x = cf_new_const(sin1_expansion);
-  cf_t conv = cf_new_convergent(x);
-  mpz_t z;
-  mpz_init(z);
-  for (int i = 0; i < 20; i++) {
-    cf_get(z, conv); gmp_printf("%d: %Zd/", i, z);
-    cf_get(z, conv); gmp_printf("%Zd\n", z);
+cf_t cf_new_sin1() {
+  return cf_new_const(sin1_expansion);
+}
+
+// cos 1 from Taylor series
+static void *cos1_expansion(cf_t cf) {
+  mpz_t p0, q0;  // p0/q0 = last convergent
+  mpz_t p1, q1;  // p1/q1 = convergent
+  mpq_t r0;      // lower bound for sin 1
+  mpq_t r1;      // upper bound for sin 1
+  mpz_t t0, t1, t2;  // Temporary variables.
+  mpq_t tq0;
+
+  mpq_init(r0); mpq_init(r1);
+  mpz_init(p0); mpz_init(p1); mpz_init(q0); mpz_init(q1);
+  mpz_init(t0); mpz_init(t1); mpz_init(t2);
+  mpq_init(tq0);
+
+  // Zero causes all sorts of problems.
+  // Output it and forget about it.
+  cf_put_int(cf, 0);
+
+  // Initialize convergents.
+  mpz_set_ui(p0, 1);
+  mpz_set_ui(q1, 1);
+
+  // cos 1 = 1 - 1/2 + ...
+  mpq_set_ui(r0, 1, 2);
+  mpq_set_ui(r1, 1, 1);
+
+  int k = 1; // Number of Taylor series terms we've computed,
+             // numbered from 0, unlike continued fraction terms.
+  int i = 2; // Number of convergent we're computing.
+
+  void more_taylor() {
+    k++;
+    mpz_mul_ui(mpq_numref(r1), mpq_numref(r0), (2 * k) * (2 * k - 1));
+    mpz_mul_ui(mpq_denref(r1), mpq_denref(r0), (2 * k) * (2 * k - 1));
+    mpz_add_ui(mpq_numref(r1), mpq_numref(r1), 1);
+    k++;
+    mpz_mul_ui(mpq_numref(r0), mpq_numref(r1), (2 * k) * (2 * k - 1));
+    mpz_mul_ui(mpq_denref(r0), mpq_denref(r1), (2 * k) * (2 * k - 1));
+    mpz_sub_ui(mpq_numref(r0), mpq_numref(r0), 1);
+  }
+  while (cf_wait(cf)) {
+    for(;;) {
+      // Find largest x such that (p1 x + p0)/(q1 x + q0)
+      // <= r0 for odd i, and >= r1 for even i
+      //  =>  x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
+      // for the appropriate r.
+      mpq_ptr r = i & 1 ? r0 : r1;
+      mpz_mul(t0, q0, mpq_numref(r));
+      mpz_mul(t1, p0, mpq_denref(r));
+      mpz_sub(t2, t0, t1);
+      mpz_mul(t0, p1, mpq_denref(r));
+      mpz_mul(t1, q1, mpq_numref(r));
+      mpz_sub(t1, t0, t1);
+      mpz_fdiv_q(t0, t2, t1);
+
+      // Then x is the next convergent provided substituting x + 1 overshoots
+      // the other bound. (If not, we need better bounds to decide.)
+      if (!mpz_sgn(t0)) {
+       more_taylor();
+      } else {
+       mpz_add_ui(t1, t0, 1);
+       mpz_set(mpq_numref(tq0), p0);
+       mpz_addmul(mpq_numref(tq0), p1, t1);
+       mpz_set(mpq_denref(tq0), q0);
+       mpz_addmul(mpq_denref(tq0), q1, t1);
+       if (i & 1) {
+         if (mpq_cmp(tq0, r1) > 0) break;
+         // Check (p1(x+1) + p0)/(q1(x+1) + q0) > r1
+       } else {
+         if (mpq_cmp(tq0, r0) < 0) break;
+         // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
+       }
+       more_taylor();
+      }
+    };
+    cf_put(cf, t0);
+    i++;
+    mpz_addmul(p0, p1, t0);
+    mpz_addmul(q0, q1, t0);
+    mpz_swap(p1, p0);
+    mpz_swap(q1, q0);
   }
-  return 0;
+
+  mpq_clear(r0); mpq_clear(r1);
+  mpz_clear(p0); mpz_clear(p1); mpz_clear(q0); mpz_clear(q1);
+  mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
+  mpq_clear(tq0);
+  return NULL;
+}
+
+cf_t cf_new_cos1() {
+  return cf_new_const(cos1_expansion);
 }
diff --git a/test.c b/test.c
new file mode 100644 (file)
index 0000000..075c714
--- /dev/null
+++ b/test.c
@@ -0,0 +1,35 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+
+void cf_expect_dec(cf_t x, char *result, char *filename, int line) {
+  mpz_t z;
+  mpz_init(z);
+  cf_t conv;
+  conv = cf_new_cf_to_decimal(x);
+  int len = strlen(result);
+  char s[len + 1];
+  for (int i = 0; i <= len; i++) {
+    cf_get(z, conv);
+    s[i] = mpz_get_ui(z) + '0';
+  }
+  s[len] = 0;
+  if (strcmp(s, result)) {
+    fprintf(stderr, "\n%s:%d: bad continued fraction decimal expansion\n",
+        filename, line);
+    fprintf(stderr, "  expected: %s\n", result);
+    fprintf(stderr, "    actual: %s\n\n", s);
+  }
+
+  cf_free(conv);
+  mpz_clear(z);
+}
+
+void cf_new_expect_dec(cf_t (*cf_new_fn)(), char *result,
+    char *filename, int line) {
+  cf_t x = cf_new_fn();
+  cf_expect_dec(x, result, filename, line);
+  cf_free(x);
+}
diff --git a/test.h b/test.h
index ca35a77..dd3cf77 100644 (file)
--- a/test.h
+++ b/test.h
@@ -1,2 +1,13 @@
 #define EXPECT(condition) \
   if (condition); else fprintf(stderr, "%s:%d: FAIL\n", __FILE__, __LINE__)
+
+#define CF_EXPECT_DEC(cf, str) \
+  cf_expect_dec(cf, str, __FILE__, __LINE__)
+
+#define CF_NEW_EXPECT_DEC(cf, str) \
+  cf_new_expect_dec(cf, str, __FILE__, __LINE__)
+
+void cf_expect_dec(cf_t x, char *result, char *filename, int line);
+
+void cf_new_expect_dec(cf_t (*cf_new_fn)(), char *result,
+    char *filename, int line);