pi.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <gmp.h>
   4 #include "cf.h"
   5
   6 static void *pi_arctan_sequence(cf_t cf) {
   7   mpz_t num, denom;
   8   mpz_init(num);
   9   mpz_init(denom);
  10
  11   mpz_set_ui(denom, 1);
  12   cf_put(cf, denom);
  13
  14   mpz_set_ui(num, 1);
  15   cf_put(cf, num);
  16
  17   while(cf_wait(cf)) {
  18     mpz_add_ui(denom, denom, 2);
  19     cf_put(cf, denom);
  20     mpz_add(num, num, denom);
  21     cf_put(cf, num);
  22   }
  23
  24   mpz_clear(num);
  25   mpz_clear(denom);
  26   return NULL;
  27 }
  28
  29 int main() {
  30   mpz_t z;
  31   mpz_t a, b, c, d;
  32   mpz_init(z);
  33   cf_t pi, conv;
  34   pi = cf_new(pi_arctan_sequence, NULL);
  35   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
  36   mpz_set_ui(a, 0); mpz_set_ui(b, 4);
  37   mpz_set_ui(c, 1); mpz_set_ui(d, 0);
  38
  39   conv = cf_new_nonregular_to_cf(pi, a, b, c, d);
  40   cf_t dec = cf_new_cf_to_decimal(conv);
  41   // Interestingly, regularizing the continued fraction and
  42   // then converting to decimal is faster than doing everything at once:
  43   //   conv = cf_new_nonregular_mobius_to_decimal(pi, a, b, c, d);
  44   for (int i = 1; i <= 5000; i++) {
  45     cf_get(z, dec);
  46     gmp_printf("%Zd",z);
  47     if (!(i%5)) printf(" ");
  48     if (!(i%50)) printf("\n");
  49   }
  50   printf("\n");
  51   cf_free(dec);
  52   cf_free(conv);
  53   cf_free(pi);
  54   mpz_clear(z);
  55   mpz_clear(d); mpz_clear(c); mpz_clear(b); mpz_clear(a);
  56   return 0;
  57 }