From a7fe8343804a424f3eb0b171359cc497bd5c66dc Mon Sep 17 00:00:00 2001 From: Ben Lynn Date: Sun, 14 Sep 2008 22:43:50 -0700 Subject: [PATCH] Started Taylor series for sin. pi program takes command line argument: number of digits to print. --- README | 31 ++++++++++++++++++++++------- cf.h | 6 ++++++ epow.c | 2 +- famous.c | 3 +-- mobius.c | 10 ++++++++-- pi.c | 10 ++++++++-- sqrt2.c | 1 + taylor.c | 69 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 118 insertions(+), 14 deletions(-) create mode 100644 taylor.c diff --git a/README b/README index 6555e91..b9d6f93 100644 --- a/README +++ b/README @@ -1,5 +1,4 @@ = Continued Fraction Library = - Ben Lynn, September 2008 == Introduction == @@ -22,11 +21,29 @@ such an interesting phenomenon on something as mundane as a rational? With continued fractions, repeated terms occur only for square roots. Other infinite but patterned sequences represent other celebrated numbers, such as e and pi. -As with floating point numbers, we truncate infinite continued fractions to compute on them. But unlike floating point numbers, when more precision is desired -we can compute more terms without redoing any work. Imagine a real-time -application where on a timer interrupt, all continued fractions simply output -the last computed convergent. Under heavy system load, approximations are -coarser but the show goes on. +Let's say we use floating-point to compute some answer to 100 decimal places. +What if we now want 200 decimal places? Even if we possess the first 100 +digits, and carefully preserved the program state, we must start from scratch +and redo all our arithmetic with 200 digit precision. + +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 +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 +get to an exact number. Is it any wonder therefore that binary operations, +whose output must conform to this wondrous representation, are clumsy and +convoluted, and that deriving continued fraction expansions for many +reals is time-consuming? + +These drawbacks are almost certainly why continued fractions remain obscure. +Yet surely there are some problems that scream for continued fractions? +Imagine a real-time application where on a timer interrupt, all continued +fractions simply output the last computed convergent. Under heavy system load, +approximations are coarser but the show goes on. == Installation == @@ -88,7 +105,7 @@ study real numbers for its own sake. But perhaps then I'll find a perfect application for continued fractions, which will push me to fix deficiencies in this code. -- The API could use serious cosmetic surgery. +- The API could use cosmetic surgery. - I want to replace the condition variable with a semaphore to silence Helgrind (whose documentation states that condition variables almost unconditionally diff --git a/cf.h b/cf.h index 0721ec4..07f223a 100644 --- a/cf.h +++ b/cf.h @@ -50,6 +50,12 @@ cf_t cf_new_tan1(); cf_t cf_new_epow(mpz_t pow); cf_t cf_new_tanh(mpz_t z); +// This won't work because my code cannot handle negative denominators, +// and also assumes the sequence of convergents alternatively overshoot +// and undershoots the target. The tan expansion leads to a sequence of +// strictly increasing convergents (for positive input). +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]); diff --git a/epow.c b/epow.c index fbbacf9..4dde0cb 100644 --- a/epow.c +++ b/epow.c @@ -9,7 +9,7 @@ int main() { mpz_init(z); cf_t x, conv; - mpz_set_ui(z, 3); + mpz_set_si(z, 3); x = cf_new_tanh(z); printf("tanh 3 = "); diff --git a/famous.c b/famous.c index 89a5db1..9e6e06b 100644 --- a/famous.c +++ b/famous.c @@ -206,7 +206,6 @@ cf_t cf_new_tanh(mpz_t z) { return cf_new_one_arg_nonreg(gauss_tanh_expansion, z); } -// TODO: Handle negative convergents so this function works. cf_t cf_new_tan(mpz_t z) { - return cf_new_one_arg_nonreg(gauss_tanh_expansion, z); + return cf_new_one_arg_nonreg(gauss_tan_expansion, z); } diff --git a/mobius.c b/mobius.c index 57d196c..68a2525 100644 --- a/mobius.c +++ b/mobius.c @@ -1,5 +1,3 @@ -// TODO: Handle negative convergents - #include #include #include @@ -27,6 +25,13 @@ void pqset_clear(pqset_t pq) { mpz_clear(pq->q); } +void pqset_neg(pqset_t pq) { + mpz_neg(pq->p, pq->p); + mpz_neg(pq->q, pq->q); + mpz_neg(pq->pold, pq->pold); + mpz_neg(pq->qold, pq->qold); +} + void pqset_print(pqset_t pq) { gmp_printf("p's: %Zd %Zd\n", pq->pold, pq->p); gmp_printf("q's: %Zd %Zd\n", pq->qold, pq->q); @@ -178,6 +183,7 @@ static void *mobius_nonregular_throughput(cf_t cf) { int recur() { pqset_nonregular_recur(pq, num, denom); pqset_remove_gcd(pq, t0, t1); + if (mpz_sgn(pq->qold)) { mpz_fdiv_qr(t1, t0, pq->pold, pq->qold); mpz_mul(t2, t1, pq->q); diff --git a/pi.c b/pi.c index 4c5a302..bbb47ab 100644 --- a/pi.c +++ b/pi.c @@ -4,19 +4,25 @@ #include #include "cf.h" -int main() { +int main(int argc, char **argv) { mpz_t z; mpz_init(z); cf_t pi, conv; pi = cf_new_pi(); + int n = 1000; + if (argc > 1) { + n = atoi(argv[1]); + if (n <= 0) n = 100; + } conv = cf_new_cf_to_decimal(pi); - for (int i = 0; i <= 5000; i++) { + for (int i = 0; i <= n; i++) { cf_get(z, conv); gmp_printf("%Zd", z); if (!(i % 5)) putchar(' '); if (!(i % 50)) putchar('\n'); } + if (n % 50) putchar('\n'); cf_free(conv); cf_free(pi); mpz_clear(z); diff --git a/sqrt2.c b/sqrt2.c index 74db949..6762460 100644 --- a/sqrt2.c +++ b/sqrt2.c @@ -8,6 +8,7 @@ static void *sqrt2(cf_t cf) { while(cf_wait(cf)) { cf_put_int(cf, 2); } + return NULL; } cf_t cf_new_sqrt2() { diff --git a/taylor.c b/taylor.c new file mode 100644 index 0000000..b20b2b9 --- /dev/null +++ b/taylor.c @@ -0,0 +1,69 @@ +#include +#include +#include + +int main(int argc, char *argv[]) { + int k = 5; + if (argc > 1) k = atoi(argv[1]); + + mpz_t r, s; + mpz_init(r); mpz_init(s); + + mpz_set_ui(r, 1); + mpz_set_ui(s, 1); + int i; + for (i = 1; i <= k ; i++ ) { + mpz_mul_ui(s, s, (2 * i) * (2 * i + 1)); + mpz_mul_ui(r, r, (2 * i) * (2 * i + 1)); + if (i & 1) { + mpz_sub_ui(r, r, 1); + } else { + mpz_add_ui(r, r, 1); + } + } + gmp_printf("%Zd / %Zd\n", r, s); + + mpz_t limit; + mpz_init(limit); + // Compute limit = (2k + 1)! / 2 + mpz_mul_ui(limit, s, k); + mpz_mul_ui(limit, s, 2 * k + 1); + gmp_printf("%Zd\n", limit); + + mpz_t p0, p1; + mpz_t q0, q1; + mpz_init(p0); mpz_init(p1); + mpz_init(q0); mpz_init(q1); + mpz_set_ui(p1, 1); + mpz_set_ui(q0, 1); + + mpz_t a, b; + mpz_init(a); mpz_init(b); + mpz_t t0, t1; + mpz_init(t0); mpz_init(t1); + + mpz_set(a, r); + mpz_set(b, s); + for (;;) { + mpz_fdiv_qr(t0, t1, a, b); + + mpz_addmul(p0, p1, t0); + mpz_addmul(q0, q1, t0); + mpz_mul(a, q0, q0); + if (mpz_cmp(a, limit) >= 0) { + printf("limit exceeded\n"); break; + } + mpz_swap(p0, p1); + mpz_swap(q0, q1); + gmp_printf("%Zd (%Zd) -> %Zd / %Zd\n", t0, t1, p1, q1); + if (!mpz_sgn(t1)) break; + mpz_set(a, b); + mpz_set(b, t1); + } + + mpz_clear(limit); + mpz_clear(r); mpz_clear(s); + mpz_clear(t0); mpz_clear(t1); + mpz_clear(a); mpz_clear(b); + return 0; +} -- 2.11.4.GIT