Computing $\mathrm{erfi}(\theta)\exp(-\theta^2)$:

Question

I'm looking to compute $f(\theta):= \mathrm{erfi}(\theta)\exp(-\theta^2)$ as efficiently as possible, to double precision, with a fairly wide radius of converge.

Computing $\mathrm{erfi}(\theta)$ and then $\exp(-\theta^2)$ is disastrous, as one is in general huge and the other tiny, leading to overflow, then underflow, and a nan result.

The Taylor series of $f$, although having an infinite radius of convergence, has low order terms that easily overflow 64 bit doubles for large $\theta$, making it a pain to use.

What are the alternatives? C, C++, Python, and pure mathematics welcome.

(I have my own algorithm, using the Taylor series, that converges using long double intermediate values which converges for $|\theta| < 7$, but this isn't enough!)

Answer 1

For $s \ge 10$, the following approximation to $f(s)$ (from a Chebyshev series of $f(1/t)$) is accurate to within about $1.468 \times 10^{-17}$.

-0.125052986270239813e-16
+.564189583547791761764485499904/s
-0.1664326883508733997486686e-10/s^2
+.282094794823404373910675400690/s^3
-0.28755349061924926808699631e-6/s^4
+.423158100817100036834674291703/s^5
-0.554080675390895562280931052e-3/s^6
+1.07041448390054521448743009959/s^7
-.186553819701635257086895352275/s^8
+5.47622071471532066505331561467/s^9
-10.0148498076973648485542119208/s^10
+43.4633370516515753294154727476/s^11