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