BigFloat erfcx (#77)

Fixes #75

Author: simonbyrne

src/erf.jl

@@ -281,3 +281,31 @@ function erfcinv(y::Float32)
 end
 
 erfcinv(x::Integer) = erfcinv(float(x))
+
+# MPFR has an open TODO item for this function
+# until then, we use [DLMF 7.12.1](https://dlmf.nist.gov/7.12.1) for the tail
+function erfcx(x::BigFloat)
+    if x <= (Clong == Int32 ? 0x1p15 : 0x1p30)
+        # any larger gives internal overflow
+        return exp(x^2)*erfc(x)
+    elseif !isfinite(x)
+        return 1/x
+    else
+        # asymptotic series
+        # starts to diverge at iteration i = 2^30 or 2^60
+        # final term will be < Γ(2*i+1)/(2^i * Γ(i+1)) / (2^(i+1))
+        # so good to (lgamma(2*i+1) - lgamma(i+1))/log(2) - 2*i - 1
+        #            ≈ 3.07e10 or 6.75e19 bits
+        # which is larger than the memory of the respective machines
+        ϵ = eps(BigFloat)/4
+        v = 1/(2*x*x)
+        k = 1
+        s = w = -k*v
+        while abs(w) > ϵ
+            k += 2
+            w *= -k*v
+            s += w
+        end
+        return (1+s)/(x*sqrt(oftype(x,pi)))
+    end
+end

test/runtests.jl

@@ -61,6 +61,10 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
 
     @test SF.erfinv(one(Int)) == SF.erfinv(1.0)
     @test SF.erfcinv(one(Int)) == SF.erfcinv(1.0)
+
+    @test SF.erfcx(1.8) ≈ SF.erfcx(big(1.8)) rtol=4*eps()
+    @test SF.erfcx(1.8e8) ≈ SF.erfcx(big(1.8e8)) rtol=4*eps()
+    @test SF.erfcx(1.8e88) ≈ SF.erfcx(big(1.8e88)) rtol=4*eps()
 end
 
 @testset "sine and cosine integrals" begin