fredrikj.net / blog /

An even speedier gamma function

February 27, 2013

I’ve just pushed some new code to the Arb git repository for computing numerical rising factorials ($x (x+1) (x+2)\cdots (x+n-1)$) faster. It works by expanding subproducts as symbolic polynomials, and evaluating them using the rectangular splitting algorithm mentioned in the previous post about computing roots of unity.

At precisions in the thousands of digits, this ultimately results in the gamma function becoming almost twice as fast!

Here is a comparison for the complex gamma function, evaluating $\Gamma(z)$ where $z = \sqrt 2 + i \sqrt 3$. The table is the same as in the previous post about the complex gamma function, but with a new column on the right (all timings are in seconds and a number in parentheses measures the first evaluation, before coefficients have been cached):

Digits Pari/GP mpmath Mathematica Arb (old) Arb (new)
10 0.000032 0.00013 0.00012 0.000077 0.000067
30 0.000064 0.00017 0.00021 0.00013 0.00013
100 0.00018 0.00043 0.00052 0.00030 0.00029
300 0.0010 0.0023 0.0022 0.00095 0.0010
1000 0.013 (0.26) 0.026 (0.37) 0.020 (0.10) 0.0080 (0.010) 0.0070 (0.010)
3000 0.19 (6.5) 0.36 (3.6) 0.31 (1.7) 0.10 (0.18) 0.068 (0.14)
10000 2.9 (235) 6.4 (95) 5.8 (44) 1.6 (3.1) 0.96 (2.4)
30000 38 (6030) 92 (3030) 80 (887) 22 (45) 11 (32)

Finally, here also is a comparison for the real gamma function, evaluating $\Gamma(x)$ where $x = \sqrt 2$. I have included MPFR instead of mpmath (of course, MPFR doesn’t have a complex gamma function, so I couldn’t include it above):

Digits Pari/GP MPFR Mathematica Arb (old) Arb (new)
30 0.000024 0.000090 0.000070 0.000050 0.000043
100 0.000068 0.00036 0.00020 0.00011 0.00011
300 0.00039 0.0028 0.00080 0.00036 0.00036
1000 0.0046 (0.25) 0.046 0.058 0.0028 0.0025
3000 0.12 (6.5) 1.2 0.76 0.031 (0.12) 0.024 (0.090)
10000 1.9 (233) 60 13 0.49 (2.1) 0.29 (1.5)
30000 13 (6154) 2680 186 6.6 (32) 3.4 (22)