taylor.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <gmp.h>
   4 #include "cf.h"
   5
   6 // sin 1 from Taylor series
   7 static void *sin1_expansion(cf_t cf) {
   8   mpz_t p0, q0;  // p0/q0 = last convergent
   9   mpz_t p1, q1;  // p1/q1 = convergent
  10   mpz_t r0, s0;  // r0/s0 = lower bound for sin 1
  11   mpz_t r1, s1;  // r1/s1 = upper bound for sin 1
  12   mpq_t lim0, lim1;  // Related to error of lower and upper bounds:
  13                      //   error = 1 / (lim)
  14   mpz_t t0, t1, t2;  // Temporary variables.
  15   mpq_t tq0, tq1;
  16
  17   mpz_init(r0); mpz_init(r1); mpz_init(s0); mpz_init(s1);
  18   mpz_init(p0); mpz_init(p1); mpz_init(q0); mpz_init(q1);
  19   mpz_init(t0); mpz_init(t1); mpz_init(t2);
  20   mpq_init(lim0); mpq_init(lim1);
  21   mpq_init(tq0); mpq_init(tq1);
  22
  23   // Zero causes all sorts of problems.
  24   // Output it and forget about it.
  25   cf_put_int(cf, 0);
  26
  27   // Initialize convergents.
  28   mpz_set_ui(p0, 1);
  29   mpz_set_ui(q1, 1);
  30
  31   // sin 1 = 1 - 1/6 + ...
  32   mpz_set_ui(r0, 5);
  33   mpz_set_ui(s0, 6);
  34   mpz_set_ui(r1, 1);
  35   mpz_set_ui(s1, 1);
  36   mpq_set_ui(lim0, 1, 120);
  37   mpq_set_ui(lim1, 1, 6);
  38
  39   int k = 1; // Number of Taylor series terms we've computed,
  40              // numbered from 0, unlike continued fraction terms.
  41   int i = 2; // Number of convergent we're computing.
  42
  43   void increase_limit() {
  44     k++;
  45     mpz_mul_ui(r1, r0, 2 * k);
  46     mpz_mul_ui(r1, r1, 2 * k + 1);
  47     mpz_mul_ui(s1, s0, 2 * k);
  48     mpz_mul_ui(s1, s1, 2 * k + 1);
  49     mpz_add_ui(r1, r1, 1);
  50     k++;
  51     mpz_mul_ui(r0, r1, 2 * k);
  52     mpz_mul_ui(r0, r0, 2 * k + 1);
  53     mpz_mul_ui(s0, s1, 2 * k);
  54     mpz_mul_ui(s0, s0, 2 * k + 1);
  55     mpz_sub_ui(r0, r0, 1);
  56     mpz_set(mpq_denref(lim1), s0);
  57     mpz_mul_ui(mpq_denref(lim0), mpq_denref(lim1), 2 * k + 2);
  58     mpz_mul_ui(mpq_denref(lim0), mpq_denref(lim0), 2 * k + 3);
  59     gmp_printf("lower: %Zd / %Zd, %Qd\n", r0, s0, lim0);
  60     gmp_printf("upper: %Zd / %Zd, %Qd\n", r1, s1, lim1);
  61   }
  62   while (cf_wait(cf)) {
  63     for(;;) {
  64       printf("find x\n");
  65       // Find largest x such that (p1 x + p0)/(q1 x + q0)
  66       // <= r0/s0 for odd i, and >= r1/s1 for even i
  67       //  =>  x < (q0 r - p0 s)/(p1 s - q1 r)
  68       // for appropriate choice of r, s.
  69       mpz_ptr r, s;
  70       mpq_ptr lim;
  71       if (i & 1) {
  72         r = r0; s = s0; lim = lim0;
  73       } else {
  74         r = r1; s = s1; lim = lim1;
  75       }
  76       mpz_mul(t0, q0, r);
  77       mpz_mul(t1, p0, s);
  78       mpz_sub(t2, t0, t1);
  79       mpz_mul(t0, p1, s);
  80       mpz_mul(t1, q1, r);
  81       mpz_sub(t1, t0, t1);
  82       mpz_fdiv_q(t0, t2, t1);
  83       gmp_printf("soln: %Zd = %Zd/%Zd\n", t0, t2, t1);
  84
  85       // Then x is the next convergent provided it is within limits.
  86       if (!mpz_sgn(t0)) {
  87         increase_limit();
  88       } else {
  89         // Check if t1/t2 is within the convergent bound
  90         mpz_mul(t1, p1, t0);
  91         mpz_add(t1, t1, p0);
  92         mpz_mul(t2, q1, t0);
  93         mpz_add(t2, t2, q0);
  94
  95         mpq_set_num(tq0, t1);
  96         mpq_set_den(tq0, t2);
  97         mpq_set_num(tq1, r);
  98         mpq_set_den(tq1, s);
  99         gmp_printf("check: %Qd - %Qd\n", tq0, tq1);
 100         mpq_sub(tq0, tq0, tq1);
 101         if (mpq_sgn(tq0) < 0) {
 102           mpq_neg(tq0, tq0);
 103         }
 104         // tq1 = 1 / 2 t2^2
 105         mpz_set_ui(mpq_numref(tq1), 1);
 106         mpz_mul(mpq_denref(tq1), t2, t2);
 107         mpz_mul_ui(mpq_denref(tq1), mpq_denref(tq1), 2);
 108
 109         gmp_printf("check: %Qd + %Qd\n", tq0, lim);
 110         mpq_add(tq0, tq0, lim);
 111         gmp_printf("%Qd > %Qd?\n", tq0, tq1);
 112         if (mpq_cmp(tq0, tq1) >= 0) {
 113           increase_limit();
 114         } else break;
 115       }
 116     };
 117     gmp_printf("term: %Zd\n", t0);
 118     cf_put(cf, t0);
 119     i++;
 120     mpz_set(p0, p1);
 121     mpz_set(q0, q1);
 122     mpz_set(p1, t1);
 123     mpz_set(q1, t2);
 124   }
 125   return NULL;
 126 }
 127
 128 int main(int argc, char *argv[]) {
 129   cf_t x = cf_new_const(sin1_expansion);
 130   cf_t conv = cf_new_convergent(x);
 131   mpz_t z;
 132   mpz_init(z);
 133   for (int i = 0; i < 20; i++) {
 134     cf_get(z, conv); gmp_printf("%d: %Zd/", i, z);
 135     cf_get(z, conv); gmp_printf("%Zd\n", z);
 136   }
 137   return 0;
 138 }