Opened 14 years ago
Closed 12 years ago
#5739 closed defect (fixed)
changed from zeta(CDF(1)) go boom! + zeta of 1 return value be consistent in different rings
Reported by: | William Stein | Owned by: | fredrik |
---|---|---|---|
Priority: | major | Milestone: | sage-4.6 |
Component: | number theory | Keywords: | |
Cc: | Robert Bradshaw, fredrik | Merged in: | sage-4.6.alpha2 |
Authors: | Mike Hansen, Robert Bradshaw | Reviewers: | Karl-Dieter Crisman, Robert Bradshaw, David Loeffler |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description
wstein@bsd:~/build/sage-3.4.1.rc1$ uname -a Darwin bsd.local 9.6.0 Darwin Kernel Version 9.6.0: Mon Nov 24 17:37:00 PST 2008; root:xnu-1228.9.59~1/RELEASE_I386 i386 wstein@bsd:~/build/sage-3.4.1.rc1$ sage ---------------------------------------------------------------------- | Sage Version 3.4.1.rc1, Release Date: 2009-04-05 | | Type notebook() for the GUI, and license() for information. | ---------------------------------------------------------------------- sage: zeta(CDF(1)) ------------------------------------------------------------ Unhandled SIGSEGV: A segmentation fault occured in SAGE. This probably occured because a *compiled* component of SAGE has a bug in it (typically accessing invalid memory) or is not properly wrapped with _sig_on, _sig_off. You might want to run SAGE under gdb with 'sage -gdb' to debug this. SAGE will now terminate (sorry). ------------------------------------------------------------
Attachments (5)
Change History (22)
comment:1 Changed 14 years ago by
Cc: | Robert Bradshaw added |
---|
comment:2 Changed 14 years ago by
Cc: | fredrik added |
---|---|
Owner: | changed from William Stein to fredrik |
We should be able to replace the calls to pari with calls to mpmath.
comment:3 Changed 13 years ago by
Report Upstream: | → N/A |
---|---|
Status: | new → needs_review |
comment:4 Changed 13 years ago by
Authors: | → Mike Hansen |
---|
Changed 13 years ago by
Attachment: | trac_5739.patch added |
---|
comment:5 Changed 13 years ago by
Status: | needs_review → needs_info |
---|
The patch looks fine, but results in zeta of a CDF being approximately fifty times slower. This seems problematic, and perhaps also something that would happen a lot if we start switching things to mpmath? Mpmath looks like a great package, but if it has the same issue as NetworkX versus C graphs, we might not want to move default behavior there quite yet.
Marking as needs_info since there does not seem to be a Sage-wide policy on mpmath at this point, and I am reluctant to give positive review to such a marked slowdown.
Changed 13 years ago by
Attachment: | 5739-CDF-zeta.patch added |
---|
comment:6 Changed 13 years ago by
Status: | needs_info → needs_review |
---|
I added an alternative patch that special cases the one pole at s=1 (returning the unsigned infinity, as gamma does).
comment:7 Changed 13 years ago by
Reviewers: | → Karl-Dieter Crisman, Robert Bradshaw |
---|---|
Status: | needs_review → needs_work |
Summary: | zeta(CDF(1)) go boom! → Have zeta of 1 return value be consistent in different rings |
I hate to send this back to the drawing board again, but let's just fix things once and for all...
sage: zeta(CDF(1)) Infinity sage: zeta(CC(1)) --------------------------------------------------------------------------- PariError Traceback (most recent call last) /Users/.../<ipython console> in <module>() /Users/.../sage-4.3.1.alpha2/local/lib/python2.6/site-packages/sage/functions/transcendental.pyc in zeta(s) 153 """ 154 try: --> 155 return s.zeta() 156 except AttributeError: 157 return ComplexField()(s).zeta() /Users/.../sage-4.3.1.alpha2/local/lib/python2.6/site-packages/sage/rings/complex_number.so in sage.rings.complex_number.ComplexNumber.zeta (sage/rings/complex_number.c:12174)() /Users/.../sage-4.3.1.alpha2/local/lib/python2.6/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:44110)() PariError: (8)
Can we think of any other places where this needs to be checked? For instance, zeta(1) returns this error too, though I think it inherits it from the CC example.
Also, regarding whether it should be Infinity or some signed infinity:
sage: zeta(RR(1)) +infinity sage: zeta(RDF(1)) +infinity
I'm not saying which one is better, just what the current behavior is. What do folks think?
comment:8 Changed 13 years ago by
kcrisman:
Starting with the next version, mpmath uses the Riemann-Siegel formula, so it should be much faster than Pari for large imaginary parts near the critical strip. Right now I even get a segmentation fault if I try to compute zeta(CDF(1/2+10000000*I)) in Sage.
For CDF, zeta could also use mpmath.fp.zeta that will be available in the next version of mpmath. This function is currently typically 10-50 times faster than mpmath.mp.zeta. However, fp.zeta loses accuracy proportional to the magnitude of the imaginary part near the critical strip, so the question is whether this loss would be acceptable. For small imaginary parts, it's quite accurate.
Both functions could be accelerated in Sage by overriding the base case of an internal function (mpmath/functions/zeta.py/_zetasum in the svn trunk, if anyone wants a go). This should require just few lines of Cython.
Other than that, I would recommend keeping Pari where it's faster.
comment:9 Changed 13 years ago by
Status: | needs_work → needs_review |
---|
I fixed CC. As to whether it should be a signed or unsigned infinity, I went with unsigned because it has a simple pole there.
10000.5772229475 sage: zeta(.9999) -9999.42279161783
When the new mpmath gets released, we could open a ticket with timings and accuracy comparison. Generally we favor correctness over speed.
comment:10 Changed 13 years ago by
Summary: | Have zeta of 1 return value be consistent in different rings → changed from zeta(CDF(1)) go boom! + zeta of 1 return value be consistent in different rings |
---|
comment:11 Changed 13 years ago by
Status: | needs_review → needs_work |
---|
Patch applies cleanly to 4.5.alpha1 and builds fine, but some doctests fail:
sage -t sage/rings/real_mpfr.pyx ********************************************************************** File "/storage/masiao/sage-4.5.alpha1/devel/sage-reviewing/sage/rings/real_mpfr.pyx", line 4487: sage: R(1).zeta() Expected: Infinity Got: +infinity ********************************************************************** 1 items had failures: 1 of 12 in __main__.example_149 ***Test Failed*** 1 failures. For whitespace errors, see the file /home/masiao/.sage//tmp/.doctest_real_mpfr.py [10.2 s] ********************************************************************** File "/storage/masiao/sage-4.5.alpha1/devel/sage-reviewing/sage/rings/complex_number.pyx", line 2093: sage: zeta(1) Expected: Infinity Got: zeta(1) ********************************************************************** 1 items had failures: 1 of 8 in __main__.example_72 ***Test Failed*** 1 failures. For whitespace errors, see the file /home/masiao/.sage//tmp/.doctest_complex_number.py [8.6 s]
Moreover, it doesn't seem to live up to the promise in the title of making the return value of zeta(1) consistent:
sage: zeta(RR(1)) +infinity sage: zeta(RDF(1)) Infinity
comment:12 Changed 13 years ago by
Since #8864, zeta(1)
returns the answer given by GiNaC, which leaves it unevaluated because it doesn't know about infinity. I'll change this in the next pynac release to return unsigned infinity.
comment:13 Changed 13 years ago by
Status: | needs_work → needs_review |
---|
It's really sad we still have this trivial-to-fix hard crash after over a year!
I've attached a patch that fixes the behavior of CDF(1).zeta() and CC(1).zeta(). It leaves the real fields alone, which I think is fine 'cause they have reasonable representations of (an) infinity, and we usually try to return something of the same type. (IN the complex case, infinity is a lot messier without a lot of special casing that's beyond the scope of this ticket...)
sage: RR(1).zeta(), RDF(1).zeta(), CC(1).zeta(), CDF(1).zeta() (+infinity, +infinity, Infinity, Infinity)
comment:14 Changed 12 years ago by
Authors: | Mike Hansen → Mike Hansen, Robert Bradshaw |
---|---|
Reviewers: | Karl-Dieter Crisman, Robert Bradshaw → Karl-Dieter Crisman, Robert Bradshaw, David Loeffler |
Status: | needs_review → positive_review |
OK, that looks fine. I'd still argue that it should be an unsigned infinity in the real field cases as well, but (as you say) more or less anything is better than the current behaviour! Let's get this in.
comment:15 Changed 12 years ago by
Status: | positive_review → needs_work |
---|
The patch is missing a Mercurial header. Could someone add this and restore the positive review?
Changed 12 years ago by
Attachment: | 5739-complex-zeta-with-header.patch added |
---|
Version with Mercurial header
comment:17 Changed 12 years ago by
Merged in: | → sage-4.6.alpha2 |
---|---|
Resolution: | → fixed |
Status: | positive_review → closed |
I think Robert Bradshaw massively optimized CDF special functions, and that's where this comes from. Changing the code for zeta in complex_double.pyx to:
somewhat fixes the problem, though doing
then raises a RuntimeError?. Unfortunately, doing this doesn't work because the _sig_on/_sig_off stuff doesn't play well with Cython exceptions:
So I'm not sure how to fix this in general in a nice way with the right exception, but definitely adding _sig_on/_sig_off's around all the calls to pari is a very very good idea.