hakmem.c

   1 // Compute the HAKMEM constant.
   2 //   (sqrt(3/pi^2 + e)/(tanh(sqrt(5))-sin(69))
   3 //
   4 // Confirm with:
   5 //  $ echo "pi=4*a(1)
   6 //          x=2*sqrt(5)
   7 //          (sqrt(3/pi^2+e(1)))/((e(x)-1)/(e(x)+1)-s(69))" | bc -l
   8 //
   9 // Almost exclusively uses techniques described by Gosper. Differences:
  10 //
  11 //  - Taylor series for cos(1) to find its continued fraction expansion.
  12 //  - Binary search instead of Newton's method when finding integer part of
  13 //  solution of quadratic.
  14
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17 #include <string.h>
  18 #include <gmp.h>
  19 #include "cf.h"
  20
  21 void *tanhsqrt5denom(cf_t cf) {
  22   int odd = 1;
  23   while(cf_wait(cf)) {
  24     cf_put_int(cf, odd);
  25     odd += 2;
  26     cf_put_int(cf, 5);
  27   }
  28   return NULL;
  29 }
  30
  31 void cf_dump(cf_t cf, int n) {
  32   mpz_t z;
  33   mpz_init(z);
  34   cf_t dec = cf_new_cf_to_decimal(cf);
  35   int i;
  36   for (i = 0; i <= n; i++) {
  37     cf_get(z, dec);
  38     gmp_printf("%Zd", z);
  39     if (!(i % 5)) putchar(' ');
  40     if (!(i % 50)) putchar('\n');
  41   }
  42   if (n % 50) putchar('\n');
  43   cf_free(dec);
  44   mpz_clear(z);
  45 }
  46
  47 int main(int argc, char **argv) {
  48   int n = 100;
  49   if (argc > 1) {
  50     n = atoi(argv[1]);
  51     if (n <= 0) n = 100;
  52   }
  53   cf_t c[7];
  54   cf_t t[6][2];
  55   mpz_t b[8];
  56   int i;
  57
  58   mpz8_init(b);
  59   mpz8_set_int(b,
  60       2, 0, 0, -1,
  61       0, 0, 0, 1);
  62
  63   c[0] = cf_new_cos1();
  64   for (i = 1; i < 7; i++) {
  65     cf_tee(t[i - 1], c[i - 1]);
  66     c[i] = cf_new_bihom(t[i - 1][0], t[i - 1][1], b);
  67   }
  68   // c[6] = cos 64
  69
  70   cf_t c1 = cf_new_cos1();
  71   cf_t ct0[2], ct1[2], ct2[2];
  72   cf_tee(ct1, c1);
  73   cf_tee(ct2, ct1[1]);
  74   cf_t c2 = cf_new_mul(ct2[0], ct2[1]);
  75   cf_tee(ct0, c2);
  76   mpz8_set_int(b,
  77       16, -20, 0, 5,  // cos 5n (5th Chebyshev polynomial)
  78       0, 0, 0, 1);
  79   cf_t c3 = cf_new_bihom(ct0[0], ct0[1], b);
  80   cf_t c5 = cf_new_mul(ct1[0], c3);
  81   // c5 = cos 5
  82
  83   cf_t c64t[2], c64t1[2];
  84   cf_tee(c64t, c[6]);
  85   cf_tee(c64t1, c64t[1]);
  86   mpz8_set_int(b,
  87       -1, 0, 0, 1,
  88        0, 0, 0, 1);
  89   cf_t s1 = cf_new_bihom(c64t1[0], c64t1[1], b);
  90   cf_t s64 = cf_new_sqrt(s1);
  91   // s64 = sin 64, c64t[0] = untouched cos 64
  92
  93   cf_t c5t[2], c5t1[2];
  94   cf_tee(c5t, c5);
  95   cf_tee(c5t1, c5t[1]);
  96   cf_t s3 = cf_new_bihom(c5t1[0], c5t1[1], b);
  97   cf_t s5 = cf_new_sqrt(s3);
  98   // s5 = sin 5, c5t[0] = untouched cos 5
  99
 100   cf_t cf0 = cf_new_mul(s5, c64t[0]);
 101   cf_t cf1 = cf_new_mul(s64, c5t[0]);
 102   cf_t s69 = cf_new_sub(cf0, cf1); // TODO: Respect signs.
 103                                    // This is supposed to be an addition,
 104                                    // and the result should be negative.
 105   // s69 = sin 69
 106
 107   cf_t sqrt5 = cf_new_sqrt5();
 108   cf_t ts5d = cf_new_const_nonregular(tanhsqrt5denom);
 109   cf_t ts5 = cf_new_div(sqrt5, ts5d);
 110
 111   cf_t den = cf_new_add(ts5, s69); // TODO: This should be a subtract, but
 112                                    // again, we don't handle signs properly.
 113
 114   cf_t e = cf_new_e();
 115   cf_t pi = cf_new_pi();
 116   cf_t pitee[2];
 117   cf_tee(pitee, pi);
 118   mpz8_set_int(b,
 119       0, 0, 0, 3,
 120       1, 0, 0, 0);
 121   // tops = Three-Over-Pi-Squared.
 122   cf_t tops = cf_new_bihom(pitee[0], pitee[1], b);
 123
 124   mpz8_set_add(b);
 125   cf_t sum = cf_new_bihom(tops, e, b);
 126   cf_t num = cf_new_sqrt(sum);
 127   cf_t hakmem_constant = cf_new_div(num, den);
 128   cf_dump(hakmem_constant, n);
 129
 130   mpz8_clear(b);
 131   return 0;
 132 }