I am trying to understand how I can numerically evaluate the imaginary error function with the incomplete gamma function for all $x \in \mathbb{R}$. I found here that the imaginary error function can be expressed as:
$$erfi(x) = \frac{\sqrt{-x^2}}{x}\big( Q( \frac{1}{2}, -x^2, \infty) -1 \big) $$
Where $Q$ is the regularized incomplete gamma function. For $x \in \mathbb{R}$ this equates to:
$$erfi(x) = i \big( Q( \frac{1}{2}, -x^2, \infty) -1 \big) $$
Now, I can evaluate the incomplete Gamma function, but using this method I won't get a real data type back. This confuses me, since for instance, another representation of the imaginary error function is:
$$erfi(x) = \frac{2}{\sqrt{\pi}} e^{x^{2}} D(x) $$
Where $D$ is the Dawson Function. This is entirely real valued! Also the scipy implementation gives me a real value too. [Yes I see the reference code, but it is pretty much impossible to follow. Besides, I want to understand how to do it myself.]
How can I implement the erfi function using the incomplete Gamma representation and still get a real value back?
There are two points to note:
The incomplete $\Gamma(1/2, -x^2)$ function is complex.
$\frac{\sqrt{-x^2}}{x}$ is $-i$ except for real $x>0$.
Here is an example for $x=-2$
$$\mathrm{erfi(-2)} = -18.56480241$$ $$\Gamma(1/2,-4) = 1.772453851-32.90525553i$$ $$Q(1/2,-4) = \Gamma(1/2,-4)/\sqrt{\pi}=1.000000000-18.56480241 i$$
$$\Big(Q(1/2,-4)-1\Big)\frac{\sqrt{-4}}{-2} = (-18.56480241 i) \times (-i) = -18.56480241$$