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

SciPy Trac scipy-tickets@scipy....
Thu Nov 1 11:00:23 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):

> Example point:
>     z = -6 - 2.9j
>     Faddeeva_erf(z) == -1 + 0j; erf(z) == -1.00000000000006795+5
.53582896697010254e-14j;

This was just a simple bug in the snippet I posted, not any inherent
problem.  My underflow check should have been
{{
if (mRe_z2 < -750)
return (x >= 0 ? 1.0 : -1.0);
}}

(Note that this whole check is just an optimization, and is not required
for correctness.)  With this fix, my code gives:

-1.000000000000068-5.535828966970078e-14i

which is correct to nearly machine precision (according to WolframAlpha).

{{
{
double x = real(z), y = imag(z);
double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of
overflow
double mIm_z2 = -2*x*y; // Im(-z^2)
if (mRe_z2 < -750) // underflow
return (x >= 0 ? 1.0 : -1.0);

complex<double> mz2(mRe_z2, mIm_z2); // -z^2
if (x >= 0) {
if (x < 1e-3) {
if (fabs(y) < 1e-3)
goto taylor;
else if (fabs(mIm_z2 < 1e-3))
goto taylor_erfi;
}
return 1.0 - (exp(mz2) * Faddeeva_w(complex<double>(-y,x)));
}
else { // x < 0
if (x > -1e-3) { // duplicate from above to avoid fabs(x) call
if (fabs(y) < 1e-3)
goto taylor;
else if (fabs(mIm_z2 < 1e-3))
goto taylor_erfi;
}
return (exp(complex<double>(mRe_z2, mIm_z2))
}

// Use Taylor series for small |z|, to avoid cancellation inaccuracy
//     erf(z) 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - ...)
taylor:
return z * (1.1283791670955125739
+ mz2 * (0.37612638903183752464
+ mz2 * 0.11283791670955125739));

/* for small |x| and small |xy|,
use Taylor series to avoid cancellation inaccuracy:
erf(x+iy) = erf(iy)
+ 2*exp(y^2)/sqrt(pi) *
[ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
- i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
where:
erf(iy) = exp(y^2) * Im[w(y)]
*/
taylor_erfi:
if (x == 0) // important special case (also avoid NaN for large y)
return complex<double>(x, // preserve sign of 0
else {
double x2 = x*x, y2 = y*y;
double expy2 = exp(y2);
double sexpy2 = 1.1283791670955125739*expy2; // 2*exp(y^2)/sqrt(pi)
return complex<double>
(sexpy2 * x * (1 - x2 * (0.33333333333333333333
+ 0.66666666666666666667*y2)
+ x2*x2 * (0.1 +
y2*(0.4+0.13333333333333333333*y2))),