inv_erf(x::real)::real = begin
if is_nan(x) || is_inf(x) || x < 0 || x > 1
reject("inv_erf: x must be finite and between 0 and 1; found x = ", x)
end
w = -log1m(square(x))
s = sqrt(w)
p = 0.0
if w < 5.0
w -= 2.5
p = 2.81022636e-08
p = fma(p, w, 3.43273939e-07)
p = fma(p, w, -3.5233877e-06)
p = fma(p, w, -4.39150654e-06)
p = fma(p, w, 0.00021858087)
p = fma(p, w, -0.00125372503)
p = fma(p, w, -0.00417768164)
p = fma(p, w, 0.246640727)
p = fma(p, w, 1.50140941)
else
w = s - 3.0
p = -0.000200214257
p = fma(p, w, 0.000100950558)
p = fma(p, w, 0.00134934322)
p = fma(p, w, -0.00367342844)
p = fma(p, w, 0.00573950773)
p = fma(p, w, -0.0076224613)
p = fma(p, w, 0.00943887047)
p = fma(p, w, 1.00167406)
p = fma(p, w, 2.83297682)
end
p * x
end