[Scipy-tickets] [SciPy] #1757: new upstream of complex erf code (special.wofz) - performance, accuracy improvements

SciPy Trac scipy-tickets@scipy....
Thu Nov 1 10:39:57 CDT 2012

#1757: new upstream of complex erf code (special.wofz) - performance, accuracy
 Reporter:  stevenj        |       Owner:  pv         
     Type:  defect         |      Status:  new        
 Priority:  normal         |   Milestone:  Unscheduled
Component:  scipy.special  |     Version:  devel      
 Keywords:                 |  

Comment(by stevenj):

 Replying to [comment:8 pv]:
 > (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).

 A better way (faster, cleaner code) of fixing this problem is to simply
 switch to the Taylor expansion for small |z|.  (In fact, if z=x+iy, there
 is a Taylor expansion of erf that works well for small |x| and small |xy|
 even if |y| is not small.)

 I'll look into the second problem you mentioned.  Note that requiring low
 relative error in Re[erf] and Im[erf] taken by themselves is much more
 stringent than low relative error in erf(z) measured in the complex plane
 (by which criterion your second example is okay).  The former requirement
 may not always be possible, but it seems like it should be possible in the
 example you gave.

 For the next upstream version, my plan is to include helper routines to
 compute erf, erfc, etcetera.  I've already implemented the Taylor-
 expansion trick for erf, and hopefully there is a clean fix for the other

 (I also notice that in your version you added checks to the Algorithm 916
 sum to make sure the loop terminates in the case of NaN and underflow.  I
 am already aware of the problem if the user passes NaN and was planning to
 add a check for that.  Can you give any other inputs where the loops don't
 terminate?  In general, if you have bug fixes I would appreciate it if you
 would forward them to me so that they can be incorporated upstream.)

Ticket URL: <http://projects.scipy.org/scipy/ticket/1757#comment:9>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.

More information about the Scipy-tickets mailing list