numerical evaluation of `erf` crashes PARI at large numbers

[benjaminfjones]

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

Upstream: Reported upstream. Developers deny it's a bug.

CC: @sagetrac-dsm @kcrisman

Component: symbolics

Keywords: erf pari sd40.5

Reviewer: Benjamin Jones

Issue created by migration from https://trac.sagemath.org/ticket/13003

[sagetrac-dsm]

comment:2

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..

[jdemeyer]

comment:3

This is an upstream PARI bug.

[jdemeyer]

Upstream: Reported upstream. Little or no feedback.

[jdemeyer]

Description changed:

--- 
+++ 
@@ -36,3 +36,5 @@
 
 PariError:  (15)
 ```
+
+Reported upstream: [http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1320](http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1320)

[jdemeyer]

Changed upstream from Reported upstream. Little or no feedback. to Reported upstream. Developers deny it's a bug.

[benjaminfjones]

comment:6

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.

[sagetrac-dsm]

comment:7

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

[jdemeyer]

comment:8

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

[jdemeyer]

Reviewer: Benjamin Jones