Does anyone know where I could find code for the "Inverse Error Function?" Freepascal/Delphi would be preferable but C/C++ would be fine too.
The TMath/DMath library did not have it :(
Here's an implementation of erfinv()
. Note that for it to work well, you also need a good implementation of erf()
.
function erfinv(const y: Double): Double;
//rational approx coefficients
const
a: array [0..3] of Double = ( 0.886226899, -1.645349621, 0.914624893, -0.140543331);
b: array [0..3] of Double = (-2.118377725, 1.442710462, -0.329097515, 0.012229801);
c: array [0..3] of Double = (-1.970840454, -1.624906493, 3.429567803, 1.641345311);
d: array [0..1] of Double = ( 3.543889200, 1.637067800);
const
y0 = 0.7;
var
x, z: Double;
begin
if not InRange(y, -1.0, 1.0) then begin
raise EInvalidArgument.Create('erfinv(y) argument out of range');
end;
if abs(y)=1.0 then begin
x := -y*Ln(0.0);
end else if y<-y0 then begin
z := sqrt(-Ln((1.0+y)/2.0));
x := -(((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0);
end else begin
if y<y0 then begin
z := y*y;
x := y*(((a[3]*z+a[2])*z+a[1])*z+a[0])/((((b[3]*z+b[3])*z+b[1])*z+b[0])*z+1.0);
end else begin
z := sqrt(-Ln((1.0-y)/2.0));
x := (((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0);
end;
//polish x to full accuracy
x := x - (erf(x) - y) / (2.0/sqrt(pi) * exp(-x*x));
x := x - (erf(x) - y) / (2.0/sqrt(pi) * exp(-x*x));
end;
Result := x;
end;
If you haven't got an implementation of erf()
then you can try this one converted to Pascal from Numerical Recipes. It's not accurate to double precision though.
function erfc(const x: Double): Double;
var
t,z,ans: Double;
begin
z := abs(x);
t := 1.0/(1.0+0.5*z);
ans := t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277)))))))));
if x>=0.0 then begin
Result := ans;
end else begin
Result := 2.0-ans;
end;
end;
function erf(const x: Double): Double;
begin
Result := 1.0-erfc(x);
end;
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With