From c00f6567f32ddaa1f628f98546f9e1f0330ee90a Mon Sep 17 00:00:00 2001 From: Ben Lynn Date: Mon, 15 Sep 2008 04:34:01 -0700 Subject: [PATCH] Improved tests. Appears to be a bug in bihom.c. --- Makefile | 10 +++--- README | 16 +++++---- cf.h | 4 +++ famous_test.c | 78 ++++++++++------------------------------- taylor.c | 111 +++++++++++++++++++++++++++++++++++++++++++++++++++------- test.c | 35 ++++++++++++++++++ test.h | 11 ++++++ 7 files changed, 183 insertions(+), 82 deletions(-) rewrite famous_test.c (84%) create mode 100644 test.c diff --git a/Makefile b/Makefile index 16338a1..5f61db1 100644 --- 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 --- 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 --- 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__ diff --git a/famous_test.c b/famous_test.c dissimilarity index 84% index 5e7adb9..ec2ac79 100644 --- a/famous_test.c +++ b/famous_test.c @@ -1,59 +1,19 @@ -// Test famous continued fraction expansions are correct. - -#include -#include -#include -#include -#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 +#include +#include +#include +#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; +} diff --git a/taylor.c b/taylor.c index a834a48..838b793 100644 --- 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 index 0000000..075c714 --- /dev/null +++ b/test.c @@ -0,0 +1,35 @@ +#include +#include +#include +#include +#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 --- 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); -- 2.11.4.GIT