famous.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <gmp.h>
   4 #include "cf.h"
   5
   6 // e = [2; 1, 2, 1, 1, 4, 1, ...]
   7 static void *e_expansion(cf_t cf) {
   8   int even = 2;
   9   cf_put_int(cf, even);
  10
  11   while(cf_wait(cf)) {
  12     cf_put_int(cf, 1);
  13     cf_put_int(cf, even);
  14     even += 2;
  15     cf_put_int(cf, 1);
  16   }
  17   return NULL;
  18 }
  19
  20 cf_t cf_new_e() {
  21   return cf_new_const(e_expansion);
  22 }
  23
  24 // 4/pi = 1 + 1/(3 + 4/(5 + 9/(7 + 16/(9 + ...))))
  25 static void *pi_arctan_sequence(cf_t cf) {
  26   mpz_t num, denom;
  27   mpz_init(num);
  28   mpz_init(denom);
  29
  30   mpz_set_ui(denom, 1);
  31   cf_put(cf, denom);
  32   mpz_set_ui(num, 1);
  33   cf_put(cf, num);
  34   while(cf_wait(cf)) {
  35     mpz_add_ui(denom, denom, 2);
  36     cf_put(cf, denom);
  37     mpz_add(num, num, denom);
  38     cf_put(cf, num);
  39   }
  40
  41   mpz_clear(num);
  42   mpz_clear(denom);
  43   return NULL;
  44 }
  45
  46 static void *regularized_pi(cf_t cf) {
  47   mpz_t a, b, c, d;
  48   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
  49   mpz_set_ui(a, 0); mpz_set_ui(b, 4);
  50   mpz_set_ui(c, 1); mpz_set_ui(d, 0);
  51   cf_t nonregpi = cf_new(pi_arctan_sequence, NULL);
  52   cf_t conv = cf_new_nonregular_to_cf(nonregpi, a, b, c, d);
  53   mpz_t z;
  54   mpz_init(z);
  55   while(cf_wait(cf)) {
  56     cf_get(z, conv);
  57     cf_put(cf, z);
  58   }
  59   mpz_clear(z);
  60   cf_free(conv);
  61   cf_free(nonregpi);
  62   mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
  63   return NULL;
  64 }
  65
  66 cf_t cf_new_pi() {
  67   return cf_new_const(regularized_pi);
  68 }
  69
  70 // tan 1 = [1; 1, 1, 3, 1, 5, ...]
  71 static void *tan1_expansion(cf_t cf) {
  72   int odd = 1;
  73
  74   while(cf_wait(cf)) {
  75     cf_put_int(cf, 1);
  76     cf_put_int(cf, odd);
  77     odd += 2;
  78   }
  79
  80   return NULL;
  81 }
  82
  83 cf_t cf_new_tan1() {
  84   return cf_new_const(tan1_expansion);
  85 }
  86
  87 // exp(z) = 1 + z/(1 - z/(2 + z/(3 - z/(2 + z/(5 - z/(2 + z/ ...))))))
  88 void *exp_expansion(cf_t cf) {
  89   mpz_ptr z = cf_data(cf);
  90   mpz_t minusz;
  91   int odd = 1;
  92   mpz_init(minusz);
  93   mpz_neg(minusz, z);
  94   cf_put_int(cf, 1);
  95   while(cf_wait(cf)) {
  96     cf_put(cf, z);
  97     cf_put_int(cf, odd);
  98     odd += 2;
  99     cf_put(cf, minusz);
 100     cf_put_int(cf, 2);
 101   }
 102   mpz_clear(minusz);
 103   mpz_clear(z);
 104   free(z);
 105   return NULL;
 106 }
 107
 108 // tanh n = z/(1 + z^2/(3 + z^2/(5 + z^2/...)))
 109 static void *gauss_tanh_expansion(cf_t cf) {
 110   mpz_ptr z = cf_data(cf);
 111   mpz_t z2;
 112   mpz_t odd;
 113   mpz_init(odd); mpz_init(z2);
 114   mpz_set_ui(odd, 1);
 115   mpz_set_ui(z2, 0);
 116   cf_put(cf, z2);
 117   cf_put(cf, z);
 118   mpz_mul(z2, z, z);
 119   while(cf_wait(cf)) {
 120     cf_put(cf, odd);
 121     mpz_add_ui(odd, odd, 2);
 122     cf_put(cf, z2);
 123   }
 124   mpz_clear(odd); mpz_clear(z2);
 125   mpz_clear(z);
 126   free(z);
 127   return NULL;
 128 }
 129
 130 // tan n = z/(1 - z^2/(3 - z^2/(5 - z^2/...)))
 131 static void *gauss_tan_expansion(cf_t cf) {
 132   mpz_ptr z = cf_data(cf);
 133   mpz_t z2;
 134   mpz_t odd;
 135   mpz_init(odd); mpz_init(z2);
 136   mpz_set_ui(odd, 1);
 137   mpz_set_ui(z2, 0);
 138   cf_put(cf, z2);
 139   cf_put(cf, z);
 140   mpz_mul(z2, z, z);
 141   mpz_neg(z2, z2);
 142   while(cf_wait(cf)) {
 143     cf_put(cf, odd);
 144     mpz_add_ui(odd, odd, 2);
 145     cf_put(cf, z2);
 146   }
 147   mpz_clear(odd); mpz_clear(z2);
 148   mpz_clear(z);
 149   free(z);
 150   return NULL;
 151 }
 152
 153 cf_t cf_new_one_arg(void *(*fun)(cf_t), mpz_t z) {
 154   mpz_ptr p = malloc(sizeof(*p));
 155   mpz_init(p);
 156   mpz_set(p, z);
 157   return cf_new(fun, p);
 158 }
 159
 160 struct funarg_s {
 161   void *(*fun)(cf_t);
 162   mpz_t arg;
 163 };
 164 typedef struct funarg_s *funarg_ptr;
 165
 166 static void *one_arg_nonreg(cf_t cf) {
 167   funarg_ptr p = cf_data(cf);
 168   mpz_t a, b, c, d;
 169   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
 170   mpz_set_ui(a, 1); mpz_set_ui(b, 0);
 171   mpz_set_ui(c, 0); mpz_set_ui(d, 1);
 172   mpz_ptr copy = malloc(sizeof(*copy));
 173   mpz_init(copy);
 174   mpz_set(copy, p->arg);
 175   cf_t nonreg = cf_new(p->fun, copy);
 176   cf_t conv = cf_new_nonregular_to_cf(nonreg, a, b, c, d);
 177   mpz_t z;
 178   mpz_init(z);
 179   while(cf_wait(cf)) {
 180     cf_get(z, conv);
 181     cf_put(cf, z);
 182   }
 183   mpz_clear(z);
 184   cf_free(conv);
 185   cf_free(nonreg);
 186   mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
 187
 188   mpz_clear(p->arg);
 189   free(p);
 190   return NULL;
 191 }
 192
 193 cf_t cf_new_one_arg_nonreg(void *(*fun)(cf_t), mpz_t z) {
 194   funarg_ptr p = malloc(sizeof(*p));
 195   p->fun = fun;
 196   mpz_init(p->arg);
 197   mpz_set(p->arg, z);
 198   return cf_new(one_arg_nonreg, p);
 199 }
 200
 201 cf_t cf_new_epow(mpz_t pow) {
 202   return cf_new_one_arg_nonreg(exp_expansion, pow);
 203 }
 204
 205 cf_t cf_new_tanh(mpz_t z) {
 206   return cf_new_one_arg_nonreg(gauss_tanh_expansion, z);
 207 }
 208
 209 cf_t cf_new_tan(mpz_t z) {
 210   return cf_new_one_arg_nonreg(gauss_tan_expansion, z);
 211 }