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) |

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