integration - Gaussian integral variant $int_{-infty}^infty frac{e^{-x^2}}{1+a e^{-x}} dx$

I have been trying to compute this integral
$$\int_{-\infty}^{\infty} \frac{e^{-x^2}}{1+a e^{-x}} dx$$
quickly and to a high-degree of accuracy.

I have some partial results, for example for $n \in \mathbb N$
$$ \frac{1}{\sqrt\pi} \int_{-\infty}^{\infty} \frac{e^{-x^2}}{1+e^{-(x+n/2)}}dx = \frac{1}{2}\sum_{i=0}^{2n}(-1)^{i} e^{-i(2n-i)/4},$$
and for $a \in (0,1)$ we have
$$ \frac{1}{\sqrt\pi} \int_{0}^\infty \frac{e^{-x^2}}{1+ae^{-x}}dx = \frac{1}{2}\sum_{n=0}^\infty (-1)^n a^n e^{n^2/4} \mathrm{erfc}(n/2)$$
where $\mathrm{erfc}$ is the error function complement fuction
$$ \mathrm{erfc}(x) = 1-\mathrm{erf}(x),$$
which we obtain by using the substitution
$$ \frac{1}{1+ae^{-x}} = \sum_{n=0}^\infty (-1)^n a^n e^{-nx}.$$
This $\mathrm{erfc}$ series does not converge very quickly unless $a$ is very small, and wolframalpha can compute the integral very accurately between say $-100$ and $100$ for various values of $a$, so I must be missing a trick.

Edit
I found an infinite series for the whole integral, but it converges slowly
$$\frac{1}{\sqrt\pi}\int_{-\infty}^\infty \frac{e^{-x^2}}{1+e^{-(x+a)}}dx = \frac{e^{-a^2}}{2}\left[\sum_{n=0}^\infty (-1)^n \mathrm{erfcx}(-a+n/2) + \sum_{n=1}^\infty (-1)^{n+1} \mathrm{erfcx}(a+n/2)\right],$$
where $\mathrm{erfcx}$ is the scaled complement error function
$$ \mathrm{erfcx}(x) = e^{x^2} \mathrm{erfcx}(x) = e^{x^2}(1-\mathrm{erf}(x)).$$

Answer

Since your integrand is bell-shaped and there are no singularities for $a>0$, then the Euler-Maclaurin expansion shows that simple trapezoidal quadrature converges exponentially fast.

Using this code, I was able to evaluate your integral to float precision in 17 integrand evaluations, to double in 33 evaluations, long double in 65, and 50 binary digits in 260. Doubling of the number of correct digits by halving $h$ is observed.

Assuming that your integrand takes 100ns to evaluate, this implies that you can compute your integral to double precision in ~3.3 microseconds.

Here's the code I used

Real a = 2;
auto f = [&](Real t) { return exp(-t*t)/(1+a*exp(-t)); };
Real t_max = sqrt(-log(std::numeric_limits::epsilon()));
Real Q = boost::math::trapezoidal(f, -t_max, t_max, tol);
std::cout << "Q = " << Q << std::endl;

You presumably are using a different programming environment, but the trapezoidal rule is simple enough to code that this shouldn't be a problem.