famous.c

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