2

I'm implementing an arbitrary precision arithmetic library in C++ and I'm pretty much stuck when implementing the gamma function.

By using the equivalences gamma(n) = gamma(n - 1) * n and gamma(n) = gamma(n + 1) / n, respectively, I can obtain a rational number r in the range (1; 2] for all real values x.

However, I don't know how to evaluate gamma(r). For the Lanczos approximation (https://en.wikipedia.org/wiki/Lanczos_approximation), I need precomputed values p which happen to calculate a factorial of a non-integer value (?!) and can't be calculated dynamically with my current knowledge... Precomputing values for p wouldn't make much sense when implementing an arbitrary precision library.

Are there any algorithms that compute gamma(r) in a reasonable amount of time with arbitrary precision? Thanks for your help.

Mark Dickinson
  • 29,088
  • 9
  • 83
  • 120
borchero
  • 5,562
  • 8
  • 46
  • 72
  • I presume you have studied http://functions.wolfram.com/GammaBetaErf/Gamma/introductions/Gamma/05/ and failed to find what you need. – High Performance Mark Nov 15 '16 at 17:27
  • Looks like it's more related to math than to programming... Maybe, it's better to post question there. – George Sovetov Nov 15 '16 at 17:33
  • I'd also look into MPFR's source to see what they use. – Mark Dickinson Nov 15 '16 at 17:34
  • 1
    @GeorgeSovetov: I think it's reasonably on-topic, given the practical constraints of an efficient implementation. – Mark Dickinson Nov 15 '16 at 17:35
  • 1
    Second 4.5 of Modern Computer Arithmetic (official PDF [here](https://www.google.co.uk/search?client=safari&rls=en&q=modern+computer+arithmetic&ie=UTF-8&oe=UTF-8&gfe_rd=cr&ei=B0grWITKAc3FaJb8i9AP)) gives one approach, which is to use Stirling's asymptotic expansion for large `x`, and the functional equation to transfer this to smaller `x`. (The reflection formula lets you extend this to negative `x` efficiently.) – Mark Dickinson Nov 15 '16 at 17:39
  • BTW, I just checked the MPFR source: it uses exactly the approach described above (Stirling's expansion + functional equation + reflection formula). Not a huge surprise, given that Paul Zimmerman is a coauthor of both MPFR and MCA. Check out `src/lngamma.c` in the MPFR source for details. – Mark Dickinson Nov 15 '16 at 17:51

2 Answers2

4

Spouge's approximation is similar to Lanczos's approximation, but probably easier to use for arbitrary precision, as you can set the desired error.

Simon Byrne
  • 7,694
  • 1
  • 26
  • 50
1

Lanczos approximation doesn't seem too bad. What exactly do you suspect?

Parts of code which calculate p, C (Chebyshev polynomials) and (a + 1/2)! can be implemented as stateful objects so that, for example, you can calculate p(i) from p(i-1) and Chebyshev coefficients and be computed once, maintaining their matrix.

George Sovetov
  • 4,942
  • 5
  • 36
  • 57