numerical bug in incomplete gamma function

[VivianePons]

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!

[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

[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
...

[buck]

The original discussion is here: ​https://groups.google.com/d/msg/sage-support/7ex6McLE7AA/mvZ5Rq-D1FAJ

[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:



[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?

[buck]

I'm on OS X 10.10.3

[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.

[nbruin]

Replying to 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?

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

[rws]

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

[buck]

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

[rws]

Replying to buck:

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?

[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

[buck]

demo of the gamma function's instability

[rws]

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).

[rws]

Replying to 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.

[rws]

The upstream tickets are:

• Pari/GP: ​http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1689PARI/GP
• mpmath: ​https://github.com/fredrik-johansson/mpmath/issues/311

Thanks for the report!

[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

[rws]

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

[buck]

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

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

[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.

​​https://en.wikipedia.org/wiki/Incomplete_gamma_function#Definition

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:

Replying to buck:

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:

Replying to 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

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.

[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?

[rws]

Replying to buck:

@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

But then, we see also

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.

[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.

[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?

[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.

[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 !

[buck]

Thanks! Building master sagemath branch now...

[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.

[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.

[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.)

build/pkgs/pari/checksums.ini

@@ -1 +1 @@
 tarball=pari-VERSION.tar.gz
-sha1=fa23e0c8b6e38a356048d19224dc9b9658d77724
-md5=c753faaa4780de5ad8d461f0ffd70ecf
-cksum=1220765312
+sha1=55299bfe042da491648897e830030083809d9967
+md5=97738f8e92bba498f7dfd723c8e9d910
+cksum=1221580104

build/pkgs/pari/package-version.txt

@@ -1 @@
-2.8-1813-g6157df4.p0
+2.8-2230-g450ce38

build/pkgs/pari/patches/public_memory_functions.patch

@@ -9 +9 @@
 +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

[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

@@ -233 +233 @@
         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:
@@ -249 +249 @@
 
         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

[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

@@ -254 +254 @@
         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"
@@ -277 +287 @@
         '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.

[rws]

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

[buck]

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

[buck]

Pari upgrade work migrated to #19905.

[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.

[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

Replying to 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:

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

[rws]

Replying to slelievre:

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]