RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc.

[Albert Chan]

It seems a general solver to get from erfc to its inverse is a bad idea.

erfc(-∞) = 2
erfc(-5) ≈ 2
erfc(0) = 1
erfc(+5) ≈ 0
erfc(+∞) = 0

The curve is very flat on both end, it is hard to get accurate inverse, without function's derivatives.

f = a - erfc(x) = (a-1) + erf(x)
f' = 2/sqrt(pi) * exp(-x*x)
f''/f' = -2x

Halley: x = x - f/(f' - (f''/f')*(f/2)) = x - f/(f' + x*f)

Code:

function ierfc_halley(a, verbal)
    local x = ierf_rough(1-a) -- guess
    repeat
        local f = a - erfc(x)
        local h = -f/(2/sqrt(pi)*exp(-x*x)+x*f)
        if verbal then print(x, h) end
        x = x+h
    until x+h*h == x
    return x
end

function ierf_rough(x)  -- max rel err = 0.002
    local b, neg = 4.33, x<0
    local y = log(1-x*x)/2
    y = sqrt(sqrt(b*b+y*((2-pi)*b+y)) - (b+y))
    y = finite(y) or 6
    return neg and -y or y
end

Again, ierfc is only rely on accurate erfc, which we have.

lua> x = 1e-9
lua> y = ierfc_halley(1-x)
lua> y
8.862269003885489e-10
lua> ierfc(1-x)
8.862269003885489e-10

Perfect accuracy. Now continue from previously, to get ierf(x)

lua> dydx = sqrt(pi)/2 * exp(y*y)
lua> h = x - (x-1+1)
lua> y + dydx * h
8.86226925452758e-10
lua> ierf(x)
8.862269254527581e-10