Accurate computation of log error function difference For information-theoretic purposes I am interested in computing $f(a, b)=-\log_2\left(\text{erf}(b) - \text{erf}(a)\right)$ where $b > a$ and $\text{erf}$ is the error function.
Out-of-the box approximations (e.g scipy.special.erf) fail because $\text{erf}$ rapidly approaches $\pm 1$ so that catastrophic cancellation occurs inside the logarithm while $f$ itself is actually really well-behaved. Does there exist a reasonably accurate and efficient approximation of $f$ directly?
To make the question even more elementary, we can also approximate
$$g(a, b) = \ln\int_a^b e^{-t^2}\mathrm{d}t$$
and find $f = cg + d$ for some $c,d \in \mathbb{R}$.
 A: You can write $f(a,b)=-\log_2(\operatorname{erfc}(a)-\operatorname{erfc}(b))$.  That will eliminate the catastrophic cancellation.  If you want to use large arguments for $a,b$ you can use the expansion of $\operatorname{erfc}$ near $\infty$, which is
$$\operatorname{erfc}(x)\approx e^{-x^2}\pi^{-1/2}\left(\frac 1x-\frac 1{2x^3}+\frac 3{4x^5}+O(x^{-6})\right)$$
Let $c$ be the midpoint of your interval, so $c=\frac {a+b}2$ and pull out the factor $e^{-c^2}$.  You can take its logarithm analytically, getting $\log_2(e^{-c^2})=-\frac {c^2}{\log(2)}$ and none of the rest will overflow.  
Putting it all together, we have 
$$c=\frac {a+b}2\\
f(a,b)=-\log_2(\operatorname{erfc}(a)-\operatorname{erfc}(b))\\
\approx -\log_2\left(e^{-a^2}\pi^{-1/2}\left(\frac 1a-\frac 1{2a^3}+\frac 3{4a^5}\right)-e^{-b^2}\pi^{-1/2}\left(\frac 1b-\frac 1{2b^3}+\frac 3{4b^5}\right)\right)\\
=\frac {c^2}{\log(2)}-\log_2\left(e^{-a^2+c^2}\pi^{-1/2}\left(\frac 1a-\frac 1{2a^3}+\frac 3{4a^5}\right)-e^{-b^2+c^2}\pi^{-1/2}\left(\frac 1b-\frac 1{2b^3}+\frac 3{4b^5}\right)\right)$$
