Opened 12 years ago

Closed 12 years ago

# Bug in log gamma evaluation

Reported by: Owned by: kcrisman AlexGhitza major sage-4.6.1 basic arithmetic burcin, zimmerma sage-4.6.1.alpha0 Flavia Stan Paul Zimmermann N/A

### GitHub link to the corresponding issue

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.

### comment:1 follow-up:  2 Changed 12 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 12 years ago by kcrisman

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 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 12 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 12 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 12 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 12 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:  9 Changed 12 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 12 years ago by kcrisman

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 zimmerma

Status: new → needs_work → 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

### comment:11 Changed 12 years ago by fstan

Authors: → Flavia Stan needs_work → needs_review 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 zimmerma

Reviewers: → Paul Zimmermann needs_review → positive_review 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 jdemeyer

Merged in: → sage-4.6.1.alpha0 → fixed positive_review → closed
Note: See TracTickets for help on using tickets.