Opened 10 years ago

Closed 10 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 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.

Attachments (1)

trac_10072-log_gamma.patch (2.8 KB) - added by fstan 10 years ago.

Download all attachments as: .zip

Change History (14)

comment:1 follow-up: Changed 10 years ago by 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

comment:2 in reply to: ↑ 1 Changed 10 years ago by 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.

comment:3 Changed 10 years ago by 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

comment:4 Changed 10 years ago by 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.

comment:5 Changed 10 years ago by 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.

comment:6 Changed 10 years ago by kcrisman

  • 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 10 years ago by 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

comment:8 follow-up: Changed 10 years ago by 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

comment:9 in reply to: ↑ 8 Changed 10 years ago by 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.

comment:10 Changed 10 years ago by zimmerma

  • Status changed from new to needs_work
  • Work issues set to 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 10 years ago by fstan

comment:11 Changed 10 years ago by fstan

  • Authors set to Flavia Stan
  • Status changed from needs_work to needs_review
  • Work issues changed from easy patch to write to patch to written

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

Flavia

comment:12 Changed 10 years ago by zimmerma

  • Reviewers set to Paul Zimmermann
  • Status changed from needs_review to positive_review
  • Work issues patch to written deleted

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

Paul

comment:13 Changed 10 years ago by jdemeyer

  • Merged in set to sage-4.6.1.alpha0
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.