hakmem.c

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