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   mpz_t odd, two;
  92   mpz_init(odd); mpz_init(two); mpz_init(minusz);
  93   mpz_set_ui(odd, 1);
  94   mpz_set_ui(two, 2);
  95   mpz_neg(minusz, z);
  96   cf_put(cf, odd);
  97   while(cf_wait(cf)) {
  98     cf_put(cf, z);
  99     cf_put(cf, odd);
 100     mpz_add_ui(odd, odd, 2);
 101     cf_put(cf, minusz);
 102     cf_put(cf, two);
 103   }
 104   mpz_clear(odd); mpz_clear(two); mpz_clear(minusz);
 105   mpz_clear(z);
 106   free(z);
 107   return NULL;
 108 }
 109
 110 // tanh n = z/(1 + z^2/(3 + z^2/(5 + z^2/...)))
 111 static void *gauss_tanh_expansion(cf_t cf) {
 112   mpz_ptr z = cf_data(cf);
 113   mpz_t z2;
 114   mpz_t odd;
 115   mpz_init(odd); mpz_init(z2);
 116   mpz_set_ui(odd, 1);
 117   mpz_set_ui(z2, 0);
 118   cf_put(cf, z2);
 119   cf_put(cf, z);
 120   mpz_mul(z2, z, z);
 121   while(cf_wait(cf)) {
 122     cf_put(cf, odd);
 123     mpz_add_ui(odd, odd, 2);
 124     cf_put(cf, z2);
 125   }
 126   mpz_clear(odd); mpz_clear(z2);
 127   mpz_clear(z);
 128   free(z);
 129   return NULL;
 130 }
 131
 132 // tan n = z/(1 - z^2/(3 - z^2/(5 - z^2/...)))
 133 static void *gauss_tan_expansion(cf_t cf) {
 134   mpz_ptr z = cf_data(cf);
 135   mpz_t z2;
 136   mpz_t odd;
 137   mpz_init(odd); mpz_init(z2);
 138   mpz_set_ui(odd, 1);
 139   mpz_set_ui(z2, 0);
 140   cf_put(cf, z2);
 141   cf_put(cf, z);
 142   mpz_mul(z2, z, z);
 143   mpz_neg(z2, z2);
 144   while(cf_wait(cf)) {
 145     cf_put(cf, odd);
 146     mpz_add_ui(odd, odd, 2);
 147     cf_put(cf, z2);
 148   }
 149   mpz_clear(odd); mpz_clear(z2);
 150   mpz_clear(z);
 151   free(z);
 152   return NULL;
 153 }
 154
 155 cf_t cf_new_one_arg(void *(*fun)(cf_t), mpz_t z) {
 156   mpz_ptr p = malloc(sizeof(*p));
 157   mpz_init(p);
 158   mpz_set(p, z);
 159   return cf_new(fun, p);
 160 }
 161
 162 struct funarg_s {
 163   void *(*fun)(cf_t);
 164   mpz_t arg;
 165 };
 166 typedef struct funarg_s *funarg_ptr;
 167
 168 static void *one_arg_nonreg(cf_t cf) {
 169   funarg_ptr p = cf_data(cf);
 170   mpz_t a, b, c, d;
 171   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
 172   mpz_set_ui(a, 1); mpz_set_ui(b, 0);
 173   mpz_set_ui(c, 0); mpz_set_ui(d, 1);
 174   mpz_ptr copy = malloc(sizeof(*copy));
 175   mpz_init(copy);
 176   mpz_set(copy, p->arg);
 177   cf_t nonreg = cf_new(p->fun, copy);
 178   cf_t conv = cf_new_nonregular_to_cf(nonreg, a, b, c, d);
 179   mpz_t z;
 180   mpz_init(z);
 181   while(cf_wait(cf)) {
 182     cf_get(z, conv);
 183     cf_put(cf, z);
 184   }
 185   mpz_clear(z);
 186   cf_free(conv);
 187   cf_free(nonreg);
 188   mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
 189
 190   mpz_clear(p->arg);
 191   free(p);
 192   return NULL;
 193 }
 194
 195 cf_t cf_new_one_arg_nonreg(void *(*fun)(cf_t), mpz_t z) {
 196   funarg_ptr p = malloc(sizeof(*p));
 197   p->fun = fun;
 198   mpz_init(p->arg);
 199   mpz_set(p->arg, z);
 200   return cf_new(one_arg_nonreg, p);
 201 }
 202
 203 cf_t cf_new_epow(mpz_t pow) {
 204   return cf_new_one_arg_nonreg(exp_expansion, pow);
 205 }
 206
 207 cf_t cf_new_tanh(mpz_t z) {
 208   return cf_new_one_arg_nonreg(gauss_tanh_expansion, z);
 209 }
 210
 211 // TODO: Handle negative convergents so this function works.
 212 cf_t cf_new_tan(mpz_t z) {
 213   return cf_new_one_arg_nonreg(gauss_tanh_expansion, z);
 214 }