Implemented tee.
authorBen Lynn <benlynn@gmail.com>
Wed, 17 Sep 2008 11:36:10 +0000 (17 04:36 -0700)
committerBen Lynn <benlynn@gmail.com>
Wed, 17 Sep 2008 11:36:10 +0000 (17 04:36 -0700)
Also realized I had forgotten to add newton* files.

Makefile
cf.c
cf.h
newton.c [new file with mode: 0644]
newton_test.c [new file with mode: 0644]
tee.c [new file with mode: 0644]
tee_test.c [new file with mode: 0644]

index 3420387..522fc9a 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,7 @@
 .PHONY: test target clean
 
-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
+CF_OBJS:=cf.o mobius.o famous.o bihom.o taylor.o test.o newton.o tee.o
+TESTS:=bihom_test cf_test famous_test mobius_test newton_test tee_test
 BINS:=pi epow
 
 target : $(BINS)
diff --git a/cf.c b/cf.c
index 780f2c5..d6515e1 100644 (file)
--- a/cf.c
+++ b/cf.c
@@ -126,6 +126,13 @@ void cf_put_int(cf_t cf, int n) {
   mpz_clear(z);
 }
 
+void cf_signal(cf_t cf) {
+  sem_post(&cf->demand_sem);
+}
+void cf_wait_special(cf_t cf) {
+  sem_wait(&cf->demand_sem);
+}
+
 void cf_get(mpz_t z, cf_t cf) {
   pthread_mutex_lock(&cf->chan_mu);
   if (!cf->chan) {
diff --git a/cf.h b/cf.h
index aab54a4..67c7f27 100644 (file)
--- a/cf.h
+++ b/cf.h
@@ -25,8 +25,15 @@ int cf_wait(cf_t cf);
 
 void *cf_data(cf_t cf);
 
-// cf_mobius.c:
+void cf_signal(cf_t cf); // For tee.
+void cf_waitspecial(cf_t cf);
 
+// From cf_tee.c:
+//
+void cf_tee(cf_t *out_array, cf_t in);
+
+// From 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_cf_convergent(cf_t x);
diff --git a/newton.c b/newton.c
new file mode 100644 (file)
index 0000000..aeb8ce8
--- /dev/null
+++ b/newton.c
@@ -0,0 +1,213 @@
+// Solve quadratic equations involving continued fractions.
+// A slight misnomer: Gosper describes how to use Newton's method, but
+// I use binary search instead, and assume unique roots in the right ranges.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <gmp.h>
+#include "cf.h"
+
+// Represents:
+//
+// a0 xy + a1 x + a2 y + a3
+// ------------------------
+// a4 xy - a0 x + a5 y - a2
+
+struct newton_data_s {
+  cf_t x;
+  mpz_t a[6];
+  mpz_t lower;
+};
+typedef struct newton_data_s newton_data_t[1];
+typedef struct newton_data_s *newton_data_ptr;
+
+// In the 3D table, we have two sets of equations. Left is a0, b0, c0,
+// Right is a1, b1, c1, and these are the terms involving x:
+//
+//   a0 y + b0     a1 y + b1
+//   ---------     ---------
+//   c0 y - a0     c1 y - a1
+//
+// I work left to right, and top to bottom, unlike Gosper.
+//
+// This notation is different to that in cf_bihom.c, and similar to cf_mobius.
+
+struct abc_s {
+  mpz_t a0, a1;
+  mpz_t b0, b1;
+  mpz_t c0, c1;
+};
+typedef struct abc_s abc_t[1];
+
+void abc_init(abc_t p) {
+  mpz_init(p->a0); mpz_init(p->a1);
+  mpz_init(p->b0); mpz_init(p->b1);
+  mpz_init(p->c0); mpz_init(p->c1);
+}
+
+void abc_clear(abc_t p) {
+  mpz_clear(p->a0); mpz_clear(p->a1);
+  mpz_clear(p->b0); mpz_clear(p->b1);
+  mpz_clear(p->c0); mpz_clear(p->c1);
+}
+
+void abc_set_coeff(abc_t p, mpz_t a[6]) {
+  mpz_set(p->a0, a[2]); mpz_set(p->b0, a[3]); mpz_set(p->c0, a[5]);
+  mpz_set(p->a1, a[0]); mpz_set(p->b1, a[1]); mpz_set(p->c1, a[4]);
+}
+
+void abc_print(abc_t p) {
+  gmp_printf("%Zd y + %Zd    %Zd y + %Zd\n", p->a0, p->b0, p->a1, p->b1);
+  gmp_printf("%Zd y - %Zd    %Zd y - %Zd\n", p->c0, p->a0, p->c1, p->a1);
+  gmp_printf("(%Zd+sqrt(%Zd^2+%Zd*%Zd))/%Zd\n", p->a0, p->a0, p->b0, p->c0, p->c0);
+  gmp_printf("%Zd*z^2-2*%Zd*z-%Zd\n", p->c0, p->a0, p->b0);
+}
+
+// Finds smallest root greater than given lower bound.
+// Assumes there exists exactly one root with this property.
+// (Impossible to solve equations if roots have same integer part at
+// the moment. Don't do that.)
+// TODO: Rewrite so you choose if you want smaller or greater root.
+// or so it returns both solutions in an array.
+static void *newton(cf_t cf) {
+  newton_data_ptr nd = cf_data(cf);
+  abc_t p;
+  abc_init(p);
+  abc_set_coeff(p, nd->a);
+  cf_t x = nd->x;
+  mpz_t z, z0, z1, one, pow2, t0, t1, t2;
+  mpz_init(z); mpz_init(z0); mpz_init(z1); mpz_init(pow2);
+  mpz_init(t0); mpz_init(t1); mpz_init(t2);
+  mpz_init(one);
+  mpz_set_ui(one, 1);
+
+  void move_right() {
+    cf_get(z, x);
+    mpz_mul(t0, z, p->b1);
+    mpz_add(t0, t0, p->b0);
+    mpz_set(p->b0, p->b1);
+    mpz_set(p->b1, t0);
+
+    mpz_mul(t0, z, p->a1);  mpz_mul(t1, z, p->c1);
+    mpz_add(t0, t0, p->a0); mpz_add(t1, t1, p->c0);
+    mpz_set(p->a0, p->a1);  mpz_set(p->c0, p->c1);
+    mpz_set(p->a1, t0);     mpz_set(p->c1, t1);
+  }
+
+  void move_down() {
+    // Recurrence relation with z, then
+    // Subtract and invert z.
+    mpz_mul(t0, z, p->a0);
+    mpz_add(t0, t0, p->b0);
+    mpz_mul(t1, z, p->c0);
+    mpz_sub(t1, t1, p->a0);
+    mpz_mul(p->a0, z, t1);
+    mpz_set(p->b0, p->c0);
+    mpz_sub(p->c0, t0, p->a0);
+    mpz_set(p->a0, t1);
+
+    mpz_mul(t0, z, p->a1);
+    mpz_add(t0, t0, p->b1);
+    mpz_mul(t1, z, p->c1);
+    mpz_sub(t1, t1, p->a1);
+    mpz_mul(p->a1, z, t1);
+    mpz_set(p->b1, p->c1);
+    mpz_sub(p->c1, t0, p->a1);
+    mpz_set(p->a1, t1);
+  }
+
+  int sign_quad() {
+    // Returns sign of c0 z^2 - 2 a0 z - b0
+    mpz_mul_si(t0, p->a0, -2);
+    mpz_addmul(t0, p->c0, z);
+    mpz_mul(t0, t0, z);
+    mpz_sub(t0, t0, p->b0);
+    return mpz_sgn(t0);
+  }
+
+  int sign_quad1() {
+    // Returns sign of c1 z^2 - 2 a1 z - b1
+    mpz_mul_si(t0, p->a1, -2);
+    mpz_addmul(t0, p->c1, z);
+    mpz_mul(t0, t0, z);
+    mpz_sub(t0, t0, p->b1);
+    return mpz_sgn(t0);
+  }
+
+  move_right();  // Get rid of pathological cases.
+  move_right();
+  move_right();
+
+  // Get integer part, starting search from given lower bound.
+  void binary_search(mpz_ptr lower) {
+    while (!mpz_sgn(p->c0)) move_right();
+    for (;;) {
+      mpz_set(z0, lower);
+      mpz_set(z, lower);
+      int sign = sign_quad();
+      mpz_set_ui(pow2, 1);
+      for (;;) {
+       mpz_add(z, z0, pow2);
+       if (sign_quad() != sign) break;
+       mpz_mul_2exp(pow2, pow2, 1);
+      }
+      mpz_set(z1, z);
+
+      for (;;) {
+       mpz_add(z, z0, z1);
+       mpz_div_2exp(z, z, 1);
+       if (!mpz_cmp(z, z0)) break;
+       if (sign_quad() == sign) {
+         mpz_set(z0, z);
+       } else {
+         mpz_set(z1, z);
+       }
+      }
+      sign = sign_quad1();
+      mpz_set(z, z1);
+      if (sign_quad1() != sign) break;
+      move_right();
+    }
+  }
+
+  binary_search(nd->lower);
+  cf_put(cf, z0);
+  mpz_set(z, z0);
+  move_down();
+
+  while(cf_wait(cf)) {
+    binary_search(one);
+    cf_put(cf, z0);
+    mpz_set(z, z0);
+    move_down();
+  }
+  abc_clear(p);
+  mpz_clear(z); mpz_clear(z0); mpz_clear(z1); mpz_clear(pow2);
+  mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
+  for (int i = 0; i < 6; i++) mpz_clear(nd->a[i]);
+  free(nd);
+  return NULL;
+}
+
+cf_t cf_new_newton(cf_t x, mpz_t a[6], mpz_t lower) {
+  newton_data_ptr p = malloc(sizeof(*p));
+  p->x = x;
+  for (int i = 0; i < 6; i++) {
+    mpz_init(p->a[i]);
+    mpz_set(p->a[i], a[i]);
+  }
+  mpz_init(p->lower);
+  mpz_set(p->lower, lower);
+  return cf_new(newton, p);
+}
+
+cf_t cf_new_sqrt(cf_t x) {
+  newton_data_ptr p = malloc(sizeof(*p));
+  p->x = x;
+  mpz_init(p->lower);
+  for (int i = 0; i < 6; i++) mpz_init(p->a[i]);
+  // Solve y = x/y.
+  mpz_set_si(p->a[1], 1);
+  mpz_set_si(p->a[5], 1);
+  return cf_new(newton, p);
+}
diff --git a/newton_test.c b/newton_test.c
new file mode 100644 (file)
index 0000000..1e0655b
--- /dev/null
@@ -0,0 +1,56 @@
+// Test Gosper's example: find coth 1/2 via Newton's method.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+#include "test.h"
+
+static void *cot1_fn(cf_t cf) {
+  int n = 1;
+  while(cf_wait(cf)) {
+    cf_put_int(cf, n);
+    n += 2;
+  }
+  return NULL;
+}
+
+int main() {
+  mpz_t a[6];
+  int i;
+  for (i = 0; i < 6; i++) mpz_init(a[i]);
+  cf_t x, y;
+
+  // Gosper's example.
+  mpz_set_si(a[0], 1);
+  mpz_set_si(a[3], -1);
+  mpz_set_si(a[5], 1);
+  x = cf_new_const(cot1_fn);
+  y = cf_new_newton(x, a, a[0]);
+  CF_EXPECT_DEC(y, "2.16395341373865284877");
+  cf_free(x);
+  cf_free(y);
+  
+  // Confirm sqrt(1-(sin 1)^2) = cos 1
+  x = cf_new_sin1();
+  y = cf_new_sin1();  // TODO: Use tee here.
+  mpz_t b[8];
+  for (i = 0; i < 8; i++) mpz_init(b[i]);
+  mpz_set_si(b[0], -1);
+  mpz_set_si(b[3], 1);
+  mpz_set_si(b[7], 1);
+
+  cf_t bi = cf_new_bihom(x, y, b);
+  cf_t n = cf_new_sqrt(bi);
+
+  CF_EXPECT_DEC(n, "0.54030230586813971740");
+  cf_free(x);
+  cf_free(y);
+  cf_free(bi);
+  cf_free(n);
+
+  for (i = 0; i < 8; i++) mpz_clear(b[i]);
+  for (i = 0; i < 6; i++) mpz_clear(a[i]);
+  return 0;
+}
diff --git a/tee.c b/tee.c
new file mode 100644 (file)
index 0000000..82d2985
--- /dev/null
+++ b/tee.c
@@ -0,0 +1,118 @@
+// Tee: read input channel and write to two output channels
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+
+struct tee_s {
+  int n;
+  cf_t parent;
+};
+typedef struct tee_s tee_t[1];
+typedef struct tee_s *tee_ptr;
+
+static void *branch_loop(cf_t cf) {
+  tee_ptr t = cf_data(cf);
+  mpz_t z;
+  mpz_init(z);
+  while(cf_wait(cf)) {
+    cf_put_int(t->parent, t->n);  // Causes parent to put to cf.
+    cf_signal(t->parent);
+  }
+  cf_put_int(t->parent, -1 - t->n);  // Notify parent of our destruction.
+  mpz_clear(z);
+  free(t);
+  return NULL;
+}
+
+struct parent_data_s {
+  cf_t kid[2];
+  cf_t in;
+};
+typedef struct parent_data_s parent_data_t[1];
+typedef struct parent_data_s *parent_data_ptr;
+
+// TODO: Destroy after children are destroyed.
+static void *parent_loop(cf_t cf) {
+  struct backlog_s {
+    mpz_t z;
+    struct backlog_s *next;
+  };
+  typedef struct backlog_s *backlog_ptr;
+
+  backlog_ptr head[2];
+  backlog_ptr last[2];
+  head[0] = NULL;
+  head[1] = NULL;
+  last[0] = NULL;
+  last[1] = NULL;
+  mpz_t z;
+  mpz_init(z);
+  parent_data_ptr pd = cf_data(cf);
+  for (;;) {
+    cf_wait_special(cf);
+    cf_get(z, cf);
+    int k = mpz_get_ui(z);
+    if (k < 0) {
+      k = -k + 1;
+      pd->kid[k] = NULL;
+      backlog_ptr pnext = head[k], p;
+      do {
+       p = pnext;
+       pnext = p->next;
+       mpz_clear(p->z);
+       free(p);
+      } while(pnext);
+      if (!pd->kid[1-k]) break;
+    } else {
+      if (head[k]) {
+       backlog_ptr p = head[k];
+       mpz_set(z, p->z);
+       mpz_clear(p->z);
+       head[k] = p->next;
+       free(p);
+       if (!head[k]) last[k] = NULL;
+      } else {
+       cf_get(z, pd->in);
+       if (pd->kid[1-k]) {
+         backlog_ptr p = malloc(sizeof(*p));
+         mpz_init(p->z);
+         mpz_set(p->z, z);
+         p->next = NULL;
+         if (last[1-k]) last[1-k]->next = p;
+         last[1-k] = p;
+         if (!head[1-k]) head[1-k] = p;
+       }
+      }
+    }
+    cf_put(pd->kid[k], z);
+  }
+  mpz_clear(z);
+  return NULL;
+}
+
+cf_t cf_new_parent(cf_t in) {
+  parent_data_ptr p = malloc(sizeof(*p));
+  p->in = in;
+  p->kid[0] = NULL;
+  p->kid[1] = NULL;
+  return cf_new(parent_loop, p);
+}
+
+cf_t cf_new_branch(cf_t parent, int n) {
+  tee_ptr t = malloc(sizeof(*t));
+  t->n = n;
+  t->parent = parent;
+  parent_data_ptr p = cf_data(parent);
+  cf_t res = cf_new(branch_loop, t);
+  p->kid[n] = res;
+  return res;
+}
+
+void cf_tee(cf_t *out_array, cf_t in) {
+  cf_t parent = cf_new_parent(in);
+  out_array[0] = cf_new_branch(parent, 0);
+  out_array[1] = cf_new_branch(parent, 1);
+}
diff --git a/tee_test.c b/tee_test.c
new file mode 100644 (file)
index 0000000..1ac2e46
--- /dev/null
@@ -0,0 +1,45 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+#include "test.h"
+
+static void *count_int_fn(cf_t cf) {
+  int n = 1;
+  while(cf_wait(cf)) {
+    cf_put_int(cf, n);
+    n++;
+  }
+  return NULL;
+}
+
+int main() {
+  mpz_t z;
+  mpz_init(z);
+  cf_t x = cf_new_const(count_int_fn);
+  cf_t out[2];
+  int i[2];
+  cf_tee(out, x);
+
+  void get(int k, int n) {
+    while(n) {
+      cf_get(z, out[k]);
+      EXPECT(!mpz_cmp_ui(z, i[k]++));
+      n--;
+    }
+  }
+
+  i[0] = 1;
+  i[1] = 1;
+  get(0, 5);
+  get(1, 10);
+  get(0, 20);
+
+  cf_free(out[0]);
+  get(1, 100);
+  cf_free(out[1]);
+
+  mpz_clear(z);
+  return 0;
+}