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   int k = 1;
   9
  10   mpz_t r, s;
  11   mpz_init(r); mpz_init(s);
  12
  13   mpz_set_ui(r, 1);
  14   mpz_set_ui(s, 1);
  15   int i;
  16   for (i = 1; i <= k ; i++ ) {
  17     mpz_mul_ui(s, s, (2 * i) * (2 * i + 1));
  18     mpz_mul_ui(r, r, (2 * i) * (2 * i + 1));
  19     if (i & 1) {
  20       mpz_sub_ui(r, r, 1);
  21     } else {
  22       mpz_add_ui(r, r, 1);
  23     }
  24   }
  25
  26   mpz_t limit;
  27   mpz_init(limit);
  28   // Compute limit = (2k + 1)! / 2
  29   mpz_mul_ui(limit, s, k);
  30   mpz_mul_ui(limit, s, 2 * k + 1);
  31
  32   mpz_t p0, p1;
  33   mpz_t q0, q1;
  34   mpz_init(p0); mpz_init(p1);
  35   mpz_init(q0); mpz_init(q1);
  36   mpz_set_ui(p1, 1);
  37   mpz_set_ui(q0, 1);
  38
  39   mpz_t a, b;
  40   mpz_init(a); mpz_init(b);
  41   mpz_t t0, t1;
  42   mpz_init(t0); mpz_init(t1);
  43
  44   mpz_set(a, r);
  45   mpz_set(b, s);
  46
  47   mpz_t t2;
  48   mpz_init(t2);
  49   // Pump out start of expansion.
  50   i = 0;
  51   for (;;) {
  52     mpz_fdiv_qr(t0, t1, a, b);
  53
  54     mpz_set(t2, q0);
  55     mpz_addmul(q0, q1, t0);
  56     mpz_mul(a, q0, q0);
  57     if (mpz_cmp(a, limit) >= 0) {
  58       mpz_set(q0, t2);
  59       break;
  60     }
  61     mpz_addmul(p0, p1, t0);
  62     mpz_swap(p0, p1);
  63     mpz_swap(q0, q1);
  64     i++;
  65     // gmp_printf("%Zd (%Zd) -> %Zd / %Zd\n", t0, t1, p1, q1);
  66     cf_put(cf, t0);
  67     // if (!mpz_sgn(t1)) die("sin(1) is rational!\n");
  68     mpz_set(a, b);
  69     mpz_set(b, t1);
  70   }
  71
  72   void increase_k() {
  73     k++;
  74     mpz_mul_ui(limit, limit, 2 * k);
  75     mpz_mul_ui(limit, limit, 2 * k + 1);
  76     mpz_mul_ui(s, s, (2 * k) * (2 * k + 1));
  77     mpz_mul_ui(r, r, (2 * k) * (2 * k + 1));
  78     if (k & 1) {
  79       mpz_sub_ui(r, r, 1);
  80     } else {
  81       mpz_add_ui(r, r, 1);
  82     }
  83 gmp_printf("r/s = %Zd/%Zd\n", r, s);
  84   }
  85   void increase_limit() {
  86     if (i & 1) {
  87       // Last output was for a lower bound.
  88       // We need an upper bound of sin 1.
  89       increase_k();
  90       if (k & 1) increase_k();
  91       printf("lower\n");
  92     } else {
  93       // Last output was for an upper bound.
  94       printf("upper\n");
  95       // We need a lower bound of sin 1.
  96       increase_k();
  97       if (!(k & 1)) increase_k();
  98     }
  99   }
 100
 101   increase_limit();
 102   while(cf_wait(cf)) {
 103     for(;;) {
 104       printf("find x\n");
 105       // Find largest x such that (p1 x + p0)/(q1 x + q0) R r/s
 106       // where R is < for even i, and > for odd i
 107       // i.e. (p1 s - q1 r)x R q0 r - p0 s
 108       mpz_mul(t0, q0, r);
 109       mpz_mul(t1, p0, s);
 110       mpz_sub(a, t0, t1);
 111       mpz_mul(t0, p1, s);
 112       mpz_mul(t1, q1, r);
 113       mpz_sub(b, t0, t1);
 114       mpz_fdiv_q(t0, a, b);
 115
 116       // Then x is the next convergent provided it is within limits.
 117       mpz_set(t2, q0);
 118       mpz_addmul(q0, q1, t0);
 119       mpz_mul(a, q0, q0);
 120       gmp_printf("%Zd > %Zd?\n", a, limit);
 121       if (mpz_cmp(a, limit) >= 0) {
 122         increase_limit();
 123         mpz_set(q0, t2);
 124       } else break;
 125     };
 126     gmp_printf("term: %Zd\n", t0);
 127     cf_put(cf, t0);
 128     i++;
 129     mpz_addmul(p0, p1, t0);
 130     mpz_swap(p0, p1);
 131     mpz_swap(q0, q1);
 132   }
 133
 134   mpz_clear(t2);
 135   return NULL;
 136 }
 137
 138 int main(int argc, char *argv[]) {
 139   cf_t x = cf_new_const(sin1_expansion);
 140   cf_t conv = cf_new_convergent(x);
 141   mpz_t z;
 142   mpz_init(z);
 143   for (int i = 0; i < 10; i++) {
 144     cf_get(z, conv); gmp_printf("%d: %Zd/", i, z);
 145     cf_get(z, conv); gmp_printf("%Zd\n", z);
 146   }
 147   return 0;
 148 }