add inverse error functions erfinv and erfcinv

[stevengj]

This patch adds functions erfinv and erfcinv that compute the inverse of the error function (erf) and the complementary error function (erfc) for real arguments, respectively.

It uses the rational approximants described by Blair et al.. This is an independent, from-scratch implementation, so there should be no copyright concerns. (I use the tables of coefficients from Blair, but mathematical data of this sort is not generally subject to copyright.)

There are separate double- and single-precision versions, since the latter can be computed more cheaply using lower-degree approximants. Complex arguments are not supported (this would be a much harder problem and I'm not sure if it has been addressed in the literature).

[JeffBezanson]

Awesome! Major bonus points for including a reference.

[stevengj]

Thanks Jeff!

PS. In a quick benchmark (x = rand(10^7); @time y = erfinv(x);), it seems to be about 4× faster than Matlab's erfinv function (using a single thread), which in Matlab's case seems to use a compiled C or Fortran implementation. As usual, nice to be able to do this in pure Julia.

[JeffBezanson]

This is also a really good example of a non-trivial macro that is useful to and accessible to those more mathematically-inclined than programming-inclined.

[stevengj]

Yes, it was pleasant to be able to use a macro for this; might be nice to have an example of this sort in the manual to illustrate expression generation.

[StefanKarpinski]

This is one of my favorite pull requests ever.

[stevengj]

Hmm, also seems to be almost 3× faster than the erfinv in SciPy (which is implemented in C using the ndtri function from Cephes).

[JeffBezanson]

It may be because it loops over arrays of coefficients, instead of unrolling everything like your macro does.

[stevengj]

Here is an oddity: the single-precision version is slower than the double-precision version.

x = rand(10^7); xf = float32(x);
@time erfinv(x);
      elapsed time: 0.479699034 seconds
@time erfinv(xf); 
      elapsed time: 0.752342857 seconds

It looks like the culprit is that typeof(erfinv(xf)) gives Array{FloatingPoint,1} rather than Array{Float32,1}. Is this a bug in the @vectorize_1arg macro?

[stevengj]

@JeffBezanson, fixing the Float64 corrected the problem. The single-precision version is now consistently about 1.5× faster than the double-precision version.