Bug in log gamma evaluation

[kcrisman]

When log_gamma is naively numerically evaluated, it seems to be wrong.

sage: log_gamma(i).n()
-0.0219785471672303 - 0.168184318273662*I

Pari gives this:

sage: log_gamma(pari(i))
-0.650923199301856 - 1.87243664726243*I

And in Sage's Maxima:

(%i3) ev(log_gamma(%i),numer);
(%o3)             - 1.872436647262428 %i - .6509231993018556

And in ginsh for Ginac 1.5.8:

> lgamma(I);
-0.65092319930185634056+4.41074865991715666*I

which is pretty clearly a different branch choice from Maxima and Pari (the complex part is exactly 2*pi away from Maxima and Pari).

The relevant bit of code is

sage: a = log_gamma(i)
sage: a._convert(RealField(53).complex_field())
-0.0219785471672303 - 0.168184318273662*I

But I'm not sure if this problem is in Pynac or in the Sage complex field.

[zimmerma]

But I'm not sure if this problem is in Pynac or in RealField?.

as it name says, RealField? only deals with real values, not complex values.

Paul

[kcrisman]

Replying to zimmerma:

But I'm not sure if this problem is in Pynac or in RealField?.

as it name says, RealField? only deals with real values, not complex values.

Okay, fair enough - so change that to whatever RealField.complex_field() becomes, which I don't know. Feel free to alter the description to wherever RealField goes under that - ComplexField, I guess.

[zimmerma]

what is strange is that the precision is not taken into account:

sage: a=log_gamma(I)               
sage: a._convert(ComplexField(53))
-0.0219785471672303 - 0.168184318273662*I
sage: a._convert(ComplexField(100))
-0.0219785471672303 - 0.168184318273662*I

Paul

[kcrisman]

Which is why I think it's a bug in Pynac/GiNaC. I am not sure how to check this, though, since _convert uses the mysterious attribute _gobj which is not accessible from the outside world.

[kcrisman]

Though the fact that log_gamma is not currently symbolic (#10075) could have something to do with this. I'm sure Burcin will eventually weigh in, but I know he's probably pretty busy right now, since he worked hard on getting the new Pynac ready just before this.

[kcrisman]

More mysterious (after MUCH trouble trying to get Ginac to compile on OS X - needed to install CLN, pkg-config, and readline!!!):

ginsh - GiNaC Interactive Shell (ginac V1.5.8)
  __,  _______  Copyright (C) 1999-2010 Johannes Gutenberg University Mainz,
 (__) *       | Germany.  This is free software with ABSOLUTELY NO WARRANTY.
  ._) i N a C | You are welcome to redistribute it under certain conditions.
<-------------' For details type `warranty;'.

Type ?? for a list of help topics.
> lgamma(I);
-0.65092319930185634056+4.41074865991715666*I

[kcrisman]

And even more perplexing is that gamma(I) is the same in all four systems.

(%i1) ev(gamma(%i),numer);
(%o1)             - .4980156681183565 %i - .1549498283018101
sage: gamma(i).n()
-0.154949828301811 - 0.498015668118356*I
sage: pari(i).gamma()
-0.154949828301811 - 0.498015668118356*I
> tgamma(I);
-0.15494982830181068507-0.49801566811835604268*I

[zimmerma]

I don't see how the wrong result is related (if any) to the correct one. The difference in the imaginary part is 1.70425232898877, which is neither 2pi nor pi nor even pi/2.

Paul

[kcrisman]

Replying to zimmerma:

I don't see how the wrong result is related (if any) to the correct one. The difference in the imaginary part is 1.70425232898877, which is neither 2pi nor pi nor even pi/2.

Correct. My comment about the 2*pi was in reference to the Ginac output, not the Sage output. But I can't for the life of me figure out how Sage is actually doing log_gamma(i).n(), because the conversion to the complex field pretty clearly just converts it to a Ginac object and then numerically evaluates.

As to the precision, I think there must be something missing in our code, because the Ginac docs state

The function evalf that was used above converts any number in GiNaC's expressions into floating point numbers. This can be done to arbitrary predefined accuracy:

     > evalf(1/7);
     0.14285714285714285714
     > Digits=150;
     150
     > evalf(1/7);
     0.1428571428571428571428571428571428571428571428571428571428571428571428
     5714285714285714285714285714285714285

So maybe we're not taking that into account, though I have no idea how I would do so.

[zimmerma]

I found the reason of the bug. In symbolic/pynac.pyx, function py_lgamma, it is written CC(x).log().gamma() instead of CC(x).gamma().log():

sage: CC(I).log().gamma()
-0.0219785471672303 - 0.168184318273662*I
sage: CC(I).gamma().log()
-0.650923199301856 - 1.87243664726243*I

If somebody provides a patch, I will review it.

Paul

[fstan]

This patch should fix the py_lgamma function using as a model the py_log

Flavia

[zimmerma]

excellent work, Flavia! Not only it fixes the reported problem, but also now we can get arbitrary precision values.

Paul