Improved Makefile, tests.
authorBen Lynn <benlynn@gmail.com>
Tue, 16 Sep 2008 11:20:08 +0000 (16 04:20 -0700)
committerBen Lynn <benlynn@gmail.com>
Tue, 16 Sep 2008 11:20:08 +0000 (16 04:20 -0700)
Makefile
cf.h
cf_test.c
famous.c
mobius.c
mobius_test.c [new file with mode: 0644]
sqrt2.c [deleted file]

dissimilarity index 99%
index 5f61db1..bc49a6c 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,13 +1,22 @@
-CF_FILES:=cf.c mobius.c famous.c bihom.c taylor.c
-TEST_FILES:=$(CF_FILES) test.c
-
-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 $(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
+.PHONY: test target clean
+
+CF_OBJS:=cf.o mobius.o famous.o bihom.o taylor.o test.o
+TESTS:=bihom_test cf_test famous_test mobius_test
+BINS:=pi pi2 sqrt2 epow
+
+target : $(BINS)
+
+libfrac.a : $(CF_OBJS)
+       ar rvs $@ $^
+
+%.o : %.c
+       gcc -g -std=c99 -Wall -c -o $@ $<
+
+% : %.c libfrac.a
+       #gcc -O3 -std=c99 -Wall -o $@ $< -lgmp -lpthread -lfrac -L .
+       gcc -g -std=c99 -Wall -o $@ $< -lgmp -lpthread -lfrac -L .
+
+test: $(TESTS)
+
+clean:
+       -rm *.o $(TESTS) $(BINS) libfrac.a
diff --git a/cf.h b/cf.h
index c319414..c607946 100644 (file)
--- a/cf.h
+++ b/cf.h
@@ -26,27 +26,33 @@ int cf_wait(cf_t cf);
 void *cf_data(cf_t cf);
 
 // cf_mobius.c:
+
 // Compute convergents of a simple continued fraction x.
 // Outputs p then q on channel, where p/q is the last convergent computed.
-cf_t cf_new_convergent(cf_t x);
+cf_t cf_new_cf_convergent(cf_t x);
+// Compute decimal representation of a simple continued fraction x.
+// Outputs integer part first, then digits one at a time.
+cf_t cf_new_cf_to_decimal(cf_t x);
 
 // Compute convergents of (a x + b)/(c x + d)
 // where x is a regular continued fraction.
 cf_t cf_new_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
+cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
+cf_t cf_new_mobius_to_cf(cf_t x, mpz_t z[4]);
 
 // Compute convergents of (a x + b)/(c x + d)
 // where x is a nonregular continued fraction.
 cf_t cf_new_nonregular_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
-
+// Input: Mobius transformation and nonregular continued fraction.
+// 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);
-
-cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
-cf_t cf_new_cf_to_decimal(cf_t x);
+// 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);
 
 // Well-known continued fraction expansions.
 // cf_famous.c:
 // e:
+cf_t cf_new_sqrt2();
 cf_t cf_new_e();
 cf_t cf_new_pi();
 cf_t cf_new_tan1();
index c83694b..c417a3e 100644 (file)
--- a/cf_test.c
+++ b/cf_test.c
@@ -22,7 +22,7 @@ int main() {
   mpz_t z;
   mpz_init(z);
   cf_t a;
-  a = cf_new(count_fn, NULL);
+  a = cf_new_const(count_fn);
   for (int i = 0; i < 100; i++) {
     cf_get(z, a);
     EXPECT(!mpz_cmp_ui(z, i));
index 9e6e06b..48e279f 100644 (file)
--- a/famous.c
+++ b/famous.c
@@ -3,6 +3,19 @@
 #include <gmp.h>
 #include "cf.h"
 
+// sqrt(2) = [1; 2, 2, ...]
+static void *sqrt2(cf_t cf) {
+  cf_put_int(cf, 1);
+  while(cf_wait(cf)) {
+    cf_put_int(cf, 2);
+  }
+  return NULL;
+}
+
+cf_t cf_new_sqrt2() {
+  return cf_new_const(sqrt2);
+}
+
 // e = [2; 1, 2, 1, 1, 4, 1, ...]
 static void *e_expansion(cf_t cf) {
   int even = 2;
index 8507449..0d00f12 100644 (file)
--- a/mobius.c
+++ b/mobius.c
@@ -122,7 +122,7 @@ cf_t cf_new_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 
 // Start a thread that, when signalled, computes the convergents of a continued
 // fraction.
-cf_t cf_new_convergent(cf_t x) {
+cf_t cf_new_cf_convergent(cf_t x) {
   mpz_t one, zero;
   mpz_init(one); mpz_init(zero);
   mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
@@ -173,6 +173,8 @@ cf_t cf_new_nonregular_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_
   return cf_new(nonregular_mobius_convergent, md);
 }
 
+// Input: Mobius transformation and nonregular continued fraction.
+// Output: Regular continued fraction. Assumes input fraction is well-behaved.
 static void *mobius_nonregular_throughput(cf_t cf) {
   mobius_data_ptr md = cf_data(cf);
   cf_t input = md->input;
@@ -234,14 +236,67 @@ cf_t cf_new_nonregular_to_cf(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
   return cf_new(mobius_nonregular_throughput, md);
 }
 
-static void *mobius_decimal(cf_t cf) {
+// This seems to be slower than regularizing the continued fraction
+// and then converting to decimal.
+static void *nonregular_mobius_decimal(cf_t cf) {
   mobius_data_ptr md = cf_data(cf);
   cf_t input = md->input;
   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
+  mpz_t num; mpz_init(num);
   mpz_t denom; mpz_init(denom);
   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
+  int recur() {
+    pqset_nonregular_recur(pq, num, denom);
+    pqset_remove_gcd(pq, t0, t1);
+
+    // If the denominator is zero, we can't do anything yet.
+    if (mpz_sgn(pq->qold)) {
+      mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
+      mpz_mul(t2, t1, pq->q);
+      if (mpz_cmp(t2, pq->p) <= 0) {
+       mpz_add(t2, t2, pq->q);
+       if (mpz_cmp(t2, pq->p) > 0) {
+         // Output a decimal digit.
+         cf_put(cf, pq->p);
+         // Subtract: remainder of p/q.
+         mpz_sub(t2, t2, pq->p);
+         mpz_sub(t2, pq->q, t2);
+         // Multiply numerator by 10.
+         mpz_mul_ui(pq->pold, t0, 10);
+         mpz_mul_ui(pq->p, t2, 10);
+         return 1;
+       }
+      }
+    }
+    return 0;
+  }
+  mpz_set_ui(num, 1);
+  cf_get(denom, input);
+  recur();
+  while(cf_wait(cf)) {
+    do {
+      cf_get(num, input);
+      cf_get(denom, input);
+    } while(!recur());
+  }
+  mpz_clear(num);
+  mpz_clear(denom);
+  pqset_clear(pq);
+  mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
+  mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
+  free(md);
+  return NULL;
+}
 
-  // Determine the sign.
+cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
+  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);
+  md->input = x;
+  return cf_new(nonregular_mobius_decimal, md);
+}
+
+static void determine_sign(cf_t cf, pqset_t pq, mpz_t denom, cf_t input) {
   cf_set_sign(cf, cf_sign(input));
   while (mpz_sgn(pq->pold) != mpz_sgn(pq->p)
       || mpz_sgn(pq->qold) != mpz_sgn(pq->q)) {
@@ -258,40 +313,39 @@ static void *mobius_decimal(cf_t cf) {
     mpz_neg(pq->p, pq->p);
     cf_flip_sign(cf);
   }
+}
+
+// Input: Mobius transformation and regular continued fraction.
+// Output: Regular continued fraction.
+static void *mobius_throughput(cf_t cf) {
+  mobius_data_ptr md = cf_data(cf);
+  cf_t input = md->input;
+  pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
+  mpz_t denom; mpz_init(denom);
+  mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
+
+  determine_sign(cf, pq, denom, input);
 
   int recur() {
     pqset_regular_recur(pq, denom);
 
-    // If the denominator is zero, we can't do anything yet.
     if (mpz_sgn(pq->qold)) {
-      // Each term except possibly the first is one of {0, ..., 9}.
-      /* Naive attempt to expoit this didn't seem faster:
-       * (and I'd have to handle the first term properly)
-      int i;
-      mpz_set(t0, pq->q);
-      for (i = 0; i <= 9; i++) {
-       if (mpz_cmp(t0, pq->p) > 0) break;
-       mpz_add(t0, t0, pq->q);
-      }
-      mpz_set_ui(pq->pold, i);
-      mpz_sub(t0, t0, pq->p);
-      mpz_sub(t0, pq->q, t0);
-      mpz_mul(pq->qold, pq->pold, pq->qnew);
-      */
-
       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
       mpz_mul(t2, t1, pq->q);
+
       if (mpz_cmp(t2, pq->p) <= 0) {
        mpz_add(t2, t2, pq->q);
        if (mpz_cmp(t2, pq->p) > 0) {
-         // Output a decimal digit.
+         // Output continued fraction term.
          cf_put(cf, t1);
-         // Compute t2 = remainder of p/q.
+         // Subtract: remainder of p/q.
          mpz_sub(t2, t2, pq->p);
          mpz_sub(t2, pq->q, t2);
-         // Multiply numerator by 10.
-         mpz_mul_ui(pq->pold, t0, 10);
-         mpz_mul_ui(pq->p, t2, 10);
+         // Invert
+         mpz_set(pq->pold, pq->qold);
+         mpz_set(pq->qold, t0);
+         mpz_set(pq->p, pq->q);
+         mpz_set(pq->q, t2);
          return 1;
        }
       }
@@ -310,46 +364,55 @@ static void *mobius_decimal(cf_t cf) {
   free(md);
   return NULL;
 }
-cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
+
+cf_t cf_new_mobius_to_cf(cf_t x, mpz_t z[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, z[0]); mpz_set(md->b, z[1]);
+  mpz_set(md->c, z[2]); mpz_set(md->d, z[3]);
   md->input = x;
-  return cf_new(mobius_decimal, md);
-}
-
-cf_t cf_new_cf_to_decimal(cf_t x) {
-  mpz_t one, zero;
-  mpz_init(one); mpz_init(zero);
-  mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
-  cf_t res = cf_new_mobius_to_decimal(x, one, zero, zero, one);
-  mpz_clear(one); mpz_clear(zero);
-  return res;
+  return cf_new(mobius_throughput, md);
 }
 
-// This seems to be slower than regularizing the continued fraction
-// and then converting to decimal.
-static void *nonregular_mobius_decimal(cf_t cf) {
+// Input: Mobius transformation and regular continued fraction.
+// Output: Decimal representation. The integer part is given first,
+// followed by one digit at a time.
+static void *mobius_decimal(cf_t cf) {
   mobius_data_ptr md = cf_data(cf);
   cf_t input = md->input;
   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
-  mpz_t num; mpz_init(num);
   mpz_t denom; mpz_init(denom);
   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
+
+  determine_sign(cf, pq, denom, input);
   int recur() {
-    pqset_nonregular_recur(pq, num, denom);
-    pqset_remove_gcd(pq, t0, t1);
+    pqset_regular_recur(pq, denom);
 
     // If the denominator is zero, we can't do anything yet.
     if (mpz_sgn(pq->qold)) {
+      // Each term except possibly the first is one of {0, ..., 9}.
+      /* Naive attempt to expoit this didn't seem faster:
+       * (and I'd have to handle the first term properly)
+      int i;
+      mpz_set(t0, pq->q);
+      for (i = 0; i <= 9; i++) {
+       if (mpz_cmp(t0, pq->p) > 0) break;
+       mpz_add(t0, t0, pq->q);
+      }
+      mpz_set_ui(pq->pold, i);
+      mpz_sub(t0, t0, pq->p);
+      mpz_sub(t0, pq->q, t0);
+      mpz_mul(pq->qold, pq->pold, pq->qnew);
+      */
+
       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
       mpz_mul(t2, t1, pq->q);
       if (mpz_cmp(t2, pq->p) <= 0) {
        mpz_add(t2, t2, pq->q);
        if (mpz_cmp(t2, pq->p) > 0) {
          // Output a decimal digit.
-         cf_put(cf, pq->p);
-         // Subtract: remainder of p/q.
+         cf_put(cf, t1);
+         // Compute t2 = remainder of p/q.
          mpz_sub(t2, t2, pq->p);
          mpz_sub(t2, pq->q, t2);
          // Multiply numerator by 10.
@@ -361,16 +424,11 @@ static void *nonregular_mobius_decimal(cf_t cf) {
     }
     return 0;
   }
-  mpz_set_ui(num, 1);
-  cf_get(denom, input);
-  recur();
   while(cf_wait(cf)) {
     do {
-      cf_get(num, input);
       cf_get(denom, input);
     } while(!recur());
   }
-  mpz_clear(num);
   mpz_clear(denom);
   pqset_clear(pq);
   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
@@ -379,10 +437,19 @@ 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_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
   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);
   md->input = x;
-  return cf_new(nonregular_mobius_decimal, md);
+  return cf_new(mobius_decimal, md);
+}
+
+cf_t cf_new_cf_to_decimal(cf_t x) {
+  mpz_t one, zero;
+  mpz_init(one); mpz_init(zero);
+  mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
+  cf_t res = cf_new_mobius_to_decimal(x, one, zero, zero, one);
+  mpz_clear(one); mpz_clear(zero);
+  return res;
 }
diff --git a/mobius_test.c b/mobius_test.c
new file mode 100644 (file)
index 0000000..150d840
--- /dev/null
@@ -0,0 +1,74 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <gmp.h>
+#include "cf.h"
+#include "test.h"
+
+static void *sqrt2(cf_t cf) {
+  cf_put_int(cf, 1);
+  while(cf_wait(cf)) {
+    cf_put_int(cf, 2);
+  }
+  return NULL;
+}
+
+int main() {
+  cf_t x, conv;
+  x = cf_new_const(sqrt2);
+  conv = cf_new_cf_convergent(x);
+
+  mpz_t p, q;
+  mpz_init(p);
+  mpz_init(q);
+
+  cf_get(p, conv);
+  EXPECT(!mpz_cmp_ui(p, 1));
+  cf_get(q, conv);
+  EXPECT(!mpz_cmp_ui(q, 1));
+
+  cf_get(p, conv);
+  EXPECT(!mpz_cmp_ui(p, 3));
+  cf_get(q, conv);
+  EXPECT(!mpz_cmp_ui(q, 2));
+
+  cf_get(p, conv);
+  EXPECT(!mpz_cmp_ui(p, 7));
+  cf_get(q, conv);
+  EXPECT(!mpz_cmp_ui(q, 5));
+
+  cf_get(p, conv);
+  EXPECT(!mpz_cmp_ui(p, 17));
+  cf_get(q, conv);
+  EXPECT(!mpz_cmp_ui(q, 12));
+
+  mpz_clear(p);
+  mpz_clear(q);
+  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);
+  mob = cf_new_mobius_to_cf(x, z);
+  CF_EXPECT_DEC(mob, "1.4142135623730");
+  cf_free(mob);
+  cf_free(x);
+  x = cf_new_const(sqrt2);
+  // (x - 2)/(-3x + 4)
+  // Works out to be 1 + sqrt(2).
+  mpz_set_si(z[0], 1);
+  mpz_set_si(z[1], -2);
+  mpz_set_si(z[2], -3);
+  mpz_set_si(z[3], 4);
+  mob = cf_new_mobius_to_cf(x, z);
+  CF_EXPECT_DEC(mob, "2.4142135623730");
+  cf_free(x);
+  cf_free(mob);
+  for (int i = 0; i < 4; i++) mpz_clear(z[i]);
+
+  return 0;
+}
diff --git a/sqrt2.c b/sqrt2.c
deleted file mode 100644 (file)
index 6109ff4..0000000
--- a/sqrt2.c
+++ /dev/null
@@ -1,54 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <gmp.h>
-#include "cf.h"
-
-static void *sqrt2(cf_t cf) {
-  cf_put_int(cf, 1);
-  while(cf_wait(cf)) {
-    cf_put_int(cf, 2);
-  }
-  return NULL;
-}
-
-cf_t cf_new_sqrt2() {
-  return cf_new_const(sqrt2);
-}
-
-int main() {
-  cf_t x, conv;
-  x = cf_new_sqrt2();
-  conv = cf_new_convergent(x);
-
-  mpz_t p, q;
-  mpz_init(p);
-  mpz_init(q);
-  for (int i = 0; i < 10; i++) {
-    cf_get(p, conv);
-    cf_get(q, conv);
-    gmp_printf("p/q = %Zd/%Zd\n", p, q);
-  }
-  cf_free(conv);
-  cf_free(x);
-
-  x = cf_new_sqrt2();
-  mpz_t a, b, c, d;
-  mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
-  mpz_set_si(a, -7);
-  mpz_set_si(b, 5);
-  mpz_set_si(c, -12);
-  mpz_set_si(d, 13);
-  cf_t dec = cf_new_mobius_to_decimal(x, a, b, c, d);
-  cf_get(p, dec);
-  putchar(cf_sign(dec) > 0 ? ' ' : '-');
-  gmp_printf("%Zd.", p);
-  for (int i = 0; i < 10; i++) {
-    cf_get(p, dec);
-    gmp_printf("%Zd", p);
-  }
-  putchar('\n');
-  cf_free(x);
-  mpz_clear(p);
-  mpz_clear(q);
-  return 0;
-}