Opened 7 years ago

Closed 7 years ago

#13003 closed defect (wontfix)

numerical evaluation of `erf` crashes PARI at large numbers

Reported by: benjaminfjones Owned by: burcin
Priority: minor Milestone: sage-duplicate/invalid/wontfix
Component: symbolics Keywords: erf pari sd40.5
Cc: dsm, kcrisman Merged in:
Authors: Reviewers: Benjamin Jones
Report Upstream: Reported upstream. Developers deny it's a bug. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by jdemeyer)

In sage-5.0:

sage: for i in range(10):
....:     print erf(45*10**i).n()
....:     
1.00000000000000
1.00000000000000
1.00000000000000
-infinity
1.00000000000000
1.00000000000000
-1.89139086613300e232327576
1.00000000000000
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (9099, 0))

---------------------------------------------------------------------------
PariError                                 Traceback (most recent call last)

/Users/jonesbe/sage/sage-5.0/<ipython console> in <module>()

/Users/jonesbe/sage/sage-5.0/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:18303)()

/Users/jonesbe/sage/sage-5.0/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._convert (sage/symbolic/expression.cpp:5394)()

/Users/jonesbe/sage/sage-5.0/local/lib/python2.7/site-packages/sage/functions/other.pyc in _evalf_(self, x, parent)
    221         except AttributeError: # not a Sage parent
    222             prec = 0
--> 223         return parent(1) - parent(pari(x).erfc(prec))
    224 
    225     def _derivative_(self, x, diff_param=None):

/Users/jonesbe/sage/sage-5.0/local/lib/python2.7/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:49843)()

PariError:  (15)

Reported upstream: http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1320

Change History (9)

comment:1 Changed 7 years ago by benjaminfjones

  • Component changed from PLEASE CHANGE to symbolics
  • Owner changed from tbd to burcin

comment:2 Changed 7 years ago by dsm

FWIW (and for the record) I get slightly different behaviour on an ubuntu 12.04 vbox running on a core duo macbook with 5.1.beta0:

sage: for i in range(10):
....:     try:   
....:         print erf(45*10**i).n()
....:     except:
....:         print 'problem'
....:         
1.00000000000000
1.00000000000000
1.00000000000000
problem
problem
problem
problem
problem
problem
problem

where the problem is the same Pari (15) error in the ticket. Pari should certainly be returning 1 here, although any error is certainly better than -1.89139086613300e232327576..

comment:3 Changed 7 years ago by jdemeyer

This is an upstream PARI bug.

comment:4 Changed 7 years ago by jdemeyer

  • Description modified (diff)
  • Report Upstream changed from N/A to Reported upstream. Little or no feedback.

comment:5 Changed 7 years ago by jdemeyer

  • Report Upstream changed from Reported upstream. Little or no feedback. to Reported upstream. Developers deny it's a bug.

comment:6 Changed 7 years ago by benjaminfjones

Here is an example that is closer to (part of) the problem:

Consider erf(10^i) where i=5,6,7:

sage: pari(10**5).erfc()
5.23488067975405 E-4342944825
sage: RR(pari(10**5).erfc())
0.000000000000000

looks good!

sage: pari(10**6).erfc()
3.15934761259943 E-434294481910
sage: RR(pari(10**6).erfc())
5.64243263573681e124617551

there's a bug.

Then,

sage: pari(10**7).erfc()
3.70390595226759 E-43429448190333
sage: RR(pari(10**7).erfc())
0.000000000000000

looks good again ?!?

If the input is too large, there is an overflow in PARI (the exponent is too large for PARI to represent):

sage: pari(10**10).erfc()
---------------------------------------------------------------------------
PariError                                 Traceback (most recent call last)

/Users/jonesbe/sage/sage-5.0/<ipython console> in <module>()

/Users/jonesbe/sage/sage-5.0/local/lib/python2.7/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:49843)()

PariError:  (15)

and that's something we should probably address separately.

comment:7 Changed 7 years ago by dsm

Might be worth noting that we get a more reasonable result by building a number from the string:

sage: pari(10**6).erfc()
3.15934761259943 E-434294481910
sage: RR(pari(10**6).erfc())
5.64243263573681e124617551
sage: RR(str(pari(10**6).erfc()))
0.000000000000000

comment:8 Changed 7 years ago by jdemeyer

The conversion PARI->RR is a different issue: #13033.

comment:9 Changed 7 years ago by jdemeyer

  • Milestone changed from sage-5.3 to sage-duplicate/invalid/wontfix
  • Resolution set to wontfix
  • Reviewers set to Benjamin Jones
  • Status changed from new to closed
Note: See TracTickets for help on using tickets.