title: Accurate computation of scaled complementary error function, erfcx()
body: <p>The (exponentially) scaled complementary error function, commonly designated by <code>erfcx</code>, is defined mathematically as erfcx(x) := e<sup>x<sup>2</sup></sup> erfc(x). It frequently occurs in diffusion problems in physics as well as chemistry. While some mathematical environments, such as <a href="https://www.mathworks.com/help/matlab/ref/erfcx.html" rel="nofollow noreferrer">MATLAB</a> and <a href="http://octave.sourceforge.net/octave/function/erfcx.html" rel="nofollow noreferrer">GNU Octave</a>, provide this function, it is absent from the C standard math library, which only provides <code>erf()</code> and <code>erfc()</code>.</p>
  
  <p>While it is possible to implement one's own <code>erfcx()</code> based directly on the mathematical definition, this only works over a limited input domain, because in the positive half-plane <code>erfc()</code> underflows for arguments of moderate magnitude, while <code>exp()</code> overflows, as noted in <a href="https://stackoverflow.com/questions/34723644/scaled-complementary-error-function-erfcxx-computation-avoiding-arithmetic-o">this question</a>, for example.</p>
  
  <p>For use with C, one could adapt some <code>erfcx()</code> open-source implementations such as the one in the <a href="http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package" rel="nofollow noreferrer">Faadeeva package</a>, as pointed to in responses to <a href="https://stackoverflow.com/questions/8962542/is-there-a-scaled-complementary-error-function-in-python-available">this question</a>. However, these implementations typically do not provide full accuracy for a given floating-point format. For example, tests with 2<sup>32</sup> test vectors show the maximum error of <code>erfcx()</code> as provided by the Faadeeva package to be 8.41 ulps in the positive half-plane and 511.68 ulps in the negative half-plane. </p>
  
  <p>A reasonable bound for an accurate implementation would be 4 ulps, corresponding to the accuracy bound of the math functions in the LA profile of <a href="https://software.intel.com/sites/products/documentation/doclib/mkl/vm/vmdata.htm" rel="nofollow noreferrer">Intel's Vector Math</a> library, which I have found to be a reasonable bound for non-trivial math function implementations that require both good accuracy and good performance.</p>
  
  <p>How could <code>erfcx()</code>, and the corresponding single-precision version, <code>erfcxf()</code>, be implemented accurately, while using only the C standard math library, and requiring no external libraries? We can assume that C's <code>float</code> nad <code>double</code> types are mapped to IEEE 754-2008 <code>binary32</code> and <code>binary64</code> floating-point types. Hardware support for the fused multiply-add operation (FMA) can be assumed, as this is supported by all major processor architectures at this time.</p>