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 void cf_dump_term(cf_t cf, int n) {
  48   mpz_t z;
  49   mpz_init(z);
  50   int i;
  51
  52   gmp_printf("[%Zd", z);
  53   for (i = 1; i <= n; i++) {
  54     cf_get(z, cf);
  55     gmp_printf(" %Zd", z);
  56   }
  57   putchar(']');
  58   putchar('\n');
  59   mpz_clear(z);
  60 }
  61
  62 int main(int argc, char **argv) {
  63   int n = 100;
  64   if (argc > 1) {
  65     n = atoi(argv[1]);
  66     if (n <= 0) n = 100;
  67   }
  68   cf_t c[7];
  69   cf_t t[6][2];
  70   mpz_t b[8];
  71   int i;
  72
  73   mpz8_init(b);
  74   mpz8_set_int(b,
  75       2, 0, 0, -1,
  76       0, 0, 0, 1);
  77
  78   c[0] = cf_new_cos1();
  79   for (i = 1; i < 7; i++) {
  80     cf_tee(t[i - 1], c[i - 1]);
  81     c[i] = cf_new_bihom(t[i - 1][0], t[i - 1][1], b);
  82   }
  83   // c[6] = cos 64
  84
  85   cf_t c1 = cf_new_cos1();
  86   cf_t ct0[2], ct1[2], ct2[2];
  87   cf_tee(ct1, c1);
  88   cf_tee(ct2, ct1[1]);
  89   cf_t c2 = cf_new_mul(ct2[0], ct2[1]);
  90   cf_tee(ct0, c2);
  91   mpz8_set_int(b,
  92       16, -20, 0, 5,  // cos 5n (5th Chebyshev polynomial)
  93       0, 0, 0, 1);
  94   cf_t c3 = cf_new_bihom(ct0[0], ct0[1], b);
  95   cf_t c5 = cf_new_mul(ct1[0], c3);
  96   // c5 = cos 5
  97
  98   cf_t c64t[2], c64t1[2];
  99   cf_tee(c64t, c[6]);
 100   cf_tee(c64t1, c64t[1]);
 101   mpz8_set_int(b,
 102       -1, 0, 0, 1,
 103        0, 0, 0, 1);
 104   cf_t s1 = cf_new_bihom(c64t1[0], c64t1[1], b);
 105   cf_t s64 = cf_new_sqrt(s1);
 106   // s64 = sin 64, c64t[0] = untouched cos 64
 107
 108   cf_t c5t[2], c5t1[2];
 109   cf_tee(c5t, c5);
 110   cf_tee(c5t1, c5t[1]);
 111   cf_t s3 = cf_new_bihom(c5t1[0], c5t1[1], b);
 112   cf_t s5 = cf_new_sqrt(s3);
 113   // s5 = sin 5, c5t[0] = untouched cos 5
 114
 115   cf_t cf0 = cf_new_mul(s5, c64t[0]);
 116   cf_t cf1 = cf_new_mul(s64, c5t[0]);
 117   cf_t s69 = cf_new_sub(cf0, cf1); // TODO: Respect signs.
 118                                    // This is supposed to be an addition,
 119                                    // and the result should be negative.
 120   // s69 = sin 69
 121
 122   cf_t sqrt5 = cf_new_sqrt5();
 123   cf_t ts5d = cf_new_const_nonregular(tanhsqrt5denom);
 124   cf_t ts5 = cf_new_div(sqrt5, ts5d);
 125
 126   cf_t den = cf_new_add(ts5, s69); // TODO: This should be a subtract, but
 127                                    // again, we don't handle signs properly.
 128
 129   cf_t e = cf_new_e();
 130   cf_t pi = cf_new_pi();
 131   cf_t pitee[2];
 132   cf_tee(pitee, pi);
 133   mpz8_set_int(b,
 134       0, 0, 0, 3,
 135       1, 0, 0, 0);
 136   // tops = Three-Over-Pi-Squared.
 137   cf_t tops = cf_new_bihom(pitee[0], pitee[1], b);
 138
 139   mpz8_set_add(b);
 140   cf_t sum = cf_new_bihom(tops, e, b);
 141   cf_t num = cf_new_sqrt(sum);
 142   cf_t hakmem_constant = cf_new_div(num, den);
 143   cf_dump(hakmem_constant, n);
 144
 145   mpz8_clear(b);
 146   return 0;
 147 }