Need help avoiding catastrophic cancellation or overflow

Specific Domains Numerics
12

The problem is that erfcx blows up exponentially for negative arguments, so if you have negative b and small c then it will quickly give Inf, and then you will get NaN from Inf - Inf.

However, there is a simple solution: for b < 0, do a change of variables x \to -x in your integral, which corresponds to flipping the signs of b, s, t and swapping the latter. That is:

function f(a,b,c,s,t)
    b < 0 && return f(a,-b,c,-t,-s)
    c == 0 && return exp(-a) * (exp(-b*s) - exp(-b*t)) / b
    return sqrt(pi/4c) * exp(-a) * (exp(-b*s-c*s^2) * erfcx((b+2c*s)/2√c) -
                                    exp(-b*t-c*t^2) * erfcx((b+2c*t)/2√c))
end

Then small c is not a problem regardless of the sign of b:

julia> f(1,2,1e-100,3,4)
0.00039423608073391826

julia> f(1,-2,1e-100,3,4)
474.1099996629409

(Note also that you only need to separate the c == 0 case, not some tricky-to-get-right precision-dependent threshold like c < 1e-16.)

4 Likes