Opened 12 years ago
Closed 12 years ago
#10072 closed defect (fixed)
Bug in log gamma evaluation
Reported by: | kcrisman | Owned by: | AlexGhitza |
---|---|---|---|
Priority: | major | Milestone: | sage-4.6.1 |
Component: | basic arithmetic | Keywords: | |
Cc: | burcin, zimmerma | Merged in: | sage-4.6.1.alpha0 |
Authors: | Flavia Stan | Reviewers: | Paul Zimmermann |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
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.
Attachments (1)
Change History (14)
comment:1 follow-up: 2 Changed 12 years ago by
comment:2 Changed 12 years ago by
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.
comment:3 Changed 12 years ago by
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
comment:4 Changed 12 years ago by
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.
comment:5 Changed 12 years ago by
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.
comment:6 Changed 12 years ago by
Description: | modified (diff) |
---|
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
comment:7 Changed 12 years ago by
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
comment:8 follow-up: 9 Changed 12 years ago by
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
comment:9 Changed 12 years ago by
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.
comment:10 Changed 12 years ago by
Status: | new → needs_work |
---|---|
Work issues: | → easy patch to write |
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
Changed 12 years ago by
Attachment: | trac_10072-log_gamma.patch added |
---|
comment:11 Changed 12 years ago by
Authors: | → Flavia Stan |
---|---|
Status: | needs_work → needs_review |
Work issues: | easy patch to write → patch to written |
This patch should fix the py_lgamma
function using as a model the py_log
Flavia
comment:12 Changed 12 years ago by
Reviewers: | → Paul Zimmermann |
---|---|
Status: | needs_review → positive_review |
Work issues: | patch to written |
excellent work, Flavia! Not only it fixes the reported problem, but also now we can get arbitrary precision values.
Paul
comment:13 Changed 12 years ago by
Merged in: | → sage-4.6.1.alpha0 |
---|---|
Resolution: | → fixed |
Status: | positive_review → closed |
as it name says, RealField? only deals with real values, not complex values.
Paul