From cd18efd1a07e04c1702f0f50ebe42b0feac133a2 Mon Sep 17 00:00:00 2001 From: Ben Lynn Date: Mon, 15 Sep 2008 03:18:26 -0700 Subject: [PATCH] Just realized my approach is flawed. --- taylor.c | 286 ++++++++++++++++++++++++++++++--------------------------------- 1 file changed, 138 insertions(+), 148 deletions(-) rewrite taylor.c (76%) diff --git a/taylor.c b/taylor.c dissimilarity index 76% index 21a20fa..0c51fe0 100644 --- a/taylor.c +++ b/taylor.c @@ -1,148 +1,138 @@ -#include -#include -#include -#include "cf.h" - -// sin 1 from Taylor series -static void *sin1_expansion(cf_t cf) { - int k = 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); - } - } - - 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); - - 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); - - mpz_t t2; - mpz_init(t2); - // Pump out start of expansion. - i = 0; - for (;;) { - mpz_fdiv_qr(t0, t1, a, b); - - mpz_set(t2, q0); - mpz_addmul(q0, q1, t0); - mpz_mul(a, q0, q0); - if (mpz_cmp(a, limit) >= 0) { - mpz_set(q0, t2); - break; - } - mpz_addmul(p0, p1, t0); - mpz_swap(p0, p1); - mpz_swap(q0, q1); - i++; - // gmp_printf("%Zd (%Zd) -> %Zd / %Zd\n", t0, t1, p1, q1); - cf_put(cf, t0); - // if (!mpz_sgn(t1)) die("sin(1) is rational!\n"); - mpz_set(a, b); - mpz_set(b, t1); - } - - void increase_k() { - k++; - mpz_mul_ui(limit, limit, 2 * k); - mpz_mul_ui(limit, limit, 2 * k + 1); - mpz_mul_ui(s, s, (2 * k) * (2 * k + 1)); - mpz_mul_ui(r, r, (2 * k) * (2 * k + 1)); - if (k & 1) { - mpz_sub_ui(r, r, 1); - } else { - mpz_add_ui(r, r, 1); - } -gmp_printf("r/s = %Zd/%Zd\n", r, s); - } - void increase_limit() { - if (i & 1) { - // Last output was for a lower bound. - // We need an upper bound of sin 1. - increase_k(); - if (k & 1) increase_k(); - printf("lower\n"); - } else { - // Last output was for an upper bound. - printf("upper\n"); - // We need a lower bound of sin 1. - increase_k(); - if (!(k & 1)) increase_k(); - } - } - - increase_limit(); - while(cf_wait(cf)) { - for(;;) { - printf("find x\n"); - // Find largest x such that (p1 x + p0)/(q1 x + q0) R r/s - // where R is < for even i, and > for odd i - // i.e. (p1 s - q1 r)x R q0 r - p0 s - mpz_mul(t0, q0, r); - mpz_mul(t1, p0, s); - mpz_sub(a, t0, t1); - mpz_mul(t0, p1, s); - mpz_mul(t1, q1, r); - mpz_sub(b, t0, t1); - mpz_fdiv_q(t0, a, b); - - // Then x is the next convergent provided it is within limits. - mpz_set(t2, q0); - mpz_addmul(q0, q1, t0); - mpz_mul(a, q0, q0); - gmp_printf("%Zd > %Zd?\n", a, limit); - if (mpz_cmp(a, limit) >= 0) { - increase_limit(); - mpz_set(q0, t2); - } else break; - }; - gmp_printf("term: %Zd\n", t0); - cf_put(cf, t0); - i++; - mpz_addmul(p0, p1, t0); - mpz_swap(p0, p1); - mpz_swap(q0, q1); - } - - mpz_clear(t2); - 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 < 10; i++) { - cf_get(z, conv); gmp_printf("%d: %Zd/", i, z); - cf_get(z, conv); gmp_printf("%Zd\n", z); - } - return 0; -} +#include +#include +#include +#include "cf.h" + +// sin 1 from Taylor series +static void *sin1_expansion(cf_t cf) { + mpz_t p0, q0; // p0/q0 = last convergent + mpz_t p1, q1; // p1/q1 = convergent + mpz_t r0, s0; // r0/s0 = lower bound for sin 1 + mpz_t r1, s1; // r1/s1 = upper bound for sin 1 + mpq_t lim0, lim1; // Related to error of lower and upper bounds: + // error = 1 / (lim) + mpz_t t0, t1, t2; // Temporary variables. + mpq_t tq0, tq1; + + mpz_init(r0); mpz_init(r1); mpz_init(s0); mpz_init(s1); + mpz_init(p0); mpz_init(p1); mpz_init(q0); mpz_init(q1); + mpz_init(t0); mpz_init(t1); mpz_init(t2); + mpq_init(lim0); mpq_init(lim1); + mpq_init(tq0); mpq_init(tq1); + + // 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); + + // sin 1 = 1 - 1/6 + ... + mpz_set_ui(r0, 5); + mpz_set_ui(s0, 6); + mpz_set_ui(r1, 1); + mpz_set_ui(s1, 1); + mpq_set_ui(lim0, 1, 120); + mpq_set_ui(lim1, 1, 6); + + 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 increase_limit() { + k++; + mpz_mul_ui(r1, r0, 2 * k); + mpz_mul_ui(r1, r1, 2 * k + 1); + mpz_mul_ui(s1, s0, 2 * k); + mpz_mul_ui(s1, s1, 2 * k + 1); + mpz_add_ui(r1, r1, 1); + k++; + mpz_mul_ui(r0, r1, 2 * k); + mpz_mul_ui(r0, r0, 2 * k + 1); + mpz_mul_ui(s0, s1, 2 * k); + mpz_mul_ui(s0, s0, 2 * k + 1); + mpz_sub_ui(r0, r0, 1); + mpz_set(mpq_denref(lim1), s0); + mpz_mul_ui(mpq_denref(lim0), mpq_denref(lim1), 2 * k + 2); + mpz_mul_ui(mpq_denref(lim0), mpq_denref(lim0), 2 * k + 3); + gmp_printf("lower: %Zd / %Zd, %Qd\n", r0, s0, lim0); + gmp_printf("upper: %Zd / %Zd, %Qd\n", r1, s1, lim1); + } + while (cf_wait(cf)) { + for(;;) { + printf("find x\n"); + // Find largest x such that (p1 x + p0)/(q1 x + q0) + // <= r0/s0 for odd i, and >= r1/s1 for even i + // => x < (q0 r - p0 s)/(p1 s - q1 r) + // for appropriate choice of r, s. + mpz_ptr r, s; + mpq_ptr lim; + if (i & 1) { + r = r0; s = s0; lim = lim0; + } else { + r = r1; s = s1; lim = lim1; + } + mpz_mul(t0, q0, r); + mpz_mul(t1, p0, s); + mpz_sub(t2, t0, t1); + mpz_mul(t0, p1, s); + mpz_mul(t1, q1, r); + mpz_sub(t1, t0, t1); + mpz_fdiv_q(t0, t2, t1); + gmp_printf("soln: %Zd = %Zd/%Zd\n", t0, t2, t1); + + // Then x is the next convergent provided it is within limits. + if (!mpz_sgn(t0)) { + increase_limit(); + } else { + // Check if t1/t2 is within the convergent bound + mpz_mul(t1, p1, t0); + mpz_add(t1, t1, p0); + mpz_mul(t2, q1, t0); + mpz_add(t2, t2, q0); + + mpq_set_num(tq0, t1); + mpq_set_den(tq0, t2); + mpq_set_num(tq1, r); + mpq_set_den(tq1, s); + gmp_printf("check: %Qd - %Qd\n", tq0, tq1); + mpq_sub(tq0, tq0, tq1); + if (mpq_sgn(tq0) < 0) { + mpq_neg(tq0, tq0); + } + // tq1 = 1 / 2 t2^2 + mpz_set_ui(mpq_numref(tq1), 1); + mpz_mul(mpq_denref(tq1), t2, t2); + mpz_mul_ui(mpq_denref(tq1), mpq_denref(tq1), 2); + + gmp_printf("check: %Qd + %Qd\n", tq0, lim); + mpq_add(tq0, tq0, lim); + gmp_printf("%Qd > %Qd?\n", tq0, tq1); + if (mpq_cmp(tq0, tq1) >= 0) { + increase_limit(); + } else break; + } + }; + gmp_printf("term: %Zd\n", t0); + cf_put(cf, t0); + i++; + mpz_set(p0, p1); + mpz_set(q0, q1); + mpz_set(p1, t1); + mpz_set(q1, t2); + } + 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); + } + return 0; +} -- 2.11.4.GIT