#13193 closed defect (fixed)
erf evaluation is wrong along imaginary axis
Reported by: | kcrisman | Owned by: | jason, jkantor |
---|---|---|---|
Priority: | critical | Milestone: | sage-5.2 |
Component: | numerical | Keywords: | |
Cc: | Merged in: | sage-5.2.rc0 | |
Authors: | Benjamin Jones | Reviewers: | Karl-Dieter Crisman |
Report Upstream: | Reported upstream. No feedback yet. | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
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.
Attachments (1)
Change History (20)
comment:1 Changed 7 years ago by
comment:2 Changed 7 years ago by
This should be reported upstream to PARI: http://pari.math.u-bordeaux.fr/Bugs/Reporting.html
comment:3 Changed 7 years ago by
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 7 years ago by
- 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 7 years ago by
Update - this is now Pari bug 1335.
comment:6 Changed 7 years ago by
- 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 7 years ago by
Does this do
erf(3*I).n()
correctly now? If you are still passing doctests, then apparently not.
comment:8 Changed 7 years ago by
Yep,
sage: erf(3*I).n() 1629.99462260157*I
where is the doctest you are referring to?
comment:9 Changed 7 years ago by
In sage/symbolic/expression.pyx
I found one. I didn't run full tests (was waiting for patchbot). I'll update the patch.
comment:10 follow-up: ↓ 11 Changed 7 years ago by
@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 7 years ago by
- 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 theerf(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 7 years ago by
- 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 7 years ago by
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 7 years ago by
Thanks for the review.
Thanks for DOING this.
comment:15 Changed 7 years ago by
- Status changed from positive_review to needs_work
Please don't do this:
# 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 7 years ago by
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:18 Changed 7 years ago by
- Merged in set to sage-5.2.rc0
- Resolution set to fixed
- Status changed from positive_review to closed
comment:19 Changed 7 years ago by
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 that up until about
1.42*i
we are fine, first noted by Jeff Denny.