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   mpq_t r0;      // lower bound for sin 1
  11   mpq_t r1;      // upper bound for sin 1
  12   mpz_t t0, t1, t2;  // Temporary variables.
  13   mpq_t tq0;
  14
  15   mpq_init(r0); mpq_init(r1);
  16   mpz_init(p0); mpz_init(p1); mpz_init(q0); mpz_init(q1);
  17   mpz_init(t0); mpz_init(t1); mpz_init(t2);
  18   mpq_init(tq0);
  19
  20   // Zero causes all sorts of problems.
  21   // Output it and forget about it.
  22   cf_put_int(cf, 0);
  23
  24   // Initialize convergents.
  25   mpz_set_ui(p0, 1);
  26   mpz_set_ui(q1, 1);
  27
  28   // sin 1 = 1 - 1/6 + ...
  29   mpq_set_ui(r0, 5, 6);
  30   mpq_set_ui(r1, 1, 1);
  31
  32   int k = 1; // Number of Taylor series terms we've computed,
  33              // numbered from 0, unlike continued fraction terms.
  34   int i = 2; // Number of convergent we're computing.
  35
  36   void increase_limit() {
  37     k++;
  38     mpz_mul_ui(mpq_numref(r1), mpq_numref(r0), (2 * k) * (2 * k + 1));
  39     mpz_mul_ui(mpq_denref(r1), mpq_denref(r0), (2 * k) * (2 * k + 1));
  40     mpz_add_ui(mpq_numref(r1), mpq_numref(r1), 1);
  41     k++;
  42     mpz_mul_ui(mpq_numref(r0), mpq_numref(r1), (2 * k) * (2 * k + 1));
  43     mpz_mul_ui(mpq_denref(r0), mpq_denref(r1), (2 * k) * (2 * k + 1));
  44     mpz_sub_ui(mpq_numref(r0), mpq_numref(r0), 1);
  45   }
  46   while (cf_wait(cf)) {
  47     for(;;) {
  48       printf("find x\n");
  49       // Find largest x such that (p1 x + p0)/(q1 x + q0)
  50       // <= r0 for odd i, and >= r1 for even i
  51       //  =>  x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
  52       // for the appropriate r.
  53       mpq_ptr r = i & 1 ? r0 : r1;
  54       mpz_mul(t0, q0, mpq_numref(r));
  55       mpz_mul(t1, p0, mpq_denref(r));
  56       mpz_sub(t2, t0, t1);
  57       mpz_mul(t0, p1, mpq_denref(r));
  58       mpz_mul(t1, q1, mpq_numref(r));
  59       mpz_sub(t1, t0, t1);
  60       mpz_fdiv_q(t0, t2, t1);
  61
  62       // Then x is the next convergent provided substituting x + 1 overshoots
  63       // the other bound. (If not, we need better bounds to decide.)
  64       if (!mpz_sgn(t0)) {
  65         increase_limit();
  66       } else {
  67         mpz_add_ui(t1, t0, 1);
  68         mpz_set(mpq_numref(tq0), p0);
  69         mpz_addmul(mpq_numref(tq0), p1, t1);
  70         mpz_set(mpq_denref(tq0), q0);
  71         mpz_addmul(mpq_denref(tq0), q1, t1);
  72         if (i & 1) {
  73           if (mpq_cmp(tq0, r1) > 0) break;
  74           // Check (p1(x+1) + p0)/(q1(x+1) + q0) > r1
  75         } else {
  76           if (mpq_cmp(tq0, r0) < 0) break;
  77           // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
  78         }
  79         increase_limit();
  80       }
  81     };
  82     cf_put(cf, t0);
  83     i++;
  84     mpz_addmul(p0, p1, t0);
  85     mpz_addmul(q0, q1, t0);
  86     mpz_swap(p1, p0);
  87     mpz_swap(q1, q0);
  88   }
  89
  90   mpq_clear(r0); mpq_clear(r1);
  91   mpz_clear(p0); mpz_clear(p1); mpz_clear(q0); mpz_clear(q1);
  92   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
  93   mpq_clear(tq0);
  94   return NULL;
  95 }
  96
  97 int main(int argc, char *argv[]) {
  98   cf_t x = cf_new_const(sin1_expansion);
  99   cf_t conv = cf_new_convergent(x);
 100   mpz_t z;
 101   mpz_init(z);
 102   for (int i = 0; i < 20; i++) {
 103     cf_get(z, conv); gmp_printf("%d: %Zd/", i, z);
 104     cf_get(z, conv); gmp_printf("%Zd\n", z);
 105   }
 106   return 0;
 107 }