[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):

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

 Please try out:

 {{
 complex<double> Faddeeva_erf(complex<double> z, double relerr)
 {
   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))
             * Faddeeva_w(complex<double>(y,-x))) - 1.0;
   }

   // 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
                            exp(y*y) * ImFaddeeva_w(y));
   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))),
        expy2 * ImFaddeeva_w(y)
        - sexpy2*x2*y * (1 - x2*(0.5+0.33333333333333333333*y2)));
   }
 }
 }}

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


More information about the Scipy-tickets mailing list