Numerically stable derivative of logerf without catastrophic cancellation

Specific Domains Numerics
,
3

In particular, with this kind of thing, the “scaled” special functions are often your friends — here, \text{erfcx}(x) = e^{x^2} \text{erfc}(x), which gives various identities

\log \text{erf}(u,v) = \log \left[ e^{-u^2} \text{erfcx}(u) - e^{-v^2} \text{erfcx}(v) \right] \\ = -u^2 + \log \left[ \text{erfcx}(u) - e^{u^2-v^2} \text{erfcx}(v) \right] \\ = -u^2 + \log \left[ e^{u^2-v^2} \text{erfcx}(-v) -\text{erfcx}(-u) \right] \\ = -v^2 + \log \left[ e^{v^2-u^2} \text{erfcx}(u) - \text{erfcx}(v) \right] \\ = -v^2 + \log \left[ \text{erfcx}(-v) - e^{v^2-u^2} \text{erfcx}(-u) \right] \\

By selecting the appropriate version of these identities, you can cancel the -u^2 or -v^2 term analytically in your sum.

I whipped something together based on this, and it seems mostly right at first glance, but there may be some corner cases I haven’t handled correctly. Hopefully it will be a useful starting point for you, at least:

function dfun2(x,a,b)
    u, v  = a+x, b+x
    u²,v² = u^2, v^2
    c = 2 / sqrt(oftype(u+v, pi))
    if u² < v²
        if v < 0
            t = log1p(-exp(u²-v²)) - log(exp(u²-v²)*erfcx(-v) - erfcx(-u))
        else
            t = log1p(-exp(u²-v²)) - log(erfcx(u) - exp(u²-v²)*erfcx(v))
        end
        return -c * exp(t)
    else
        if u < 0
            t = log1p(-exp(v²-u²)) - log(erfcx(-v) - exp(v²-u²)*erfcx(-u))
        else
            t = log1p(-exp(v²-u²)) - log(exp(v²-u²)*erfcx(u) - erfcx(v))
        end
        return c * exp(t)
    end
end
5 Likes