taylor.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <gmp.h>
   4
   5 int main(int argc, char *argv[]) {
   6   int k = 5;
   7   if (argc > 1) k = atoi(argv[1]);
   8
   9   mpz_t r, s;
  10   mpz_init(r); mpz_init(s);
  11
  12   mpz_set_ui(r, 1);
  13   mpz_set_ui(s, 1);
  14   int i;
  15   for (i = 1; i <= k ; i++ ) {
  16     mpz_mul_ui(s, s, (2 * i) * (2 * i + 1));
  17     mpz_mul_ui(r, r, (2 * i) * (2 * i + 1));
  18     if (i & 1) {
  19       mpz_sub_ui(r, r, 1);
  20     } else {
  21       mpz_add_ui(r, r, 1);
  22     }
  23   }
  24   gmp_printf("%Zd / %Zd\n", r, s);
  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   gmp_printf("%Zd\n", limit);
  32
  33   mpz_t p0, p1;
  34   mpz_t q0, q1;
  35   mpz_init(p0); mpz_init(p1);
  36   mpz_init(q0); mpz_init(q1);
  37   mpz_set_ui(p1, 1);
  38   mpz_set_ui(q0, 1);
  39
  40   mpz_t a, b;
  41   mpz_init(a); mpz_init(b);
  42   mpz_t t0, t1;
  43   mpz_init(t0); mpz_init(t1);
  44
  45   mpz_set(a, r);
  46   mpz_set(b, s);
  47   for (;;) {
  48     mpz_fdiv_qr(t0, t1, a, b);
  49
  50     mpz_addmul(p0, p1, t0);
  51     mpz_addmul(q0, q1, t0);
  52     mpz_mul(a, q0, q0);
  53     if (mpz_cmp(a, limit) >= 0) {
  54       printf("limit exceeded\n"); break;
  55     }
  56     mpz_swap(p0, p1);
  57     mpz_swap(q0, q1);
  58     gmp_printf("%Zd (%Zd) -> %Zd / %Zd\n", t0, t1, p1, q1);
  59     if (!mpz_sgn(t1)) break;
  60     mpz_set(a, b);
  61     mpz_set(b, t1);
  62   }
  63
  64   mpz_clear(limit);
  65   mpz_clear(r); mpz_clear(s);
  66   mpz_clear(t0); mpz_clear(t1);
  67   mpz_clear(a); mpz_clear(b);
  68   return 0;
  69 }