Karatsuba’s algorithm for the gamma function

September 17, 2012

Another brief illustration of computing things to very high precision using my new library Arb: I have implemented the asymptotically fast FEE algorithm of E. Karatsuba for evaluating $\Gamma(a)$ along with the first few derivatives $\Gamma’(a), \Gamma^{(2)}(a), \ldots$ at rational $a$.

As a benchmark, I evaluated the first 10 derivatives at $a = 666/1337$ to a precision of 100,000 digits. Computing just a single value with a method that is not asymptotically fast (e.g. Stirling’s series) would likely take several hours. The version of the algorithm described in Karatsuba’s paper (after fixing a small error), using exact binary splitting, takes 1135 seconds. After improving the parameter selection and enabling rounding in the binary splitting stage, it takes 287 seconds. That’s a nice 4x speedup, and by using ball arithmetic, we still get the output with a rigorous error bound.

Computing one million digits of $\Gamma(1/2) = \sqrt \pi$ takes 281 seconds with the original algorithm, and 88 seconds with the improvements enabled. That’s still 75 times slower than the Chudnovsky algorithm for $\pi$, but not too shabby for an algorithm that solves a much more general problem.

An outline of the algorithm can be found in the documentation.

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor