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:

Status badges

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)

trac_5739.patch (1.1 KB) - added by Mike Hansen 13 years ago.
5739-CDF-zeta.patch (705 bytes) - added by Robert Bradshaw 13 years ago.
5739-CC-RR.patch (2.5 KB) - added by Robert Bradshaw 13 years ago.
apply this and 5739-CDF-zeta.patch
5739-complex-zeta.patch (1.3 KB) - added by Robert Bradshaw 13 years ago.
apply only this patch
5739-complex-zeta-with-header.patch (1.5 KB) - added by David Loeffler 12 years ago.
Version with Mercurial header

Download all attachments as: .zip

Change History (22)

comment:1 Changed 14 years ago by William Stein

Cc: Robert Bradshaw added

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:

        cdef pari_sp sp
        sp = avma
        _sig_on
        x = self._new_from_gen_c(  gzeta(self._gen(), PREC),   sp)
        _sig_off
        return x

somewhat fixes the problem, though doing

sage: zeta(CDF(0))

then raises a RuntimeError?. Unfortunately, doing this doesn't work because the _sig_on/_sig_off stuff doesn't play well with Cython exceptions:

        cdef pari_sp sp
        sp = avma
        try:
           _sig_on
           x = self._new_from_gen_c(  gzeta(self._gen(), PREC),   sp)
           _sig_off
        except:
           raise ValueError
        return x

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.

comment:2 Changed 14 years ago by Mike Hansen

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 Mike Hansen

Report Upstream: N/A
Status: newneeds_review

comment:4 Changed 13 years ago by Mike Hansen

Authors: Mike Hansen

Changed 13 years ago by Mike Hansen

Attachment: trac_5739.patch added

comment:5 Changed 13 years ago by Karl-Dieter Crisman

Status: needs_reviewneeds_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 Robert Bradshaw

Attachment: 5739-CDF-zeta.patch added

comment:6 Changed 13 years ago by Robert Bradshaw

Status: needs_infoneeds_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 Karl-Dieter Crisman

Reviewers: Karl-Dieter Crisman, Robert Bradshaw
Status: needs_reviewneeds_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 Fredrik Johansson

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.

Changed 13 years ago by Robert Bradshaw

Attachment: 5739-CC-RR.patch added

apply this and 5739-CDF-zeta.patch

comment:9 Changed 13 years ago by Robert Bradshaw

Status: needs_workneeds_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 Robert Bradshaw

Summary: Have zeta of 1 return value be consistent in different ringschanged from zeta(CDF(1)) go boom! + zeta of 1 return value be consistent in different rings

comment:11 Changed 13 years ago by David Loeffler

Status: needs_reviewneeds_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 Burcin Erocal

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.

Changed 13 years ago by Robert Bradshaw

Attachment: 5739-complex-zeta.patch added

apply only this patch

comment:13 Changed 13 years ago by Robert Bradshaw

Status: needs_workneeds_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 David Loeffler

Authors: Mike HansenMike Hansen, Robert Bradshaw
Reviewers: Karl-Dieter Crisman, Robert BradshawKarl-Dieter Crisman, Robert Bradshaw, David Loeffler
Status: needs_reviewpositive_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 Mitesh Patel

Status: positive_reviewneeds_work

The patch is missing a Mercurial header. Could someone add this and restore the positive review?

Changed 12 years ago by David Loeffler

Version with Mercurial header

comment:16 Changed 12 years ago by David Loeffler

Status: needs_workpositive_review

Done

comment:17 Changed 12 years ago by Mitesh Patel

Merged in: sage-4.6.alpha2
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.