famous.c

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