numerical bug in incomplete gamma function #18210
Comments
|
comment:1
More funny stuff: |
|
comment:2
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 |
|
comment:3
The original discussion is here: https://groups.google.com/d/msg/sage-support/7ex6McLE7AA/mvZ5Rq-D1FAJ |
|
Attachment: working-graph-2012.png Attachment: broken-graph-2015.png |
|
comment:4
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:6
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
I'm on OS X 10.10.3 |
|
comment:8
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
Replying to @rwst:
(on fedora) 6.7beta1 doesn't segfault anymore either. I haven't checked if the precision issues remain. |
|
comment:10
No crash on OSX/6.7.beta1 either. Any crashes on 6.7beta1 at all? |
|
comment:11
I also get no segfault in 6.7, but the numerical issue persists. |
|
comment:12
Replying to @sagetrac-buck:
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
@rwst: 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. Attached is a notebook showing the same in graphical form. If I feed in numbers of a higher precision, I get the correct result: 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: (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: |
|
Attachment: unstable incomplete gamma.ipynb.json.gz demo of the gamma function's instability |
This comment has been minimized.
This comment has been minimized.
rwst
changed the title
|
comment:15
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: 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). |
|
Upstream: Not yet reported upstream; Will do shortly. |
|
comment:16
Replying to @rwst:
Ah no, that is wrong as well. I'll report to Johansson too. |
|
Changed upstream from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet. |
|
comment:17
The upstream tickets are:
Thanks for the report! |
|
comment:18
No, mpmath is correct. I confused lower/upper. So, as said, use the #16697 branch or better, review it: |
|
comment:19
I was wrong referring to #16697 because this is a different bug. This issue is still open upstream. |
|
Changed upstream from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug. |
|
comment:20
rws, could you please copy the relevant bits from #16697? I'm going to work to close this in coming weeks. |
|
comment:21
Discusion migrated from #16697 by buck: ... this patch doesn't fix the numerical instability in incomplete gamma as promised elsewhere. 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. by buck: Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken? by buck: Same story if I ask mpmath directly. I can't imagine what's wrong. by rws: Replying to buck:
'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:
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 by rws: Actually the 53-bit restriction applies only to larger values, so all is not lost. |
|
comment:22
@rwst: 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
Replying to @sagetrac-buck:
Okay, this ticket is the while But then, we see also so there is a mismatch between the two interfaces, and the gp interface is sometimes right where libpari is wrong.
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
The bug reported in the OP seems gone; I can't reproduce the segfault. |
|
comment:25
@rwst: 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
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 |
|
comment:27
Karim Belabas of Pari just sent: |
|
comment:28
Thanks! Building master sagemath branch now... |
|
comment:29
I still get the same broken graph in sage master. Sage is at a1b60f2 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
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
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 pariChanges 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
-md5=c753faaa4780de5ad8d461f0ffd70ecf
-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.cI 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: The sage console fails to start however, with |
|
comment:32
Subsequent I note that the build gets into a weird state afterward because the file being generated ( diff --git a/src/sage_setup/autogen/pari/generator.py b/src/sage_setup/autogen/pari/generator.py
index 7aa9990..9b78d88 100644
--- a/src/sage_setup/autogen/pari/generator.py
+++ b/src/sage_setup/autogen/pari/generator.py
@@ -233,9 +233,9 @@ 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:
@@ -249,3 +249,7 @@ 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
I believe I've fixed the diff --git a/src/sage_setup/autogen/pari/args.py b/src/sage_setup/autogen/pari/args.py
index 57356b4..a2749df 100644
--- a/src/sage_setup/autogen/pari/args.py
+++ b/src/sage_setup/autogen/pari/args.py
@@ -254,6 +254,16 @@ 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"
@@ -277,6 +287,7 @@ pari_arg_types = {
'U': PariArgumentULong,
'n': PariArgumentVariable,
'p': PariArgumentPrec,
+ 'b': PariArgumentBitPrec,
'P': PariArgumentSeriesPrec,
# Codes which are known but not actually supported for SageThere's again a further issue that I don't understand yet: On brief inspection, I don't see that the number of arguments should have changed. |
|
comment:34
While this is great work I believe you should open a separate Pari upgrade ticket for this. |
|
comment:35
@rwst Yes I didn't anticipate this would be such work. Will do. |
|
comment:36
Pari upgrade work migrated to #19905. |
|
comment:37
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: I will report this upstream, and attach a worksheet showing more detail. |
|
Attachment: unstable incomplete gamma.pdf.gz |
|
Attachment: unstable incomplete gamma.ipynb.gz Attachment: broken-graph-2016.png |
|
comment:38
This ticket seems to be solved now in Sage 7.3.beta9. Trying the code from the ticket description... Replying to @sagetrac-buck:
Trying that now seems to work for me, I get nice graphs with no spikes.
What is the correct value? I get the following in Sage 7.3.beta9. |
|
comment:39
Replying to @slel:
Agree.
|
|
Changed upstream from Reported upstream. Developers acknowledge bug. to Fixed upstream, in a later stable release. |
Wrong values from incomplete gamma:
Reason is a bug in Pari:
Previous description:
The following code:
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:
(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!
Upstream: Fixed upstream, in a later stable release.
Component: symbolics
Keywords: sd67
Issue created by migration from https://trac.sagemath.org/ticket/18210
The text was updated successfully, but these errors were encountered: