add inverse error functions erfinv and erfcinv #2987

Merged
merged 1 commit into from May 1, 2013

Conversation

Projects
None yet
3 participants
Owner

stevengj commented May 1, 2013

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).

Owner

JeffBezanson commented May 1, 2013

Awesome! Major bonus points for including a reference.

JeffBezanson merged commit ccceddc into JuliaLang:master May 1, 2013

1 check was pending

default The Travis build is in progress
Details
Owner

stevengj commented May 1, 2013

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.

stevengj referenced this pull request May 1, 2013

Closed

0.2 release notes #2581

Owner

JeffBezanson commented May 1, 2013

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.

Owner

stevengj commented May 1, 2013

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.

Owner

StefanKarpinski commented May 1, 2013

This is one of my favorite pull requests ever.

Owner

stevengj commented May 1, 2013

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

Owner

JeffBezanson commented May 1, 2013

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

Owner

stevengj commented May 1, 2013

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?

Should be Float32.

Owner

stevengj commented May 1, 2013

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

stevengj deleted the stevengj:erfinv branch May 1, 2013

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment