From 9dad1e95571c9122184c4fefd37c032a6e101d94 Mon Sep 17 00:00:00 2001 From: Ben Lynn Date: Tue, 16 Sep 2008 04:20:08 -0700 Subject: [PATCH] Improved Makefile, tests. --- Makefile | 35 +++++++----- cf.h | 16 ++++-- cf_test.c | 2 +- famous.c | 13 +++++ mobius.c | 169 ++++++++++++++++++++++++++++++++++++++++------------------ mobius_test.c | 74 +++++++++++++++++++++++++ sqrt2.c | 54 ------------------- 7 files changed, 239 insertions(+), 124 deletions(-) rewrite Makefile (99%) create mode 100644 mobius_test.c delete mode 100644 sqrt2.c diff --git a/Makefile b/Makefile dissimilarity index 99% index 5f61db1..bc49a6c 100644 --- 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 --- 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(); diff --git a/cf_test.c b/cf_test.c index c83694b..c417a3e 100644 --- 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)); diff --git a/famous.c b/famous.c index 9e6e06b..48e279f 100644 --- a/famous.c +++ b/famous.c @@ -3,6 +3,19 @@ #include #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; diff --git a/mobius.c b/mobius.c index 8507449..0d00f12 100644 --- 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 index 0000000..150d840 --- /dev/null +++ b/mobius_test.c @@ -0,0 +1,74 @@ +#include +#include +#include +#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 index 6109ff4..0000000 --- a/sqrt2.c +++ /dev/null @@ -1,54 +0,0 @@ -#include -#include -#include -#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; -} -- 2.11.4.GIT