[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
improvements
---------------------------+------------------------------------------------
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
problem.
(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