[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