The HAKMEM Constant

HAKMEM features the expression:

\[ \frac{\sqrt{3/\pi^2 + e}}{\tanh(\sqrt{5}) - \sin(69)} \]

and states it is easily computed using continued fractions.

Observe \(\tanh\sqrt{5}\) can be found with one division, one square root, and one application of the continued fraction of Gauss:

\[ \tanh\sqrt{5} = \frac{\sqrt{5}}{1 + \frac{5}{3 + \frac{5}{5 + ...}}} \]

The most troublesome subexpression is \(\sin(69)\), for the \(\tan\) continued fraction expansion behaves poorly, and the terms in the Taylor series for \(\sin(69)\) about 0 grow for quite some time before shrinking to acceptable levels for continued fractions.

As previously suggested, one possibility is to compute the continued fraction expansion of \(\cos 1\) using the Taylor series, then evaluate Chebyshev polynomials and other trigonometric identities to find \(\sin 69\).

I followed this route in a test program for the frac library. Starting from \(\cos 1\), I used double-angle formulas to find \(\cos 64\), and the fifth Chebyshev polynomial (of the first kind) to find \(\cos 5\). Then since we may square root continued fractions, we can use \(\sin x = \sqrt{1 - \cos ^2 x}\) to find \(\sin 64\) and \(\sin 5\). Lastly, the angle addition formula yields \(\sin 69\).

Thus using techniques Gosper describes, and the Taylor series, we find the first 100 decimal places of the HAKMEM constant:

1
59170 96974 31217 53554 22849 04695 38245 87294 24160 11857
90051 76587 45519 76440 69814 05924 94706 86439 45987 93518

and the first 100 terms of its continued fraction expansion:

\[ [0; 1, 1, 1, 2, 4, 2, 2, 1, 4, 1, 6, 2, 9, 13, 1, 1, 8, 3, 7, 1, 10, 6, 11, 2, 2, 2, 3, 3, 6, 12, 2, 1, 12, 1, 2, 2, 5, 1, 2, 1, 1, 24, 6, 1, 1, 2, 11, 1, 8, 50, 1, 1, 2, 3, 1, 2, 1, 130, 1, 29, 3, 1, 9, 1, 1, 1, 1, 16, 1, 2, 1, 1, 4, 1, 3, 1, 1, 7, 1, 1, 5, 1, 1168, 1, 6, 1, 10, 1, 3, 1, 8, 1, 12, 21, 1, 1, 7, 1, 1, 1, ...] \]

For computing 1000 digits of the HAKMEM constant, my program was twice as fast as running these bc -l commands:

scale=1000
pi=4*a(1)
x=2*sqrt(5)
(sqrt(3/pi^2+e(1)))/((e(x)-1)/(e(x)+1)-s(69))

This is gratifying, especially as I could easily optimize further. (For example, I should be able to remove three integer divisions in my quadratic algorithm implementation; better ways of finding \(\sin 69\) exist).


Ben Lynn blynn@cs.stanford.edu 💡