implement numerical evaluation of erf at complex arguments via mpmath algorithm

[williamstein]

When implemented, this would work:

sage: erf(2.0, algorithm='mpmath')
...

---

Apply patch trac_1173_complex_erf_v3.patch attached to this ticket to the Sage library.

CC: @benjaminfjones

Component: calculus

Keywords: sd35.5

Reviewer: Karl-Dieter Crisman

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

[williamstein]

comment:1

See also
#1174

[boothby]

comment:2

The following has a GPL'd implementation in c:

http://velveeta.che.wisc.edu/octave/lists/archive/octave-sources.1998/msg00032.html

[williamstein]

comment:3

This looks like the wya to go:

Josh Kantor:

scipy has an error function that takes complex arguments

sage:  import numpy, scipy
sage:  from scipy import special
sage:  j=numpy.complex(0,1)
sage: -j*float(sqrt(pi))*special.erf(2*j)/2
(16.45262776550727+0j)

Unfortunately numpy and sage's complex numbers are not compatible yet.

[kcrisman]

comment:4

Numpy and Sage numbers now seem to work well, at least to some extent:

sage: import numpy
sage: CC(numpy.complex(0,1))
1.00000000000000*I
sage: CC(numpy.complex(1,1))
1.00000000000000 + 1.00000000000000*I
sage: import scipy
sage: CC(scipy.special.erf(numpy.complex(1,1)))
1.31615128169795 + 0.190453469237835*I

Is it time to wrap this now, two years after opening?

[aghitza]

comment:5

This can also be done by mpmath in arbitrary precision, see

sage: import mpmath
sage: mpmath.erf?

[kcrisman]

comment:6

Replying to @aghitza:

This can also be done by mpmath in arbitrary precision, see

sage: import mpmath
sage: mpmath.erf?

And erf already has an _eval_f method, so maybe this could be changed to use mpmath? Even a loss in speed would be worth having the complex values. This could apply to error_fcn/erfc and others as well.

[kcrisman]

comment:7

Update - yes, we should definitely do this because of the complex values - just had a support request about not being able to do

sage: integrate(sin(x^2),(x,2,6))

and then use n() because of this (partly). See also #9044.

[sagetrac-dsm]

comment:8

Okay, here's v0.prealpha0.0_1a of a patch which uses mpmath to get complex and arbitrary-precision behaviour for erf. (If it works out we can do error_fcn the same way, as noted by kcrisman).

Unfortunately it does introduce a speed regression:

# 4.6.1, via pari
sage: timeit('float(erf(0))',number=1000) # must be float because 4.6.1 doesn't evaluate
1000 loops, best of 3: 109 µs per loop
sage: timeit('erf(0.1)',number=1000)
1000 loops, best of 3: 98.4 µs per loop
sage: timeit('erf(0.99)',number=1000)
1000 loops, best of 3: 98.5 µs per loop

# 4.6.2rc0 with patch
sage: timeit('erf(0)',number=1000)
1000 loops, best of 3: 68.3 µs per loop
sage: timeit('erf(0.1)',number=1000)
1000 loops, best of 3: 154 µs per loop
sage: timeit('erf(0.99)',number=1000)
1000 loops, best of 3: 165 µs per loop

I attempted to fix this by using the old approach for the default precision, but it usually wound up being more expensive to do the tests. Someone with better speed-fu is encouraged to look at it. I haven't finished a long testall yet, so there's probably something somewhere which breaks, but here's the first attempt.

[burcin]

Author: D. S. McNeil

[burcin]

comment:9

Attachment: trac_1173_add_complex_argument_erf.patch.gz

The patch looks really good. The only problem I spotted after a quick reading is that the code blocks in the documentation should be after ::. If the TESTS title is followed by text, you don't need the :: after that.

AFAIR, mpmath now supports extracting the precision from the input type in Sage. It is not necessary to do this in each calling function any more. I don't have an example handy though.

For examples of fast methods to check the type of the input you can look at sage/symbolic/pynac.pyx. If you replace the PY_TYPE_CHECK() calls with isinstance() you should get a reasonable speed from python.

[sagetrac-dsm]

comment:10

@burcin:

Okay, so just to be clear: EXAMPLES and TESTS don't need trailing colons, whether double or single (although double does seem pretty common; is it okay to use it?), but I should definitely use a double colon after test descriptions, e.g.

TESTS[::?]

Test that addition works::

    sage: 2+3
    5

Test that subtraction works::

    sage: 3-2
    1

instead of

TESTS::

Test that addition works:

    sage: 2+3
    5

Test that subtraction works:

    sage: 3-2
    1

That's easy enough to fix. The other two will take a bit more work, but I'll see what I can find. Examples to pattern after are very welcome. :-)

[kcrisman]

comment:11

Usually we use one colon after TEST if there is text after it. The double colon needs to be before any actual test block. Also, doing ./sage -docbuild reference html should give warning messages if it's wrong, if the file is in the reference manual (not all are).

[sagetrac-dsm]

comment:12

It seems like most of the speed decrease isn't due to anything I'm doing in _evalf_ but in _eval_, namely that I have one as opposed to the default BuiltinFunction.eval. When I switch back to the default one I recover the old speeds, except that I no longer know how to get erf(0)=0 without wrapping the whole thing, and it has strange ideas about which arguments shouldn't be evaluated (such as ComplexField). I also didn't manage to find whatever mpmath magic will allow me to avoid the if isinstance(parent, Parent) and hasattr(parent, 'prec') or try: parent.prec() approach.

I'm at a loss for ways to proceed that don't involve modifying function.pyx, cythonizing, or wrapping the entire function to patch the holes.

[burcin]

Reviewer: Burcin Erocal

[burcin]

comment:13

Replying to @sagetrac-dsm:

I'm at a loss for ways to proceed that don't involve modifying function.pyx, cythonizing, or wrapping the entire function to patch the holes.

The speed problems can be solved by replacing the zero test x == 0 with the example in #11143 comment:17. The patch attached to that ticket also contains an example call to mpmath without the prec clutter. AFAICT, the only remaining issue with this ticket is the docstring formatting.

I'd be happy to give positive review to a patch with these changes... :)

[kcrisman]

comment:14

I'd say that it should also add one very minor additional doctest, for erf(sqrt(2)). That would allow us to close #9044 without not having a doctest. I realize that symbolic input is checked, but it would be nice to have for completeness and addressing specific user requests :)

[kcrisman]

comment:15

Here's another thing that should be added, reported by one of the PREP workshop participants before this patch was in:

erf(4500).n()
1.0000000
erf(45000).n()
<boom>

So Pari was not handling big number in erf before. With this patch, though

sage: erf(4500).n()
1.00000000000000
sage: erf(45000).n()
1.00000000000000
sage: erf(4500000000).n()
1.00000000000000

Since none of the doctests in the current patch seem to test "big" real input to erf, we should add this too.

[kcrisman]

Work Issues: add erf(sqrt(2)) and erf(45000).n(), formatting, perhaps speed

[kcrisman]

comment:16

This should also solve #11626, I think?

[zimmermann6]

comment:17

Replying to @kcrisman:

This should also solve #11626, I think?

right, thus #11626 can be closed as soon as this ticket is closed
(or maybe right now?).

Douglas, can you implement the work issues in your patch?

Paul

[kcrisman]

Dependencies: #11513

[kcrisman]

comment:18

So, to review:

• Burcin and I want a small formatting change so the doc would build correctly.
• Burcin wants the speed to be fixed using the thing referenced at define symbolic functions for exponential integrals #11143, using the patch at add is_trivial_zero() method to symbolic expressions #11513.
• I want a doctest that checks big numbers will work (as in comment 15).
• Paul and I want a doctest for make the error function work in arbitrary precision #11626, to verify it is fixed.
• Burcin recommends calling mpmath without prec as at define symbolic functions for exponential integrals #11143.

This would make this depend on #11513. I'm still a little skeptical on this; will we really not get any wrong answers with that?

[zimmermann6]

comment:19

Replying to @kcrisman:

So, to review:

• I want a doctest that checks big numbers will work (as in comment 15).
• Paul and I want a doctest for make the error function work in arbitrary precision #11626, to verify it is fixed.

the current patch already contains examples with prec=100, both for real and complex numbers,
and thus is fine to me.

Paul

[zimmermann6]

comment:47

no I'm not on IRC. Ok, then I'll use 4.8.alpha6 instead.

Paul

[jdemeyer]

comment:48

Replying to @zimmermann6:

Hi Benjamin,

note that 4.8.rc1 is already out, should be newer than 4.8.alpha6.

It certainly isn't out yet (and maybe it will never even exist if rc0 solves all our problems). The easiest way to find out the latest development release is http://www.sagemath.org/download-latest.html, accessible as www.sagemath.org -> Download -> Development Release.
Alternatively, look at the sage-release announcements.

[benjaminfjones]

comment:49

OK, the feeling here at SD35.5 is to hijack this ticket (change the ticket description) to rebase the work done by DSM on top of Sage-4.8.alpha6 (rather than making a new ticket) so the work and documentation improvements here are preserved and compatible with #11948. Please pipe up if you're following this ticket and have an opinion.

[jdemeyer]

comment:50

Fine by me.

[kcrisman]

comment:52

See comments at #11948 for why we probably want to do this or #13050 as soon as possible.

[kcrisman]

comment:53

With respect to the move call patch, see also #13933.

[kcrisman]

comment:54

Apply trac_1173_complex_erf_v3.patch

[kcrisman]

comment:55

Pretty much everything here was done at #11948 or #13001. So let's close this (just a little I moved to #13050.

[kcrisman]

Changed dependencies from #11513 to none

[kcrisman]

Changed author from D. S. McNeil to none

[kcrisman]

Changed reviewer from Burcin Erocal, Benjamin Jones to Karl-Dieter Crisman