# Gamma function computation efficiency?

I wonder what kind of algorithm is used to compute the values for the gamma function. Specifically, I am interested in how the computational load increases when the complexity of the input grows. So, for example, evaluating

Gamma[x + I y]


for integer x and y and gradually increasing their number of digits will make the computation time longer. What is the corresponding power law for the increase?

EDIT: evaluating at 4 Pi + I 8 GoldenRatio:

data = Monitor[
Table[{k,
N[Gamma[4 Pi + I 8 GoldenRatio], k 10^3]; // Timing //
First}, {k, 1, 10, 5/10}], k];
ListLinePlot[data, Mesh -> All,
MeshStyle -> Directive[Red, PointSize[Large]]]


Now shows

and I get

line = Fit[Log@data, {1, x}, x]
-3.350 + 2.723 x


not sure if this is relevant or just an artifact.

-
Try Exit or Quit before you run the benchmark. It will quite the kernel. Als jus Gamma[1 + I ] will work. – Vitaliy Kaurov Dec 5 '13 at 4:50
Thank you! Result seems to be fairly consistent across different points in complex plane. Good to know that it is polynomial. – Kagaratsch Dec 5 '13 at 4:55

From this Documentation page: Some Notes on Internal Implementation

• For machine precision most special functions use Mathematica-derived rational minimax approximations. The notes that follow apply mainly to arbitrary precision.
• Gamma uses recursion, functional equations, and the Binet asymptotic formula.

On my machine the power exponent is 2.634. I wonder how much it depends on the setup. Benchmark it in the same number range and let me know.

data = Monitor[
Table[{k,
N[Gamma[Pi + I GoldenRatio], k 10^3]; // Timing // First},
{k, 1, 10, 5/10}], k];

ListLinePlot[data, Mesh -> All,
MeshStyle -> Directive[Red, PointSize[Large]]]


line = Fit[Log@data, {1, x}, x]
(* -2.672 + 2.634 x *)


It is true I think that Pi and GoldenRatio added a bit complexity. So for just Gamma[1 + I ] I got exponent 2.546 consistently.

line = Fit[Log@data, {1, x}, x]
(* -2.720 + 2.546 x *)

Show[ListPlot[Log@data, PlotStyle -> Red], Plot[line, {x, 0, 5}]]


-
I am not sure that I understand what the line Gamma[Pi + I GoldenRatio] in definition for data does. Are you evaluating gamma several times at the same point? – Kagaratsch Dec 5 '13 at 4:19
@Kagaratsch yes but with increasing precision driven by growing k*10^3. Mathematica has arbitrary precision, so I picked two numbers that have infinite aperiodic decimal part. – Vitaliy Kaurov Dec 5 '13 at 4:23
I get -2.409 + 2.677 x for line = Fit[Log@data, {1, x}, x]. Looks fairly consistent. And I am currently on a laptop. – Kagaratsch Dec 5 '13 at 4:31
@Kagaratsch Hah, interesting, I would be thrilled to see other folks let us know the exponent. Some neat experiments. – Vitaliy Kaurov Dec 5 '13 at 4:36
In a small machine i.stack.imgur.com/CEOIa.png, -3.70 + 3.761 x – Dr. belisarius Dec 5 '13 at 5:46

I'm not sure you are correct about it taking longer and longer. Here's a test which takes the first 1000 terms 2^x insider the gamma function and takes the timing:

all = N[Timing[Gamma[2^# + I 2^#]]] & /@ Range[1000];
ListPlot[Transpose[all][[1]]]


This does not show a generally increasing time.

-
I think in this case precision is too small to produce TIming growth above error fluctuations. – Vitaliy Kaurov Dec 5 '13 at 4:20
Looks like you are right, but what did Vitaliy Kaurov do different? – Kagaratsch Dec 5 '13 at 4:21
@Kagaratsch I used many more digits in precision so times will become significant, above error. See my comment right above yours. – Vitaliy Kaurov Dec 5 '13 at 4:25
It looks like Vitaly Kaurov's answer measures the calculation of gamma as the number of desired digits in the answer grows, whereas the above calculates the gamma values for larger values of the arguments. – bill s Dec 5 '13 at 4:26
Oh, yes, I guess going with just more digits he has a clever way to avoid an overflow and probe with much larger input. – Kagaratsch Dec 5 '13 at 4:28