using ApproxFunFourier, ApproxFun, SpecialFunctions, SingularIntegralEquations, RiemannHilbert
Here we will consider the problem of computing the complementary error function $\mathrm{erfc}(z)$ for all $z \in \mathbb C$. Our approach is just a reinterpretation of $\newcommand{\E}{\mathrm{e}} \newcommand{\D}{\mathrm{d}} \newcommand{\I}{\mathrm{i}}$
A. C. Weideman, Computation of the Complex Error Function, SIAM Journal on Numerical Analysis, Vol. 31, pp. 1497-1518 (1994)
The complementary error function is defined by \begin{align*} \mathrm{erfc}(z) = \frac{2}{\sqrt \pi} \int_z^\infty \E^{-s^2} \D s. \end{align*}
The large $z$ asymptotics for this can be determined using integration by parts. After this, we define \begin{align*} \Phi(z) = \begin{cases} - \E^{z^2} \mathrm{erfc}(z) & \mathrm{Re}\, z > 0, \\ \E^{z^2} \left( 2 - \mathrm{erfc}(z) \right) & \mathrm{Re}\, z < 0. \end{cases} \end{align*}
And then it follows that $\Phi$ solves the following Riemann-Hilbert problem \begin{align*} \Phi_+(s) - \Phi_-(s) = 2 \E^{s^2}, \quad s \in \I \mathbb R, \quad \Phi(\infty) = 0. \end{align*} The orientation of $\I \mathbb R$ is that of increasing imaginary part.
The package ApproxFunFourier.jl includes functionality to perform an expansion of a function in the basis
$$ \begin{equation} f(x) \approx \sum_{j=-N}^N c_j \left( \frac{ L i + x}{ L i-x} \right)^j, \quad c_j = c_j(L). \end{equation} $$
And it makes sense to optimize over $L$ so that we minimize the number of terms required to reach machine accuracy.
f = x -> 2.0exp(-x^2)
n = 100_000
L = .1
for i = 1:50
global n,L
dom = PeriodicLine{false,Float64}(0.,L+.1)
F = Fun(f,Laurent(dom))
if length(F.coefficients) > n
break
end
n = length(F.coefficients)
L += .1
end
L
The package SingularIntegralEquations.jl implements the cauchy transform:
dom = PeriodicLine{false,Float64}(0.,L)
F = Fun(f,Laurent(dom))
my_erfc(z) = real(z) > 0 ? -exp(-z^2)*cauchy(F,-1im*z) : 2 - exp(-z^2)*cauchy(F,-1im*z)
x = -1.0; println(abs(my_erfc(x) - erfc(x))/abs(erfc(x)))
x = -100000000000.0+10000000im; println(abs(my_erfc(x) - erfc(x))/abs(erfc(x)))