[Scipy-tickets] [SciPy] #1757: new upstream of complex erf code (special.wofz) - performance, accuracy improvements
SciPy Trac
scipy-tickets@scipy....
Thu Nov 1 04:10:13 CDT 2012
#1757: new upstream of complex erf code (special.wofz) - performance, accuracy
improvements
---------------------------+------------------------------------------------
Reporter: stevenj | Owner: pv
Type: defect | Status: new
Priority: normal | Milestone: Unscheduled
Component: scipy.special | Version: devel
Keywords: |
---------------------------+------------------------------------------------
Comment(by pv):
@stevenj: first, a lot of thanks for your spending time on this. I'd be
really interested in seeing eps-relative accuracy versions of also erf and
erfc based on your code, as those can probably be implemented without too
much trouble.
The problems with the simplest approach of reusing results from Faddeeva W
directly are:
(i) Loss of relative precision occurs at |z| < eps where in this
implementation erf(z) == 0. This can be worked around relatively
painlessly by pushing the "1" in "1-erfc" inside the Faddeeva evaluation
routine in the Zaghloul algorithm (see the code linked).
(ii) The second loss of precision occurs in the imaginary part near real
axis in the expression (both in erf and erfc)
w = wofz(i z)
Im erf = exp(-x^2 + y^2) (-real(w) sin(2 x y) + imag(w) cos(2 x y))
Here, even if `w` is accurate down to machine epsilon, the subtraction
producing erf loses relative precision in the imaginary part completely.
Example point:
z = -6 - 2.9j
Faddeeva_erf(z) == -1 + 0j; erf(z) == -1.00000000000006795+5
.53582896697010254e-14j;
Using the mirror symmetry around imaginary axis indeed seems to make using
a separate expansion to avoid overflows unnecessary, I missed that, and it
seems your implementation indeed doesn't have problems with that.
--
Ticket URL: <http://projects.scipy.org/scipy/ticket/1757#comment:8>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
More information about the Scipy-tickets
mailing list