Which methods/algorithms for computation of the function $F$, where
$$F(a,b) = \int_a^b e^{-t^2}dt,\quad a\leq b,$$ are the best, i. e. fast and accurate?
I need to compute those integrals (multiplied by $1/\sqrt{2\pi}$) many times, therefore trapezoidal-like rules for numeric integration are rather slow.
The methods I am intrested in, do not have to be very accurate ($3$ or $4$ decimal places will do), so we can trade off some accuracy against speed.
Let us assume, you want to comput $$F(a,b) = \frac{2}{\sqrt \pi}\int_a^b e^{-t^2}dt,\quad a\leq b.$$ Then $$ F(0,b) = \operatorname{erf}(b),$$ where $\operatorname{erf}(b)$ is the error function.
Now simply split your domain into subintervalls and compute the value $F(0,x_i)$ for all values of your discretization, then $F(a,b) = F(0,x_i)-F(0,x_j)$ with $x_i$ close to $b$ and $x_j$ close to a.
So the question is how to quickly evalute the error function.
Wikipedia quotes numerical recipes for this. Quote:
Over the complete range of values, there is an approximation with a maximal error of $1.2\times 10^{-7}$, as follows: $$\operatorname{erf}(x)=\begin{cases} 1-\tau & \text{for }x\ge 0\\ \tau-1 & \text{for }x < 0 \end{cases}$$ with :$$\begin{align} \tau = {} & t\cdot\exp\left(-x^2-1.26551223+1.00002368 t+0.37409196 t^2+0.09678418 t^3\right.\\ & \left.{}-0.18628806 t^4+0.27886807 t^5-1.13520398 t^6+1.48851587\cdot t^7\right. \\ & \left.{}-0.82215223 t^8+0.17087277 t^9\right) \end{align}$$
Maybe you can use a lower order expansion if you don't need high accuracy. But precomputing the values seems for me the way to go.
Considering $$F(a,b) = \int_a^b e^{-t^2}dt=\frac{1}{2} \sqrt{\pi } \Big(\text{erf}(b)-\text{erf}(a)\Big)$$ There are many ways to approximate the error functions; I particularly enjoy $$\text{erf}(x)\approx\text{sgn}(x)\sqrt{1-e^{-\frac{x^2 \left(a x^2+\frac{4}{\pi }\right)}{a x^2+1}}}$$ Using $$a=\frac{8 (\pi -3)}{3 (4-\pi ) \pi }\approx 0.140012$$ leads to an error which is less than $0.00035$ for all $x$.
This can be improved using $a=0.1480092$ which leads to an error which is less than $0.00010$ for all $x$. I obtained this value minimizing $$\Phi(a)=\int_0^\infty \Big(\text{erf}(x)-f(x,a)\Big)^2\, dx$$
Edit
Inspired by the approximation I gave in my answer, I thought that we could improve it writing $$\text{erf}(x)\approx\text{sgn}(x)\sqrt{1-e^{\phi(x)}}$$ and build for $\phi(x)$ a Pade approximant at $x=0$. $$\phi(x)=\frac{12 \left(\pi ^2-10\right) x^2-60 (\pi -3) \pi }{\pi \left((120+\pi (7 \pi -60)) x^2+15 (\pi -3) \pi \right)}x^2$$ which is purely theoretical seems to be quite good.
A much better one (the literal values of the coefficients being too long to fit the page) could be $$\phi(x)=-\frac{0.0167527 x^4+0.160257 x^2+1.27324 }{0.0151778 x^4+0.155912 x^2+1.}x^2$$ It leads to a maximum error equal to $0.000012$ for all $x$.
