Opened 8 years ago

Closed 8 years ago

# erf evaluation is wrong along imaginary axis

Reported by: Owned by: kcrisman jason, jkantor critical sage-5.2 numerical sage-5.2.rc0 Benjamin Jones Karl-Dieter Crisman Reported upstream. No feedback yet.

### Description

This was first noted at this ask.sagemath.org question. Here are some of the symptoms.

• Being off by a real part of 1 along the imaginary axis
```sage: for z in [3,33,333,3333]:
....:     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
```
• Starting to be wrong in the imaginary component as well as the number gets closer to 0
```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')
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')
```

Ideally this would be fixed by #1173 or #13050, but if there is something else we can do (and soon) that would be fine too.

### comment:1 Changed 8 years ago by kcrisman

Note that up until about `1.42*i` we are fine, first noted by Jeff Denny.

```sage: mpmath.erf(1.4*i)
mpc(real='0.0', imag='3.6569574831625342')
sage: erf(1.4*i)
3.65695748316253*I
```

### comment:2 Changed 8 years ago by jdemeyer

This should be reported upstream to PARI: http://pari.math.u-bordeaux.fr/Bugs/Reporting.html

### comment:3 Changed 8 years ago by kcrisman

This should be reported upstream to PARI: http://pari.math.u-bordeaux.fr/Bugs/Reporting.html

I'd been ... discouraged from this in the past. But I'll try. Part of the problem is that I don't have access to a more recent Pari, in case it's been fixed since our version.

### comment:4 Changed 8 years ago by kcrisman

• Report Upstream changed from N/A to Reported upstream. No feedback yet.

Okay, I tried this. I am not sure where to look for it yet, but when I find it online I'll put up a link.

### comment:5 Changed 8 years ago by kcrisman

Update - this is now Pari bug 1335.

### comment:6 Changed 8 years ago by benjaminfjones

• Authors set to Benjamin Jones
• Status changed from new to needs_review

The madness with `erf` has gone on long enough! Here is a patch to tide us over until the `algorithm=...` parameter to `_evalf_` can be implemented.

I only had to make one (tiny) change to one doctest after switching the numerical eval to mpmath. I also added a couple of doctests that reference this ticket.

### comment:7 Changed 8 years ago by kcrisman

Does this do

```erf(3*I).n()
```

correctly now? If you are still passing doctests, then apparently not.

### comment:8 Changed 8 years ago by benjaminfjones

Yep,

```sage: erf(3*I).n()
1629.99462260157*I
```

where is the doctest you are referring to?

### comment:9 Changed 8 years ago by benjaminfjones

In `sage/symbolic/expression.pyx` I found one. I didn't run full tests (was waiting for patchbot). I'll update the patch.

### Changed 8 years ago by benjaminfjones

update doctest in sage/symbolic/expression.pyx

### comment:10 follow-up: ↓ 11 Changed 8 years ago by benjaminfjones

@kcrisman, it looks like your reviewer patch on #11948 was merged in 4.8.alpha6, but then my patch at #13001 reorganized the doctests (moved a bunch from `__init__` to the class docstring) and the `erf(3*I).n()` doctest got removed in the process (oops!).

Should I add it back in do you think? I have `erf(3.0*I)` which calls `_evalf_` just the same.

### comment:11 in reply to: ↑ 10 Changed 8 years ago by kcrisman

• Reviewers set to Karl-Dieter Crisman

@kcrisman, it looks like your reviewer patch on #11948 was merged in 4.8.alpha6, but then my patch at #13001 reorganized the doctests (moved a bunch from `__init__` to the class docstring) and the `erf(3*I).n()` doctest got removed in the process (oops!).

That's hilarious! You must have done this at the Sage Days, and not mentioned it to me :) So you removed an erroneous doctest. Glad you nabbed the one in symbolic/expression.pyx, though.

Should I add it back in do you think? I have `erf(3.0*I)` which calls `_evalf_` just the same.

No, the one you changed is sufficient, since it does integer times I instead of real times I. Nice.

Should this be a blocker for Sage 5.2, since it is mathematically wrong? Anyway, on a first glance this seems like a good patch.

Before:

```sage: timeit('erf(320).n()')
625 loops, best of 3: 147 µs per loop
sage: timeit('erf(320*I).n()')
125 loops, best of 3: 5.49 ms per loop
sage: timeit('erf(-100).n()')
625 loops, best of 3: 160 µs per loop
sage: timeit('erf(100-100*i).n()')
125 loops, best of 3: 1.81 ms per loop
sage: timeit('erf(1+i).n()')
125 loops, best of 3: 1.65 ms per loop

```

After:

```sage: timeit('erf(320).n()')
625 loops, best of 3: 138 µs per loop
sage: timeit('erf(320*I).n()')
125 loops, best of 3: 2.15 ms per loop
sage: timeit('erf(-100).n()')
625 loops, best of 3: 137 µs per loop
sage: timeit('erf(100-100*i).n()')
125 loops, best of 3: 2.36 ms per loop
sage: timeit('erf(1+i).n()')
125 loops, best of 3: 1.73 ms per loop
```

Sometimes slower, sometimes faster, except on the one that was wrong before, where it's faster.

Looks pretty good. I'll see what the patchbot says.

### comment:12 Changed 8 years ago by kcrisman

• Status changed from needs_review to positive_review

I think this looks good. How on earth does patchbot test the entire library on Volker's computer in 540 seconds? My computer's a dinosaur in comparison - or, better, my computer's not a quick, lean, agile dinosaur in comparison.

### comment:13 Changed 8 years ago by benjaminfjones

I was wondering that myself :) Why bother running doctests locally when Volker's machine will power through them in under 10 minutes! Thanks for the review.

### comment:14 Changed 8 years ago by kcrisman

Thanks for the review.

Thanks for DOING this.

### comment:15 Changed 8 years ago by jdemeyer

• Status changed from positive_review to needs_work

```# avoid name conflicts with `parent` as a function parameter
from sage.structure.coerce import parent as s_parent
```

If you need to change something to avoid a conflict, change the function parameter instead.

### comment:16 Changed 8 years ago by kcrisman

I thought this hack was used in several other place, though, so that we're being consistent with other stuff? Also, in that case we would need some deprecation period, and we'd still need the hack to actually fix the problem.

### comment:17 Changed 8 years ago by jdemeyer

• Status changed from needs_work to positive_review

Agreed.

### comment:18 Changed 8 years ago by jdemeyer

• Merged in set to sage-5.2.rc0
• Resolution set to fixed
• Status changed from positive_review to closed

### comment:19 Changed 8 years ago by kcrisman

Incidentally, the Pari bug is now fixed, so if/when we ever implement various backends for special functions, we can use Pari for this one!

Note: See TracTickets for help on using tickets.