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 more_taylor() {
  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       // Find largest x such that (p1 x + p0)/(q1 x + q0)
  49       // <= r0 for odd i, and >= r1 for even i
  50       //  =>  x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
  51       // for the appropriate r.
  52       mpq_ptr r = i & 1 ? r0 : r1;
  53       mpz_mul(t0, q0, mpq_numref(r));
  54       mpz_mul(t1, p0, mpq_denref(r));
  55       mpz_sub(t2, t0, t1);
  56       mpz_mul(t0, p1, mpq_denref(r));
  57       mpz_mul(t1, q1, mpq_numref(r));
  58       mpz_sub(t1, t0, t1);
  59       mpz_fdiv_q(t0, t2, t1);
  60
  61       // Then x is the next convergent provided substituting x + 1 overshoots
  62       // the other bound. (If not, we need better bounds to decide.)
  63       if (!mpz_sgn(t0)) {
  64         more_taylor();
  65       } else {
  66         mpz_add_ui(t1, t0, 1);
  67         mpz_set(mpq_numref(tq0), p0);
  68         mpz_addmul(mpq_numref(tq0), p1, t1);
  69         mpz_set(mpq_denref(tq0), q0);
  70         mpz_addmul(mpq_denref(tq0), q1, t1);
  71         if (i & 1) {
  72           if (mpq_cmp(tq0, r1) > 0) break;
  73           // Check (p1(x+1) + p0)/(q1(x+1) + q0) > r1
  74         } else {
  75           if (mpq_cmp(tq0, r0) < 0) break;
  76           // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
  77         }
  78         more_taylor();
  79       }
  80     };
  81     cf_put(cf, t0);
  82     i++;
  83     mpz_addmul(p0, p1, t0);
  84     mpz_addmul(q0, q1, t0);
  85     mpz_swap(p1, p0);
  86     mpz_swap(q1, q0);
  87   }
  88
  89   mpq_clear(r0); mpq_clear(r1);
  90   mpz_clear(p0); mpz_clear(p1); mpz_clear(q0); mpz_clear(q1);
  91   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
  92   mpq_clear(tq0);
  93   return NULL;
  94 }
  95
  96 cf_t cf_new_sin1() {
  97   return cf_new_const(sin1_expansion);
  98 }
  99
 100 // cos 1 from Taylor series
 101 static void *cos1_expansion(cf_t cf) {
 102   mpz_t p0, q0;  // p0/q0 = last convergent
 103   mpz_t p1, q1;  // p1/q1 = convergent
 104   mpq_t r0;      // lower bound for sin 1
 105   mpq_t r1;      // upper bound for sin 1
 106   mpz_t t0, t1, t2;  // Temporary variables.
 107   mpq_t tq0;
 108
 109   mpq_init(r0); mpq_init(r1);
 110   mpz_init(p0); mpz_init(p1); mpz_init(q0); mpz_init(q1);
 111   mpz_init(t0); mpz_init(t1); mpz_init(t2);
 112   mpq_init(tq0);
 113
 114   // Zero causes all sorts of problems.
 115   // Output it and forget about it.
 116   cf_put_int(cf, 0);
 117
 118   // Initialize convergents.
 119   mpz_set_ui(p0, 1);
 120   mpz_set_ui(q1, 1);
 121
 122   // cos 1 = 1 - 1/2 + ...
 123   mpq_set_ui(r0, 1, 2);
 124   mpq_set_ui(r1, 1, 1);
 125
 126   int k = 1; // Number of Taylor series terms we've computed,
 127              // numbered from 0, unlike continued fraction terms.
 128   int i = 2; // Number of convergent we're computing.
 129
 130   void more_taylor() {
 131     k++;
 132     mpz_mul_ui(mpq_numref(r1), mpq_numref(r0), (2 * k) * (2 * k - 1));
 133     mpz_mul_ui(mpq_denref(r1), mpq_denref(r0), (2 * k) * (2 * k - 1));
 134     mpz_add_ui(mpq_numref(r1), mpq_numref(r1), 1);
 135     k++;
 136     mpz_mul_ui(mpq_numref(r0), mpq_numref(r1), (2 * k) * (2 * k - 1));
 137     mpz_mul_ui(mpq_denref(r0), mpq_denref(r1), (2 * k) * (2 * k - 1));
 138     mpz_sub_ui(mpq_numref(r0), mpq_numref(r0), 1);
 139   }
 140   while (cf_wait(cf)) {
 141     for(;;) {
 142       // Find largest x such that (p1 x + p0)/(q1 x + q0)
 143       // <= r0 for odd i, and >= r1 for even i
 144       //  =>  x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
 145       // for the appropriate r.
 146       mpq_ptr r = i & 1 ? r0 : r1;
 147       mpz_mul(t0, q0, mpq_numref(r));
 148       mpz_mul(t1, p0, mpq_denref(r));
 149       mpz_sub(t2, t0, t1);
 150       mpz_mul(t0, p1, mpq_denref(r));
 151       mpz_mul(t1, q1, mpq_numref(r));
 152       mpz_sub(t1, t0, t1);
 153       mpz_fdiv_q(t0, t2, t1);
 154
 155       // Then x is the next convergent provided substituting x + 1 overshoots
 156       // the other bound. (If not, we need better bounds to decide.)
 157       if (!mpz_sgn(t0)) {
 158         more_taylor();
 159       } else {
 160         mpz_add_ui(t1, t0, 1);
 161         mpz_set(mpq_numref(tq0), p0);
 162         mpz_addmul(mpq_numref(tq0), p1, t1);
 163         mpz_set(mpq_denref(tq0), q0);
 164         mpz_addmul(mpq_denref(tq0), q1, t1);
 165         if (i & 1) {
 166           if (mpq_cmp(tq0, r1) > 0) break;
 167           // Check (p1(x+1) + p0)/(q1(x+1) + q0) > r1
 168         } else {
 169           if (mpq_cmp(tq0, r0) < 0) break;
 170           // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
 171         }
 172         more_taylor();
 173       }
 174     };
 175     cf_put(cf, t0);
 176     i++;
 177     mpz_addmul(p0, p1, t0);
 178     mpz_addmul(q0, q1, t0);
 179     mpz_swap(p1, p0);
 180     mpz_swap(q1, q0);
 181   }
 182
 183   mpq_clear(r0); mpq_clear(r1);
 184   mpz_clear(p0); mpz_clear(p1); mpz_clear(q0); mpz_clear(q1);
 185   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 186   mpq_clear(tq0);
 187   return NULL;
 188 }
 189
 190 cf_t cf_new_cos1() {
 191   return cf_new_const(cos1_expansion);
 192 }