Opened 7 years ago
Closed 5 years ago
#18210 closed defect (invalid)
numerical bug in incomplete gamma function
Reported by:  VivianePons  Owned by:  

Priority:  major  Milestone:  sageduplicate/invalid/wontfix 
Component:  symbolics  Keywords:  sd67 
Cc:  Merged in:  
Authors:  Reviewers:  
Report Upstream:  Fixed upstream, in a later stable release.  Work issues:  
Branch:  Commit:  
Dependencies:  Stopgaps: 
Description (last modified by )
Wrong values from incomplete gamma:
$ ./sage v SageMath Version 6.7.beta1, Release Date: 20150415 $ ./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!
Attachments (6)
Change History (47)
comment:1 Changed 7 years ago by
comment:2 Changed 7 years ago by
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/sagegit/local/lib/libpynac.so.1 #12574 0x00007fffc80f170b in GiNaC::ex::construct_from_basic(GiNaC::basic const&) () from /usr/local/sage/sagegit/local/lib/libpynac.so.1 #12575 0x00007fffc80acb59 in GiNaC::ex::ex(GiNaC::basic const&) () from /usr/local/sage/sagegit/local/lib/libpynac.so.1 #12576 0x00007fffc8185e3d in GiNaC::mul::eval(int) const () from /usr/local/sage/sagegit/local/lib/libpynac.so.1 #12577 0x00007fffc80f170b in GiNaC::ex::construct_from_basic(GiNaC::basic const&) () from /usr/local/sage/sagegit/local/lib/libpynac.so.1 #12578 0x00007fffc80acb59 in GiNaC::ex::ex(GiNaC::basic const&) () from /usr/local/sage/sagegit/local/lib/libpynac.so.1 ...
comment:3 Changed 7 years ago by
The original discussion is here: https://groups.google.com/d/msg/sagesupport/7ex6McLE7AA/mvZ5RqD1FAJ
Changed 7 years ago by
Changed 7 years ago by
comment:4 Changed 7 years ago by
comment:5 Changed 7 years ago by
 Type changed from PLEASE CHANGE to defect
comment:6 followup: ↓ 9 Changed 7 years ago by
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
I'm on OS X 10.10.3
comment:8 Changed 7 years ago by
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
comment:10 Changed 7 years ago by
No crash on OSX/6.7.beta1 either. Any crashes on 6.7beta1 at all?
comment:11 followup: ↓ 12 Changed 7 years ago by
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
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 sagedevel 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
@rws: This is a reposting of my sagesupport 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: 20150415 $ ./sage sage: gamma(60, 30).numerical_approx() 1.28306738270893e74
comment:14 Changed 7 years ago by
 Description modified (diff)
 Summary changed from Symbolic computation of Gamma related function crashes Sage to numerical bug in incomplete gamma function
comment:15 followup: ↓ 16 Changed 7 years ago by
 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
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.
comment:17 Changed 7 years ago by
 Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.
The upstream tickets are:
 Pari/GP: http://pari.math.ubordeaux.fr/cgibin/bugreport.cgi?bug=1689PARI/GP
 mpmath: https://github.com/fredrikjohansson/mpmath/issues/311
Thanks for the report!
comment:18 Changed 7 years ago by
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
 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
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
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.ubordeaux.fr/cgibin/bugreport.cgi?bug=1689PARI/GP
by rws:
Actually the 53bit restriction applies only to larger values, so all is not lost.
comment:22 followup: ↓ 23 Changed 6 years ago by
@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 negativenumber return value?
comment:23 in reply to: ↑ 22 Changed 6 years ago by
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 negativenumber 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
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
@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
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
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.ubordeaux1.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
Thanks! Building master sagemath branch now...
comment:29 Changed 6 years ago by
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 git6157df4)'.
I can't tell immediately if that's before or after the referenced pari commit above.
comment:30 Changed 6 years ago by
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
Procedure:
cd $PARI_ROOT TAG=$(git describe) git archive HEAD o $TAG.tar.gz prefix=$TAG/ mv pari2.82230g450ce38.tar.gz $SAGE_ROOT/upstream/ cd $SAGE_ROOT ./sage fixpkgchecksums ./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
diff git a/build/pkgs/pari/checksums.ini b/build/pkgs/pari/checksums.ini index c62c530..aabe08d 100644
a b 1 1 tarball=pariVERSION.tar.gz 2 sha1= fa23e0c8b6e38a356048d19224dc9b9658d777243 md5= c753faaa4780de5ad8d461f0ffd70ecf4 cksum=122 07653122 sha1=55299bfe042da491648897e830030083809d9967 3 md5=97738f8e92bba498f7dfd723c8e9d910 4 cksum=1221580104 
build/pkgs/pari/packageversion.txt
diff git a/build/pkgs/pari/packageversion.txt b/build/pkgs/pari/packageversion.txt index 2b25bd1..5fdd71e 100644
a b 1 2.8 1813g6157df4.p01 2.82230g450ce38 
build/pkgs/pari/patches/public_memory_functions.patch
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 b index 7067183..4ede6ed 100644 9 9 +void * pari_mainstack_malloc(size_t size); 10 10 +void pari_mainstack_mfree(void *s, size_t size); 11 11 +void pari_mainstack_free(struct pari_mainstack *st); 12 void paristack_alloc(size_t rsize, size_t vsize);13 12 void paristack_newrsize(ulong newsize); 14 13 void paristack_resize(ulong newsize); 14 void paristack_setsize(size_t rsize, size_t vsize); 15 15 diff git a/src/language/init.c b/src/language/init.c 16 16 index 7b5922d..2a578d7 100644 17 17  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 pari2.82230g450ce38; The PARI selftests all passed
The sage console fails to start however, with ImportError: sage/local/lib/python2.7/sitepackages/sage/libs/pari/gen.so: undefined symbol: gp_algcenter
comment:32 Changed 6 years ago by
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 b class PariFunctionGenerator(object): 233 233 D = sorted(D.values(), key=lambda d: d['function']) 234 234 sys.stdout.write("Generating PARI functions:") 235 235 236 self.gen_file = open(self.gen_filename , 'w')236 self.gen_file = open(self.gen_filename + '.tmp', 'w') 237 237 self.gen_file.write(gen_banner) 238 self.instance_file = open(self.instance_filename , 'w')238 self.instance_file = open(self.instance_filename + '.tmp', 'w') 239 239 self.instance_file.write(instance_banner) 240 240 241 241 for v in D: … … class PariFunctionGenerator(object): 249 249 250 250 self.gen_file.close() 251 251 self.instance_file.close() 252 253 # All done? Let's commit. 254 os.rename(self.gen_filename + '.tmp', self.gen_filename) 255 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
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 b class PariArgumentPrec(PariArgumentClass): 254 254 s = " {name} = prec_bits_to_words({name})\n" 255 255 return s.format(name=self.name) 256 256 257 class PariArgumentBitPrec(PariArgumentClass): 258 def _typerepr(self): 259 return "bitprec" 260 def ctype(self): 261 return "long" 262 def always_default(self): 263 return "0" 264 def get_argument_name(self, namesiter): 265 return "bitprecision" 266 257 267 class PariArgumentSeriesPrec(PariArgumentClass): 258 268 def _typerepr(self): 259 269 return "serprec" … … pari_arg_types = { 277 287 'U': PariArgumentULong, 278 288 'n': PariArgumentVariable, 279 289 'p': PariArgumentPrec, 290 'b': PariArgumentBitPrec, 280 291 'P': PariArgumentSeriesPrec, 281 292 282 293 # 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
While this is great work I believe you should open a separate Pari upgrade ticket for this.
comment:35 Changed 6 years ago by
@rws Yes I didn't anticipate this would be such work. Will do.
comment:36 Changed 6 years ago by
Pari upgrade work migrated to #19905.
comment:37 followup: ↓ 38 Changed 6 years ago by
I was finally able to integrate pari masterbranch 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.
Changed 6 years ago by
Changed 6 years ago by
Changed 6 years ago by
comment:38 in reply to: ↑ 37 ; followup: ↓ 39 Changed 5 years ago by
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: 20160722 $ ./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: <ctrlD> > 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 masterbranch 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.71843697211100e151I 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
 Milestone changed from sage6.7 to sageduplicate/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
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]
comment:40 Changed 5 years ago by
 Status changed from needs_review to positive_review
comment:41 Changed 5 years ago by
 Resolution set to invalid
 Status changed from positive_review to closed
More funny stuff: