From 5c7432317ea676f70642e3dec4d67f28168d7f1c Mon Sep 17 00:00:00 2001 From: Ben Lynn Date: Wed, 17 Sep 2008 04:36:10 -0700 Subject: [PATCH] Implemented tee. Also realized I had forgotten to add newton* files. --- Makefile | 4 +- cf.c | 7 ++ cf.h | 9 ++- newton.c | 213 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ newton_test.c | 56 +++++++++++++++ tee.c | 118 ++++++++++++++++++++++++++++++++ tee_test.c | 45 +++++++++++++ 7 files changed, 449 insertions(+), 3 deletions(-) create mode 100644 newton.c create mode 100644 newton_test.c create mode 100644 tee.c create mode 100644 tee_test.c diff --git a/Makefile b/Makefile index 3420387..522fc9a 100644 --- 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 --- 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 --- 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 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 +#include +#include +#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 index 0000000..1e0655b --- /dev/null +++ b/newton_test.c @@ -0,0 +1,56 @@ +// Test Gosper's example: find coth 1/2 via Newton's method. + +#include +#include +#include +#include +#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 index 0000000..82d2985 --- /dev/null +++ b/tee.c @@ -0,0 +1,118 @@ +// Tee: read input channel and write to two output channels + +#include +#include +#include +#include +#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 index 0000000..1ac2e46 --- /dev/null +++ b/tee_test.c @@ -0,0 +1,45 @@ +#include +#include +#include +#include +#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; +} -- 2.11.4.GIT