Opened 7 years ago

Closed 5 years ago

# numerical bug in incomplete gamma function

Reported by: Owned by: VivianePons major sage-duplicate/invalid/wontfix symbolics sd67 Fixed upstream, in a later stable release.

Wrong values from incomplete gamma:

```\$ ./sage -v
SageMath Version 6.7.beta1, Release Date: 2015-04-15

\$ ./sage
sage: gamma(60, 30).numerical_approx()
-1.28306738270893e74
```

Reason is a bug in Pari:

```? \p 50
realprecision = 57 significant digits (50 digits displayed)
? incgam(60,30)
%2 = 1.3868299023788799504133735635863795825413935892945 E80
? \p 30
realprecision = 38 significant digits (30 digits displayed)
? incgam(60,30)
%3 = -1.23084064495468737276106696496 E74
```

Previous description:

The following code:

```sage: var('R k')
(R, k)
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41), (R,0,6))
```

gives a Seg fault error and crashes Sage entirely:

line 134: 4216 Segmentation fault (core dumped)

I believe it comes from the symbolic computation because the following code is working:

```sage: var('R k')
(R, k)
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41.), (R,0,6))
```

(Note the 41. instead of 41)

I have no idea where this comes from but it's really bad, because whatever is going wrong, crashing Sage is bad!

### comment:1 Changed 7 years ago by VivianePons

More funny stuff:

```sage: a = 319830986772877770815625*sqrt(pi)
sage: b = 1048576*gamma(41/2, 41/2*x)
sage: f = a - b
sage: plot(f, (R,0,1))  # works fine
sage: plot(f / sqrt(pi), (R,0,1))  # BOOM
```

### comment:2 Changed 7 years ago by nbruin

This is recent because on 6.5beta1 the error doesn't occur.

It looks like some infinite recursion in ginac. Excerpt from traceback obtained from `sage -gdb` after segv has occurred:

```...
#12573 0x00007fffc8185e3d in GiNaC::mul::eval(int) const ()
from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12574 0x00007fffc80f170b in GiNaC::ex::construct_from_basic(GiNaC::basic const&) () from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12575 0x00007fffc80acb59 in GiNaC::ex::ex(GiNaC::basic const&) ()
from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12576 0x00007fffc8185e3d in GiNaC::mul::eval(int) const ()
from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12577 0x00007fffc80f170b in GiNaC::ex::construct_from_basic(GiNaC::basic const&) () from /usr/local/sage/sage-git/local/lib/libpynac.so.1
#12578 0x00007fffc80acb59 in GiNaC::ex::ex(GiNaC::basic const&) ()
from /usr/local/sage/sage-git/local/lib/libpynac.so.1
...
```

### comment:4 Changed 7 years ago by buck

I think there's more to it than that even. Using the limited-precision numerics, I get this plot:

￼￼

In 2012, when the exact numerics still worked, I got this graph:

￼￼

### comment:5 Changed 7 years ago by kcrisman

• Type changed from PLEASE CHANGE to defect

### comment:6 follow-up: ↓ 9 Changed 7 years ago by rws

Interestingly, neither the first snippet in the description nor the BOOM snippet in comment:1 will crash on 6.7beta1/OpenSuSE. What machine/system is that, Viviane/buck/nbruin?

### comment:7 Changed 7 years ago by buck

I'm on OS X 10.10.3

### comment:8 Changed 7 years ago by VivianePons

Right, I just tried on 6.7 and indeed it is not crashing!

Buck I'd be happy to know if you obtain plot that "make sense" on 6.7. Sorry for the extra compiling, to get 6.7, just checkout the develop ranch and do a git pull. Let me know if you have problems.

### comment:9 in reply to: ↑ 6 Changed 7 years ago by nbruin

Interestingly, neither the first snippet in the description nor the BOOM snippet in comment:1 will crash on 6.7beta1/OpenSuSE. What machine/system is that, Viviane/buck/nbruin?

(on fedora) 6.7beta1 doesn't segfault anymore either. I haven't checked if the precision issues remain.

### comment:10 Changed 7 years ago by rws

No crash on OSX/6.7.beta1 either. Any crashes on 6.7beta1 at all?

### comment:11 follow-up: ↓ 12 Changed 7 years ago by buck

I also get no segfault in 6.7, but the numerical issue persists. Separate issue I suppose.

### comment:12 in reply to: ↑ 11 Changed 7 years ago by rws

I also get no segfault in 6.7, but the numerical issue persists. Separate issue I suppose.

Since no one can reproduce the crash on a recent Sage we will use this ticket for the numeric issue. I tried the code in the sage-devel posting but had no glitches in the graph. Can you please give a fully contained code snipped that shows the glitches on a recent Sage system?

### comment:13 Changed 7 years ago by buck

@rws: This is a reposting of my sage-support thread here where I gave some small reproductions of the problem.

The incomplete gamma function is resulting in negative outputs with positive real input, which shouldn't happen.

```sage: gamma(60, 30).numerical_approx()
-1.28306738270893e74
```

Attached is a notebook showing the same in graphical form.

If I feed in numbers of a higher precision, I get the correct result:

```sage: R10 = RealField(200)
sage: gamma(R10(60), R10(30)).numerical_approx()
1.38682990237888e80
```

My final problem is that there doesn't seem to be any way to get higher precision numbers to feed all the way through the plot function:

```sage: show(plot3d(log(gamma(k, k*x)), (x, R10(0), R10(5/2)), (k, R10(1), R10(101))))
Launched jmol viewer for Graphics3d Object
```

(you would see a 3d graph with a bite out of it where the precision issue exists)

I'm able to reproduce it today on the latest sage:

```\$ ./sage -v
SageMath Version 6.7.beta1, Release Date: 2015-04-15

\$ ./sage
sage: gamma(60, 30).numerical_approx()
-1.28306738270893e74
```

### Changed 7 years ago by buck

demo of the gamma function's instability

### comment:14 Changed 7 years ago by rws

• Description modified (diff)
• Summary changed from Symbolic computation of Gamma related function crashes Sage to numerical bug in incomplete gamma function

### comment:15 follow-up: ↓ 16 Changed 7 years ago by rws

• Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.

BTW this will be fixed with #16697 (needs review) because there we switched to mpmath for numerics. You can already use this in Sage via:

```sage: from sage.libs.mpmath import utils as mpmath_utils
sage: import mpmath
sage: mpmath_utils.call(mpmath.gammainc, 60, 0, 30)
1.28307801824120e74
```

This ticket may be closed therefore (after the Pari devs have ack'ed the bug) as duplicate of #16697 (augmenting its description with the bug).

### comment:16 in reply to: ↑ 15 Changed 7 years ago by rws

```sage: from sage.libs.mpmath import utils as mpmath_utils
sage: import mpmath
sage: mpmath_utils.call(mpmath.gammainc, 60, 0, 30)
1.28307801824120e74
```

Ah no, that is wrong as well. I'll report to Johansson too.

### comment:17 Changed 7 years ago by rws

• Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.

The upstream tickets are:

Thanks for the report!

### comment:18 Changed 7 years ago by rws

No, mpmath is correct. I confused lower/upper. So, as said, use the #16697 branch or better, review it:

```sage: gamma(60,30).n()
1.38682990237888e80
```

### comment:19 Changed 6 years ago by rws

• Report Upstream changed from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.

I was wrong referring to #16697 because this is a different bug. This issue is still open upstream.

### comment:20 Changed 6 years ago by buck

rws, could you please copy the relevant bits from #16697?

I'm going to work to close this in coming weeks.

### comment:21 Changed 6 years ago by buck

Discusion migrated from #16697

by buck:

... this patch doesn't fix the numerical instability in incomplete gamma as promised elsewhere.

```sage: gamma_inc_lower(25, 14.5).n()
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='mpmath')
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='pari')
-6.63538851954338e22
```

by buck:

Of note, gamma (even incomplete) is strictly positive in the reals; it's an integral of a strictly positive function.

by buck:

Wolfram gives the value of the same as 4.7e21, which seems correct. ​http://www.wolframalpha.com/input/?i=Gamma%5B25%2C+0%2C+14.5%5D

by buck:

Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken?

```sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21
```

by buck:

Same story if I ask mpmath directly. I can't imagine what's wrong.

```sage: call(mpmath.gammainc, 25, 0, 14.5)
4.71197173824555e21
```

by rws:

Same story if I ask mpmath directly. I can't imagine what's wrong.

'algorithm' isn't propagated to evalf so we always get the buggy pari values. As we would have to switch the default anyway as soon as a later pari version fixes the bug, I'll now simply disable pari. We'll have to live with mpmath's 53 bits for now.

by rws:

Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken?

```sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21
```

This works because it uses the gp interface (the libpari might be broken on our side). Ah, that's two different bugs. However, even if it is in our interface to libpari and we fix this, or if we use the gp interface, there is still the unfixed pari bug at ​http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1689PARI/GP

by rws:

Actually the 53-bit restriction applies only to larger values, so all is not lost.

### comment:22 follow-up: ↓ 23 Changed 6 years ago by buck

@rws: Could you expand on the "two different bugs" and link them (and/or create tickets for them)?

Also, could you explain why "because it uses the gp interface" would obscure the negative-number return value?

### comment:23 in reply to: ↑ 22 Changed 6 years ago by rws

@rws: Could you expand on the "two different bugs" and link them (and/or create tickets for them)?

Okay, this ticket is the `(60,30)` bug from comment:13, i.e. `gamma(60, 30).numerical_approx()`. The other one is the fact that `gamma(25,14.5)` and `CC(25).gamma_inc(14.5)`, which both use libpari, give

```sage: CC(25).gamma_inc(14.5)
6.86802286928673e23
sage: gamma(25,14.5)
6.86802286928673e23
```

while

```sage: pari('incgam(25, 14.5)')
6.15736429994994 E23
```

```sage: pari('incgam(60,30)')
-1.23084064495469 E74
sage: gamma(60,30).n()
-1.28306738270893e74
```

so there is a mismatch between the two interfaces, and the gp interface is sometimes right where libpari is wrong.

Also, could you explain why "because it uses the gp interface" would obscure the negative-number return value?

Because the gp interface is sometimes right, apparently.

However, I don't think it useful to open two tickets as long as we don't know if there really are two different bugs.

### comment:24 Changed 6 years ago by buck

The bug reported in the OP seems gone; I can't reproduce the segfault. Can we change the topic of the bug to strictly about the numerical issue? That seems to be what most of the discussion was about.

### comment:25 Changed 6 years ago by buck

@rws: I'll spend some time digging into this myself, and produce a pari and/or gp patch if I can, but I'm a bit lost as to where to start. Could you give some clues as to how I might pointpoint the issue, what code I should look at patching?

### comment:26 Changed 6 years ago by rws

In my opinion the Pari bug (which shows itself when using the GP/Pari program without Sage) should be fixed first. This means looking into the Pari code and reporting to its developers on that ticket already given in comment:17. After that it remains to be seen what's still wrong. For completeness Sage's libpari interface is in `libs/pari`.

### comment:27 Changed 6 years ago by rws

Karim Belabas of Pari just sent:

```I believe the problem is now fixed in master by the following commit
(rewritten from Henri Cohen's proposed patch)

commit 4276cc667123a92373c4aff7195eb728dd39f6e1
Author: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Date:   Fri Jan 15 15:58:17 2016 +0100

HC140- incgam(30,60) < 0. More generally, wrong results for s >> 1 [#1689]

Sorry it took so long. Please test and report problems !
```

### comment:28 Changed 6 years ago by buck

Thanks! Building master sagemath branch now...

### comment:29 Changed 6 years ago by buck

I still get the same broken graph in sage master.

Sage is at a1b60f2c3e06171caa4aa044a73ea58236686d30 and pari is at 'GP/PARI CALCULATOR Version 2.8.0 (development git-6157df4)'.

I can't tell immediately if that's before or after the referenced pari commit above.

### comment:30 Changed 6 years ago by buck

Ah after cloning pari, it's apparent that the installed sage version is four months behind the fixing commit.

I'll try to find in the docs how to get sage to install this newer version.

### comment:31 Changed 6 years ago by buck

Procedure:

```cd \$PARI_ROOT
TAG=\$(git describe)
git archive HEAD -o \$TAG.tar.gz --prefix=\$TAG/
mv pari-2.8-2230-g450ce38.tar.gz \$SAGE_ROOT/upstream/
cd \$SAGE_ROOT
./sage --fix-pkg-checksums
./sage -i -f -c pari
```

Changes to sage: (I had to touch up one of the patches to get it to apply.)

```diff --git a/build/pkgs/pari/checksums.ini b/build/pkgs/pari/checksums.ini
index c62c530..aabe08d 100644
--- a/build/pkgs/pari/checksums.ini
+++ b/build/pkgs/pari/checksums.ini
@@ -1,4 +1,4 @@
tarball=pari-VERSION.tar.gz
-sha1=fa23e0c8b6e38a356048d19224dc9b9658d77724
diff --git a/build/pkgs/pari/checksums.ini b/build/pkgs/pari/checksums.ini
index c62c530..aabe08d 100644
--- a/build/pkgs/pari/checksums.ini
+++ b/build/pkgs/pari/checksums.ini
@@ -1,4 +1,4 @@
tarball=pari-VERSION.tar.gz
-sha1=fa23e0c8b6e38a356048d19224dc9b9658d77724
-cksum=1220765312
+sha1=55299bfe042da491648897e830030083809d9967
+md5=97738f8e92bba498f7dfd723c8e9d910
+cksum=1221580104
diff --git a/build/pkgs/pari/package-version.txt b/build/pkgs/pari/package-version.txt
index 2b25bd1..5fdd71e 100644
--- a/build/pkgs/pari/package-version.txt
+++ b/build/pkgs/pari/package-version.txt
@@ -1 +1 @@
-2.8-1813-g6157df4.p0
+2.8-2230-g450ce38
diff --git a/build/pkgs/pari/patches/public_memory_functions.patch b/build/pkgs/pari/patches/public_memory_functions.patch
index b3726d7..ee14fa4 100644
--- a/build/pkgs/pari/patches/public_memory_functions.patch
+++ b/build/pkgs/pari/patches/public_memory_functions.patch
@@ -9,9 +9,9 @@ index 7067183..4ede6ed 100644
+void *  pari_mainstack_malloc(size_t size);
+void    pari_mainstack_mfree(void *s, size_t size);
+void    pari_mainstack_free(struct pari_mainstack *st);
- void    paristack_alloc(size_t rsize, size_t vsize);
void    paristack_newrsize(ulong newsize);
void    paristack_resize(ulong newsize);
+ void    paristack_setsize(size_t rsize, size_t vsize);
diff --git a/src/language/init.c b/src/language/init.c
index 7b5922d..2a578d7 100644
--- a/src/language/init.c
```

I also had to install bison to get the build to run, which I hadn't installed before, but I *had* built pari before (I think!) so I believe that's a new build dependency there?

Result: `Successfully installed pari-2.8-2230-g450ce38; The PARI self-tests all passed`

The sage console fails to start however, with `ImportError: sage/local/lib/python2.7/site-packages/sage/libs/pari/gen.so: undefined symbol: gp_algcenter`

Version 0, edited 6 years ago by buck (next)

### comment:32 Changed 6 years ago by buck

Subsequent `make` fails during pari autoregen:

```python -c "from sage_setup.autogen.pari import rebuild; rebuild()"
Generating PARI functions: (!_) (#_) (%) (%#) (+_) (-_) Catalan Col Colrev (DEBUGLEVEL) Euler I List Map Mat Mod (O)
...
galoispermtopol galoissubcyclo galoissubfields galoissubgroups gamma gammah gammamellininvTraceback (most recent call last):
File "<string>", line 1, in <module>
File "sage_setup/autogen/pari/__init__.py", line 5, in rebuild
G()
File "sage_setup/autogen/pari/generator.py", line 247, in __call__
self.handle_pari_function(**v)
File "sage_setup/autogen/pari/generator.py", line 173, in handle_pari_function
args, ret = parse_prototype(prototype, help)
File "sage_setup/autogen/pari/parser.py", line 204, in parse_prototype
raise ValueError('unknown prototype character %r' % c)
ValueError: unknown prototype character 'b'
```

I note that the build gets into a weird state afterward because the file being generated (`src/sage/libs/pari/auto_gen.pxi`) is written to even though its generation failed. We could avoid this state by making a patch like so:

• ## src/sage_setup/autogen/pari/generator.py

```diff --git a/src/sage_setup/autogen/pari/generator.py b/src/sage_setup/autogen/pari/generator.py
index 7aa9990..9b78d88 100644```
 a class PariFunctionGenerator(object): D = sorted(D.values(), key=lambda d: d['function']) sys.stdout.write("Generating PARI functions:") self.gen_file = open(self.gen_filename, 'w') self.gen_file = open(self.gen_filename + '.tmp', 'w') self.gen_file.write(gen_banner) self.instance_file = open(self.instance_filename, 'w') self.instance_file = open(self.instance_filename + '.tmp', 'w') self.instance_file.write(instance_banner) for v in D: class PariFunctionGenerator(object): self.gen_file.close() self.instance_file.close() # All done? Let's commit. os.rename(self.gen_filename + '.tmp', self.gen_filename) os.rename(self.instance_filename + '.tmp', self.instance_filename)

I'm out of my depth at this point. I think I need help from someone that actually know what pari is at all =X

### comment:33 Changed 6 years ago by buck

I believe I've fixed the `ValueError: unknown prototype character 'b'` issue. Pari recently added the concept of "bit precision", which (I believe) is the natural precision format in sage, so needs no conversion:

• ## src/sage_setup/autogen/pari/args.py

```diff --git a/src/sage_setup/autogen/pari/args.py b/src/sage_setup/autogen/pari/args.py
index 57356b4..a2749df 100644```
 a class PariArgumentPrec(PariArgumentClass): s = "        {name} = prec_bits_to_words({name})\n" return s.format(name=self.name) class PariArgumentBitPrec(PariArgumentClass): def _typerepr(self): return "bitprec" def ctype(self): return "long" def always_default(self): return "0" def get_argument_name(self, namesiter): return "bitprecision" class PariArgumentSeriesPrec(PariArgumentClass): def _typerepr(self): return "serprec" pari_arg_types = { 'U': PariArgumentULong, 'n': PariArgumentVariable, 'p': PariArgumentPrec, 'b': PariArgumentBitPrec, 'P': PariArgumentSeriesPrec, # Codes which are known but not actually supported for Sage

There's again a further issue that I don't understand yet:

```Compiling sage/libs/pari/pari_instance.pyx because it depends on sage/libs/pari/auto_instance.pxi.
[1/2] Cythonizing sage/libs/pari/gen.pyx

Error compiling Cython file:
------------------------------------------------------------
...

[0 0 -2 -1]
"""
cdef GEN _al = al.g
pari_catch_sig_on()
cdef GEN _ret = algmultable(_al)
^
------------------------------------------------------------

sage/libs/pari/auto_gen.pxi:1558:35: Call with wrong number of arguments (expected 2, got 1)
```

On brief inspection, I don't see that the number of arguments should have changed.

### comment:34 Changed 6 years ago by rws

While this is great work I believe you should open a separate Pari upgrade ticket for this.

### comment:35 Changed 6 years ago by buck

@rws Yes I didn't anticipate this would be such work. Will do.

### comment:36 Changed 6 years ago by buck

Pari upgrade work migrated to #19905.

### comment:37 follow-up: ↓ 38 Changed 6 years ago by buck

I was finally able to integrate pari master-branch with sage, and while the problem is much improved, it's not fixed. Here's the graph from the original bug report, using newest pari:

￼￼

Those discontinuities are false. An example bad value:

```sage: gamma(100, 7.01)
-2.71843697211100e151
```

I will report this upstream, and attach a worksheet showing more detail.

Last edited 6 years ago by buck (previous) (diff)

### comment:38 in reply to: ↑ 37 ; follow-up: ↓ 39 Changed 5 years ago by slelievre

This ticket seems to be solved now in Sage 7.3.beta9.

Trying the code from the ticket description...

```\$ ./sage -v
SageMath version 7.3.beta9, Release Date: 2016-07-22
\$ ./sage -q
sage: gamma(60, 30).numerical_approx()
1.38682990237888e80
```
```sage: %gp

--> Switching to PARI/GP interpreter <--

pari: \p 50
realprecision = 57 significant digits (50 digits displayed)
pari: incgam(60,30)
1.3868299023788801161747839921239559320835799009004 E80
pari: \p 30
realprecision = 38 significant digits (30 digits displayed)
pari: incgam(60,30)
1.38682990237888011617478399212 E80
pari: <ctrl-D>

--> Exiting back to Sage <--

```
```sage: _ = var('R k')
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41), (R,0,6))
Launched png viewer for Graphics object consisting of 1 graphics primitive
sage: integrated = -(gamma(1/2*k, 1/2*R*k) - gamma(1/2*k))/gamma(1/2*k)
sage: plot(integrated.subs(k=41.), (R,0,6))
Launched png viewer for Graphics object consisting of 2 graphics primitives
```
```sage: from sage.libs.mpmath import utils as mpmath_utils
sage: import mpmath
sage: mpmath_utils.call(mpmath.gammainc, 60, 0, 30)
1.28307801824120e74
```

I was finally able to integrate pari master-branch with sage, and while the problem is much improved, it's not fixed. Here's the graph from the original bug report, using newest pari:

Trying that now seems to work for me, I get nice graphs with no spikes.

```sage: sum(
....:     plot(
....:         integrated(k=k),
....:         (x, 0, 2.5),
....:         ymax=1.6,
....:         color=color,
....:         legend_label='k=%i' % k,
....:         figsize=6,
....:         aspect_ratio=1.0,
....:     )
....:     for k, color in zip(
....:         [21, 22, 23, 26, 31, 41],
....:         ['blue', 'brown', 'red', 'green', 'black', 'yellow', 'orange'],
....:     )
....: )
Launched png viewer for Graphics object consisting of 6 graphics primitives
```

Those discontinuities are false. An example bad value:

```sage: gamma(100, 7.01)
-2.71843697211100e151
```

I will report this upstream, and attach a worksheet showing more detail.

What is the correct value? I get the following in Sage 7.3.beta9.

```sage: gamma(100, 7.01)
9.33262154439442e155
```

### comment:39 in reply to: ↑ 38 Changed 5 years ago by rws

• Milestone changed from sage-6.7 to sage-duplicate/invalid/wontfix
• Report Upstream changed from Reported upstream. Developers acknowledge bug. to Fixed upstream, in a later stable release.
• Status changed from new to needs_review

This ticket seems to be solved now in Sage 7.3.beta9.

Agree.

What is the correct value? I get the following in Sage 7.3.beta9.

```sage: gamma(100, 7.01)
9.33262154439442e155
```
```sage: ComplexBallField(200)(100).gamma(ComplexBallField(200)(701/100))
[9.3326215443944152681699238856266700490715968264381621468593e+155 +/- 4.56e+96]
```

### comment:40 Changed 5 years ago by slelievre

• Status changed from needs_review to positive_review

### comment:41 Changed 5 years ago by vbraun

• Resolution set to invalid
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.