Fix numeric evaluation of error function

[jdemeyer]

Because the argument is converted to float first, the error function erf() cannot numerically be evaluated for complex arguments:

sage: erf(pi - 1/2*I).n()
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', (1108, 0))

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

/usr/local/src/sage-4.7.2.alpha2/<ipython console> in <module>()

/usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:17793)()

/usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._convert (sage/symbolic/expression.cpp:5004)()

/usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/functions/other.pyc in _evalf_(self, x, parent)
     79         if prec > 53:
     80             raise NotImplementedError, "erf not implemented for precision higher than 53"
---> 81         return parent(1 - pari(float(x)).erfc())
     82
     83     def _derivative_(self, x, diff_param=None):

/usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/rings/complex_number.so in sage.rings.complex_number.ComplexNumber.__float__ (sage/rings/complex_number.c:7233)()

TypeError: Unable to convert 3.14159265358979 - 0.500000000000000*I to float; use abs() or real_part() as desired

Depends on #11130
Depends on #11321
Depends on #11854
Depends on #11891
Depends on #11890
Depends on #11836
Depends on #11952

Component: numerical

Keywords: erf erfc

Author: Jeroen Demeyer

Reviewer: Karl-Dieter Crisman

Merged: sage-4.8.alpha6

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

[kcrisman]

comment:1

#1173 is closely related, and #8983 is somewhat related. I haven't tried whether #1173 solves this issue or not.

[jdemeyer]

comment:2

#1173 and #11948 both seem to solve the problem is a different way (this ticket here using the new PARI and #1173 using mpmath).

[jdemeyer]

Author: Jeroen Demeyer

[jdemeyer]

Changed dependencies from #11130 to #11130, #11321, #11854, #11891, #11890, #11836, #11904, #11952

[kcrisman]

comment:6

Thanks for that info; I didn't realize it was the new Pari that made this possible.

I'd say that if someone reviews this, it should go in. #1173 is probably not going to get finished immediately, and we can always switch to mpmath if necessary or convenient or speedy at that time.

[jdemeyer]

comment:7

Replying to @kcrisman:

I'd say that if someone reviews this

this and #11130 and #11321 and #11854 and #11904 and #11952 unfortunately...

[jdemeyer]

Changed dependencies from #11130, #11321, #11854, #11891, #11890, #11836, #11904, #11952 to #11130, #11321, #11854, #11891, #11890, #11836, #11952

[jdemeyer]

Attachment: 11948.patch.gz

[kcrisman]

Apply after main patch

[kcrisman]

comment:9

Attachment: trac_11948-reviewer.patch.gz

Looks good. I've added a tiny bit of doc, because erf? gives something very short otherwise now, and fixed some grammar in the place where Jeroen did something I can't figure out with whitespace. Docs still look fine there, so all seems well.

Positive review, unless the author finds something he doesn't like about the review patch.

[kcrisman]

Reviewer: Karl-Dieter Crisman

[jdemeyer]

comment:10

Replying to @kcrisman:

Jeroen did something I can't figure out with whitespace.

Probably changing TABS to spaces.

[jdemeyer]

Merged: sage-4.8.alpha6

[kcrisman]

comment:14

Unfortunately, this causes a nasty problem.

This example in the documentation is fine - comparing mpmath:

sage: mpmath.erf(pi-1/2*i)
mpc(real='1.0000111669099367', imag='1.6332655417638451e-6')

But anything along the imaginary axis seems to be off by exactly 1:

sage: for z in [3,33,333,3333,33333]:
....:     mpmath.erf(i*z); erf(n(z)*i)
....:     
mpc(real='0.0', imag='1629.9946226015657')
1.00000000000000 + 1629.86732385786*I
mpc(real='0.0', imag='1.5128697751040891e+471')
1.00000000000000 + 1.51286977510409e471*I
mpc(real='0.0', imag='5.1260939089106243e+48155')
1.00000000000000 + 5.12609390891062e48155*I
mpc(real='0.0', imag='2.6385510598470926e+4824525')
1.00000000000000 + 2.63855105984709e4824525*I

But other values seem ok.

sage: for z in [3,33,333,3333]:
    mpmath.erf(1+i*z); erf(1.+n(z)*i)
....:     
mpc(real='-330.81538696857206', imag='443.38888183939281')
-330.815386892947 + 443.388881909712*I
mpc(real='2.0957487368415288e+468', imag='-5.5629367605580166e+470')
2.09574873684153e468 - 5.56293676055802e470*I
mpc(real='-3.8930178706420656e+48153', imag='1.8853741770265906e+48155')
-3.89301787064206e48153 + 1.88537417702659e48155*I
mpc(real='-4.3084905090066053e+4824524', imag='8.6980843586535772e+4824524')
-4.30849050900660e4824524 + 8.69808435865358e4824524*I
sage: for z in [3,33,333,3333,33333]:
    mpmath.erf(-1-i*z); erf(-1.-n(z)*i)
....:     
mpc(real='330.81538696857206', imag='-443.38888183939281')
330.815386968572 - 443.388881839393*I
mpc(real='-2.0957487368415288e+468', imag='5.5629367605580166e+470')
-2.09574873684153e468 + 5.56293676055802e470*I
mpc(real='3.8930178706420656e+48153', imag='-1.8853741770265906e+48155')
3.89301787064207e48153 - 1.88537417702659e48155*I
mpc(real='4.3084905090066053e+4824524', imag='-8.6980843586535772e+4824524')
4.30849050900661e4824524 - 8.69808435865358e4824524*I

[kcrisman]

comment:15

More precisely:

sage: pari(3*i).erfc()
-1.76710569338983 E-16 - 1629.86732385786*I
sage: mpmath.erfc(3*i)
mpc(real='1.0', imag='-1629.9946226015657')
sage: 1-pari(3*i).erfc()
1.00000000000000 + 1629.86732385786*I
sage: mpmath.erf(3*i)
mpc(real='0.0', imag='1629.9946226015657')

[kcrisman]

comment:16

Now, in 5.1.beta6's gp I get

$ Downloads/sage-5.1.beta6/sage -gp
           GP/PARI CALCULATOR Version 2.5.1 (development git-5c5e253)
          i386 running darwin (x86-64/GMP-5.0.2 kernel) 64-bit version
                    compiled: Jun 25 2012, gcc-4.6.3 (GCC) 
                 (readline v6.2 enabled, extended help enabled)

                     Copyright (C) 2000-2011 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes 
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500509
? erfc(3*I)
%1 = 0.E-35 - 1629.5516567497094550267455532288372861*I

so I see that this is still in "our" version of Pari. Is this fixed upstream? Is it possible that this is not actually a bug in Pari, but in our use of it? Wolfram Alpha gives the mpmath answer, for what it's worth.

[kcrisman]

comment:17

I should also point out that we even doctest that it is wrong - due to my reviewer patch on this ticket :( so this is certainly a mea culpa among others.

[kcrisman]

comment:18

And apparently even weirder, courtesy of Jeff Denny.

sage: erf(i*1.42)

1.00000000000000 + 4.03986343036907*I

sage: import mpmath
sage: mpmath.erf(i*1.42)

mpc(real='0.0', imag='3.8217653554366318')

I'm opening a ticket for this, even if it gets fixed by one of the logical others (#1173 or #13050) because it is just mathematically wrong. This is #13193; continue discussion there.