From cd33177fec3335720748d520ba45d27553906a7d Mon Sep 17 00:00:00 2001 From: ziyixi Date: Wed, 13 Jun 2018 17:12:37 +0800 Subject: [PATCH] translate most functions in erf.jl in Julia's form --- src/erf.jl | 12 +- src/erf_libm.jl | 116 +++++ src/erf_openspecfun.jl | 1007 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 1129 insertions(+), 6 deletions(-) create mode 100644 src/erf_libm.jl create mode 100644 src/erf_openspecfun.jl diff --git a/src/erf.jl b/src/erf.jl index 99b707e7..04609d91 100644 --- a/src/erf.jl +++ b/src/erf.jl @@ -3,10 +3,12 @@ using Base.Math: @horner, libm using Base.MPFR: ROUNDING_MODE +include("erf_openspecfun.jl") +include("erf_libm.jl") + for f in (:erf, :erfc) @eval begin - ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x) - ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x) + ($f)(x::Float32) = Float32(($f)(Float64(x))) ($f)(x::Real) = ($f)(float(x)) ($f)(a::Float16) = Float16($f(Float32(a))) ($f)(a::Complex{Float16}) = Complex{Float16}($f(Complex{Float32}(a))) @@ -22,8 +24,7 @@ end for f in (:erf, :erfc, :erfcx, :erfi, :Dawson) fname = (f === :Dawson) ? :dawson : f @eval begin - ($fname)(z::Complex{Float64}) = Complex{Float64}(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) - ($fname)(z::Complex{Float32}) = Complex{Float32}(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex{Float64}(z), Float64(eps(Float32)))) + ($fname)(z::Complex{Float32}) = Complex{Float32}(($fname)(Float64(z))) ($fname)(z::Complex) = ($fname)(Complex{Float64}(z)) end end @@ -31,8 +32,7 @@ end for f in (:erfcx, :erfi, :Dawson) fname = (f === :Dawson) ? :dawson : f @eval begin - ($fname)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), x) - ($fname)(x::Float32) = Float32(ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), Float64(x))) + ($fname)(x::Float32) = Float32(($fname)(Float64(x))) ($fname)(x::Integer) = ($fname)(float(x)) end end diff --git a/src/erf_libm.jl b/src/erf_libm.jl new file mode 100644 index 00000000..42a96376 --- /dev/null +++ b/src/erf_libm.jl @@ -0,0 +1,116 @@ +""" + ErrFunApprox(x,iserfc,result) + +Compute the erf and erfc for x::Float64, translated from apple's libm. +""" +function ErrFunApprox(x::Float64, iserfc::Bool, result::Float64) + if isinf(x) + if x > 0 + return iserfc ? 0.0 : 1.0 + else + return iserfc ? 2.0 : -1.0 + end + end + if isnan(x) + return NaN + end + + # the largest representable double value + _HUGE = 6.71e+7 + + InvSqrtPI = 5.6418958354775628695e-1 + + a = [3.16112374387056560e+0, + 1.13864154151050156e+2, + 3.77485237685302021e+2, + 3.20937758913846947e+3, + 1.85777706184603153e-1] + + b = [2.36012909523441209e+1, + 2.44024637934444173e+2, + 1.28261652607737228e+3, + 2.84423683343917062e+3] + + ccc = [5.64188496988670089e-1, + 8.88314979438837594e+0, + 6.61191906371416295e+1, + 2.98635138197400131e+2, + 8.81952221241769090e+2, + 1.71204761263407058e+3, + 2.05107837782607147e+3, + 1.23033935479799725e+3, + 2.15311535474403846e-8] + + d = [1.57449261107098347e+1, + 1.17693950891312499e+2, + 5.37181101862009858e+2, + 1.62138957456669019e+3, + 3.29079923573345963e+3, + 4.36261909014324716e+3, + 3.43936767414372164e+3, + 1.23033935480374942e+3] + + pp = [3.05326634961232344e-1, + 3.60344899949804439e-1, + 1.25781726111229246e-1, + 1.60837851487422766e-2, + 6.58749161529837803e-4, + 1.63153871373020978e-2] + + qq = [2.56852019228982242e+0, + 1.87295284992346047e+0, + 5.27905102951428412e-1, + 6.05183413124413191e-2, + 2.33520497626869185e-3] + + y = abs(x) + + if y <= 0.46875e+0 + if y > 1.11e-16 + ysquared = y^2 + numerator=@horner(ysquared,0.,a[3],a[2],a[1],a[5]) + denominator=@horner(ysquared,0.,b[3],b[2],b[1],1.) + + result = y * (numerator + a[4]) / (denominator + b[4]); + else + result = y * a[4] / b[4] + end + if iserfc + result = 1.0 - result + end + return result + elseif y <= 4.0 + numerator = ccc[9] * y + denominator = y + for i in 1:7 + numerator = (numerator + ccc[i]) * y + denominator = (denominator + d[i]) * y + end + # refer to apple's libm, I don't exactly know the algorithm + result = (numerator + ccc[8]) / (denominator + d[8]) + ysquared = trunc(y * 16.0) / 16.0 + del = (y - ysquared) * (y + ysquared) + result = exp(-ysquared^2) * exp(-del) * result + else + if y >= _HUGE + result = InvSqrtPI / y + return result + end + ysquared = 1.0 / (y^2) + numerator = pp[6] * ysquared + denominator = ysquared + for i in 1:4 + numerator = (numerator + pp[i]) * ysquared + denominator = (denominator + qq[i]) * ysquared + end + result = ysquared * (numerator + pp[5]) / (denominator + qq[5]) + result = (InvSqrtPI - result) / y + ysquared = trunc(y * 16.0) / 16.0 + del = (y - ysquared) * (y + ysquared) + result = exp(-ysquared^2) * exp(-del) * result + end + return iserfc ? result : (0.5 - result) + 0.5 +end + +erf(x::Float64) = x >= 0 ? ErrFunApprox(x, false, 1.0) : -ErrFunApprox(x, false, 1.0) +erfc(x::Float64) = x >= 0 ? ErrFunApprox(x, true, 0.0) : 2.0 - ErrFunApprox(x, true, 0.0) diff --git a/src/erf_openspecfun.jl b/src/erf_openspecfun.jl new file mode 100644 index 00000000..fbaf872e --- /dev/null +++ b/src/erf_openspecfun.jl @@ -0,0 +1,1007 @@ +""" + erfcx(z) + +Compute the scaled complementary error function erfcx(z::Complex{Float64}). + +... +# Algorithm +- erfcx(z)=e^(z^2)*erfc(z)=w(im*z). +... +""" +function erfcx(z::Complex{Float64}) + return w(complex(-imag(z), real(z))) +end + + +""" + erf(x) + +Compute the error function erf(x::Float64). + +... +# Algorithm +- When |x| is relatively big, + - erf(x)=1-e^(-z^2)*w(im*z) if x>0. + - erf(x)=e^(-z^2)*w(im*z)-1 if x<0. +- When |x| is relatively small, use Taylor series for small |x|, to avoid cancellation inaccuracy. + - erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...) +... +""" +function erf_(x::Float64) + if isnan(x) + return NaN + elseif x^2 > 750 && x < 0 + return -1.0 + elseif x < -8e-2 + return exp(-x^2) * erfcx(-x) - 1.0 + elseif x < 8e-2 + return x * @horner(x^2, + 0.0052239776254421878422, + 0.026866170645131251760, + 0.11283791670955125739, + 0.37612638903183752464, + 1.1283791670955125739) + elseif x^2 < 750 + return 1.0 - exp(-x^2) * erfcx(x) + else + return 1.0 + end +end + + +""" + erf(z) + +Compute the error function erf(z::Complex{Float64}). + +... +# Algorithm +- When |z| is relatively big, + - erf(z)=1-e^(-z^2)*w(im*z) if real(z)>0 + - erf(z)=e^(-z^2)*w(im*z)-1 if real(z)<0. +- Use Taylor series for small |z|, to avoid cancellation inaccuracy. + - erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...). +- for small |x| and small |xy|, use Taylor series to avoid cancellation inaccuracy: + - erf(x+iy) = erf(iy) + 2*exp(y^2)/sqrt(pi) * + [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... + -i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ] + + where: erf(iy) = exp(y^2) * imag[w(y)] +... +""" +function erf(z::Complex{Float64}) + x, y = real(z), imag(z) + # isnan + if isnan(x) + return y == 0 ? NaN + 0.0im : NaN + NaN * im + elseif isnan(y) + return x == 0 ? 0.0 + NaN * im : NaN + NaN * im + end + # isinf + if isinf(y) + return x == 0 ? (y > 0 ? 0.0 + Inf * im : 0.0 - Inf * im) : NaN + NaN * im + elseif isinf(x) + return x > 0 ? 1.0 + 0.0im : -1.0 + 0.0im + end + # iszero + if y == 0 + return complex(erf_(x), y) + elseif x == 0 + return complex(x, y * y > 720 ? (y > 0 ? Inf : -Inf) : exp(y * y) * w_im(y)) + end + # normal value + mRe_z2 = (y - x) * (x + y) + mIm_z2 = -2 * x * y + if mRe_z2 < -750 + return x >= 0 ? 1.0 : -1.0 + end + #define taylor and taylor_erfi + @inline function taylor(z::Complex{Float64}) + mz2 = complex(mRe_z2, mIm_z2) + return z * @horner(mz2, + 1.1283791670955125739, + 0.37612638903183752464, + 0.11283791670955125739, + 0.026866170645131251760, + 0.0052239776254421878422) + end + @inline function taylor_erfi(z::Complex{Float64}) + x2, y2 = x^2, y^2 + expy2 = exp(y2) + return complex(expy2 * x * (1.1283791670955125739 + - x2 * (0.37612638903183752464 + + 0.75225277806367504925 * y2) + + x2^2 * (0.11283791670955125739 + + y2 * (0.45135166683820502956 + + 0.15045055561273500986 * y2))), + expy2 * (w_im(y) + - x2 * y * (1.1283791670955125739 + - x2 * (0.56418958354775628695 + + 0.37612638903183752464 * y2)))) + end + # start to compare + if x >= 0 + if x < 8e-2 + if(abs(y) < 1e-2) + return taylor(z) + elseif abs(mIm_z2 < 5e-3 && x < 5e-3) + return taylor_erfi(z) + end + end + return isinf(mIm_z2) ? NaN + NaN * im : 1.0 - exp(mRe_z2) * complex(cos(mIm_z2), sin(mIm_z2)) * w(complex(-y, x)) + else + if x > -8e-2 + if abs(y) < 1e-2 + return taylor(z) + elseif abs(mIm_z2) < 5e-3 && x > -5e-3 + return taylor_erfi(z) + end + end + return isinf(mIm_z2) ? NaN + NaN * im : exp(mRe_z2) * (complex(cos(mIm_z2), sin(mIm_z2)) * w(complex(y, -x))) - 1.0 + end +end + + +""" + erfi(x) + +Compute the imaginary error function erfi(x::Float64). + +... +# Algorithm +- erfi(x) = -i*erf(im*x). for real x, erfi(x)=e^(x^2)*imag(w(x)). +... +""" +function erfi(x::Float64) + return x * x > 720 ? (x > 0 ? Inf : -Inf) : exp(x * x) * w_im(x) +end + + +""" + erfi(z) + +Compute the imaginary error function erfi(z::Complex{Float64}). + +... +# Algorithm +- erfi(z) = -im*erf(im*z). +... +""" +function erfi(z::Complex{Float64}) + e = erf(complex(-imag(z), real(z))) + return complex(imag(e), -real(e)) +end + + +""" + erfc(x) + +Compute the complementary error function erfc(x::Float64). + +... +# Algorithm +- erfc(x)=e^(-x^2)*w(im*x) if x>=0 +- erfc(x)=2-e^(-x^2)*w(-im*x) if x<0 +... +""" +function erfc_(x::Float64) + if x^2 > 750 + return x >= 0 ? 0.0 : 2.0 + end + return x >= 0 ? exp(-x * x) * erfcx(x) : 2.0 - exp(-x * x) * erfcx(-x) +end + +""" + erfc(z) + +Compute the complementary error function erfc(x::Complex{Float64}). + +... +# Algorithm +- erfc(z)=e^(-z^2)*w(im*z) if real(z)>=0 +- erfc(z)=2-e^(-z^2)*w(-im*z) if real(z)<0 +... +""" +function erfc(z::Complex{Float64}) + x, y = real(z), imag(z) + + if x == 0.0 + return complex(1.0, y * y > 720 ? (y > 0 ? -Inf : Inf) : -exp(y * y) * w_im(y)) + end + if y == 0.0 + if x^2 > 750 + return complex(x >= 0 ? 0.0 : 2.0, -y) + end + return complex(x >= 0 ? exp(-x * x) * erfcx(x) : 2. - exp(-x * x) * erfcx(-x), -y) + end + + mRe_z2 = (y - x) * (x + y) + mIm_z2 = -2 * x * y + if mRe_z2 < -750 + return x >= 0 ? 0.0 : 2.0 + end + if x >= 0 + return exp(complex(mRe_z2, mIm_z2)) * w(complex(-y, x)) + else + return 2.0 - exp(complex(mRe_z2, mIm_z2)) * w(complex(y, -x)) + end +end + + +""" + dawson(x) + +Compute the Dawson function dawson(x::Float64) + +... +# Algorithm +- for real x, dawson(x)=sqrt(pi)/2*imag(w(x)) +... +""" +function dawson(x::Float64) + return sqrt(pi) / 2 * w_im(x) +end + + +""" + dawson(z) + +Compute the Dawson function dawson(z::Complex{Float64}) + +... +# Algorithm +- When |z| is relatively big, dawson(z)=sqrt(pi)/2*e^(-z^2)*erfi(z), which is equal to: + - dawson(z)=im*sqrt(pi)/2*(e^(-z^2)-w(z)) + - dawson(z)=im*sqrt(pi)/2*(w(-z)-e^(-z^2)) +- Use Taylor series for small |z|, to avoid cancellation inaccuracy + - dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ... +- For small |y| and small |xy|, use Taylor series to avoid cancellation inaccuracy: + - dawson(x + iy)= D + y^2 (D + x - 2Dx^2)+ y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3) + +im*y [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3) + +y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ... + + where D = dawson(x) + +- However, for large |x|, 2Dx -> 1 which gives cancellation problems in + this series (many of the leading terms cancel). So, for large |x|, + we need to substitute a continued-fraction expansion for D. + - dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...)))))) + + The 6 terms shown here seems to be the minimum needed to be + accurate as soon as the simpler Taylor expansion above starts + breaking down. Using this 6-term expansion, factoring out the + denominator, and simplifying with Maple, we obtain: + - real(dawson(x + iy)) * (-15 + 90x^2 - 60x^4 + 8x^6) / x + = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4 + - imag(dawson(x + iy)) * (-15 + 90x^2 - 60x^4 + 8x^6) / y + = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4 + +- Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction + expansion for the real part, and a 2-term expansion for the imaginary + part. (This avoids overflow problems for huge |x|.) This yields: + - real(dawson(x + iy)) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x) + - imag(dawson(x + iy)) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1) +... +""" +function dawson(z::Complex{Float64}) + const spi2 = 0.8862269254527580136490837416705725913990 + x, y = real(z), imag(z) + # isnan + if isnan(x) + return y == 0 ? NaN + 0.0im : NaN + NaN * im + elseif isnan(y) + return x == 0 ? 0.0 + NaN * im : NaN + NaN * im + end + # isinf + if isinf(x) && isinf(y) + return NaN + NaN * im + elseif isinf(x) + return 0.0 + 0.0im + elseif isinf(y) + return x == 0 ? 0.0 + y * im : NaN + y * im + end + # iszero + if y == 0 + return complex(spi2 * w_im(x), -y) + end + if x == 0 + y2 = y^2 + if y2 < 2.5e-5 + return complex(x, y * @horner(y2,1.,2 ./ 3.,4 ./ 15.)) + end + return complex(x, spi2 * (y >= 0 ? (exp(y2) - erfcx(y)) : (erfcx(-y) - exp(y2)))) + end + + mRe_z2 = (y - x) * (x + y) + mIm_z2 = -2 * x * y + mz2 = complex(mRe_z2, mIm_z2) + # define taylor and taylor_realaxis + @inline function taylor(z::Complex{Float64}) + return z * @horner(mz2,1.,2 ./ 3.,4 ./ 15.) + end + @inline function taylor_realaxis(z::Complex{Float64}) + x2 = x^2 + y2 = y^2 + if x2 > 1600 + if x2 > 25e14 + xy2 = (x * y)^2 + return complex((0.5 + y2 * (0.5 + 0.25 * y2 - 1 ./ 6. * xy2)) / x, + y * (-1 + y2 * (-2 ./ 3 .+ 2 ./ 15. * xy2 - 4 ./ 15. * y2)) / (2 * x2 - 1)) + end + return (1. / (-15 + x2 * (90 + x2 * (-60 + 8 * x2)))) * + complex(x * (33 + x2 * (-28 + 4 * x2) + y2 * (18 - 4 * x2 + 4 * y2)), + y * (-15 + x2 * (24 - 4 * x2) + y2 * (4 * x2 - 10 - 4 * y2))) + else + D = spi2 * w_im(x) + return complex(D + y2 * (D + x - 2 * D * x2) + y2 * y2 * (D * (0.5 - x2 * (2 - 2 ./ 3. * x2)) + x * (5 ./ 6 .- 1 ./ 3. * x2)), + y * (1 - 2 * D * x + y2 * 2 ./ 3. * (1 - x2 - D * x * (3 - 2 * x2)) + y2 * y2 * (4 ./ 15. - x2 * (0.6 - 2 ./ 15. * x2) - D * x * (1 - x2 * (4 ./ 3 .- 4 ./ 15. * x2))))) + end + end + # start to cpmpare + if y >= 0 + if y < 5e-3 + if abs(x) < 5e-3 + return taylor(z) + elseif abs(mIm_z2 < 5e-3) + return taylor_realaxis(z) + end + end + res = exp(mz2) - w(z) + return spi2 * complex(-imag(res), real(res)) + else + if y > -5e-3 + if abs(x) < 5e-3 + return taylor(z) + elseif abs(mIm_z2) < 5e-3 + return taylor_realaxis(z) + end + end + res = w(-z) - exp(mz2) + return spi2 * complex(-imag(res), real(res)) + end +end + + + +########################################################################################################### +###################### part to compute the Faddeeva function ############################################## +########################################################################################################### + +# precomputed table of expa2n2[n] = exp(-a2*n*n) +# for double-precision a2 = 0.26865... in function w(z::Complex{Float64}), below. +const expa2n2 = [ + 7.64405281671221563e-01, + 3.41424527166548425e-01, + 8.91072646929412548e-02, + 1.35887299055460086e-02, + 1.21085455253437481e-03, + 6.30452613933449404e-05, + 1.91805156577114683e-06, + 3.40969447714832381e-08, + 3.54175089099469393e-10, + 2.14965079583260682e-12, + 7.62368911833724354e-15, + 1.57982797110681093e-17, + 1.91294189103582677e-20, + 1.35344656764205340e-23, + 5.59535712428588720e-27, + 1.35164257972401769e-30, + 1.90784582843501167e-34, + 1.57351920291442930e-38, + 7.58312432328032845e-43, + 2.13536275438697082e-47, + 3.51352063787195769e-52, + 3.37800830266396920e-57, + 1.89769439468301000e-62, + 6.22929926072668851e-68, + 1.19481172006938722e-73, + 1.33908181133005953e-79, + 8.76924303483223939e-86, + 3.35555576166254986e-92, + 7.50264110688173024e-99, + 9.80192200745410268e-106, + 7.48265412822268959e-113, + 3.33770122566809425e-120, + 8.69934598159861140e-128, + 1.32486951484088852e-135, + 1.17898144201315253e-143, + 6.13039120236180012e-152, + 1.86258785950822098e-160, + 3.30668408201432783e-169, + 3.43017280887946235e-178, + 2.07915397775808219e-187, + 7.36384545323984966e-197, + 1.52394760394085741e-206, + 1.84281935046532100e-216, + 1.30209553802992923e-226, + 5.37588903521080531e-237, + 1.29689584599763145e-247, + 1.82813078022866562e-258, + 1.50576355348684241e-269, + 7.24692320799294194e-281, + 2.03797051314726829e-292, + 3.34880215927873807e-304, + 0.0 # underflow (also prevents reads past array end, below) +] + +""" + w(z) + +Compute the Faddeeva function, based on the algorithm in http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package . +""" +function w(z::Complex{Float64}) + @inline function sinc(x::Float64, sinx::Float64) + return abs(x) < 1e-4 ? 1 - (1 ./ 6.) * x * x : sinx / x + end + + @inline function sinh_taylor(x::Float64) + return x * (1 + (x * x) * (1 ./ 6. + + 1 ./ 120. * (x * x))) + end + + if real(z) == 0.0 + return complex(erfcx(imag(z)), real(z)) + elseif imag(z) == 0.0 + return complex(exp(-(real(z))^2), w_im(real(z))) + end + a = 0.518321480430085929872 + c = 0.329973702884629072537 + a2 = 0.268657157075235951582 + x = abs(real(z)) + y = imag(z) + ya = abs(y) + ret::Complex{Float64} = 0 + 0im + sum1, sum2, sum3, sum4, sum5 = 0.0, 0.0, 0.0, 0.0, 0.0 + + if ya > 7 || (x > 6 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28)) + ispi::Float64 = 0.56418958354775628694807945156 + xs = y < 0 ? -real(z) : real(z) + if x + ya > 4000 + if x + ya > 1e7 + if x > ya + yax = ya / xs + denom = ispi / (xs + yax * ya) + ret = complex(denom * yax, denom) + elseif isinf(ya) + return ((isnan(x) || y < 0) ? complex(NaN, NaN) : complex(0, 0)) + else + xya = xs / ya + denom = ispi / (xya * xs + ya) + ret = complex(denom, denom * xya) + end + else + dr = xs * xs - ya * ya - 0.5 + di = 2 * xs * ya + denom = ispi / (dr * dr + di * di) + ret = complex(denom * (xs * di - ya * dr), denom * (xs * dr + ya * di)) + end + else + c0, c1, c2, c3, c4 = 3.9, 11.398, 0.08254, 0.1421, 0.2023 + nu = floor(c0 + c1 / (c2 * x + c3 * ya + c4)) + wr, wi = xs, ya + nu = 0.5 * (nu - 1) + while nu > 0.4 + denom = nu / (wr * wr + wi * wi) + wr = xs - wr * denom + wi = ya + wi * denom + nu -= 0.5 + end + denom = ispi / (wr * wr + wi * wi) + ret = complex(denom * wi, denom * wr) + end + if y < 0 + return 2.0 * exp(complex((ya - xs) * (xs + ya), 2 * xs * y)) - ret + else + return ret + end + elseif x < 10 + prod2ax, prodm2ax, expx2 = 1.0, 1.0, 0.0 + + if isnan(y) + return complex(y, y) + end + + if x < 5e-4 + x2 = x^2 + expx2 = 1 - x2 * (1 - 0.5 * x2) + const ax2 = 1.036642960860171859744 * x + const exp2ax = 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667 * ax2)) + const expm2ax = 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667 * ax2)) + + n = 1 + while true + coef = expa2n2[n] * expx2 / (a2 * (n * n) + y * y) + prod2ax *= exp2ax + prodm2ax *= expm2ax + sum1 += coef + sum2 += coef * prodm2ax + sum3 += coef * prod2ax + sum5 += coef * (2 * a) * n * sinh_taylor((2 * a) * n * x) + if coef * prod2ax < eps(Float64) * sum3 + break + end + n += 1 + end + else + expx2 = exp(-x * x) + const exp2ax = exp((2 * a) * x) + const expm2ax = 1 / exp2ax + n = 1 + while true + const coef = expa2n2[n] * expx2 / (a2 * (n * n) + y * y) + prod2ax *= exp2ax + prodm2ax *= expm2ax + sum1 += coef + sum2 += coef * prodm2ax + sum4 += (coef * prodm2ax) * (a * n) + sum3 += coef * prod2ax + sum5 += (coef * prod2ax) * (a * n) + if (coef * prod2ax) * (a * n) < eps(Float64) * sum5 + break + end + n += 1 + end + end + const expx2erfcxy::Float64 = y > -6 ? expx2 * erfcx(y) : 2 * exp(y * y - x * x) + if y > 5 + sinxy::Float64 = sin(x * y) + ret = (expx2erfcxy - c * y * sum1) * cos(2 * x * y) + (c * x * expx2) * sinxy * sinc(x * y, sinxy) + else + xs = real(z) + const sinxy = sin(xs * y) + const sin2xy = sin(2 * xs * y) + const cos2xy = cos(2 * xs * y) + const coef1 = expx2erfcxy - c * y * sum1 + const coef2 = c * xs * expx2 + ret = complex(coef1 * cos2xy + coef2 * sinxy * sinc(xs * y, sinxy), + coef2 * sinc(2 * xs * y, sin2xy) - coef1 * sin2xy) + end + else + if isnan(x) + return complex(x, x) + end + if isnan(y) + return complex(y, y) + end + ret = exp(-x * x) + + n0 = floor(x / a + 0.5) + dx = a * n0 - x + sum3 = exp(-dx * dx) / (a2 * (n0 * n0) + y * y) + sum5 = a * n0 * sum3 + exp1 = exp(4 * a * dx) + exp1dn = 1 + dn = 1 + while n0 - dn > 0 + np = n0 + dn + nm = n0 - dn + tp = exp(-(a * dn + dx)^2) + tm = tp * (exp1dn *= exp1) + tp /= (a2 * (np * np) + y * y) + tm /= (a2 * (nm * nm) + y * y) + sum3 += tp + tm + sum5 += a * (np * tp + nm * tm) + if a * (np * tp + nm * tm) < eps(Float64) * sum5 + return ret + complex((0.5 * c) * y * (sum2 + sum3), + (0.5 * c) * copysign(sum5 - sum4, real(z))) + end + dn += 1 + end + while true + dn += 1 + np = n0 + dn + tp = exp(-(a * dn + dx)^2) / (a2 * (np * np) + y * y) + sum3 += tp + sum5 += a * np * tp + if a * np * tp < eps(Float64) * sum5 + return ret + complex((0.5 * c) * y * (sum2 + sum3), + (0.5 * c) * copysign(sum5 - sum4, real(z))) + end + end + end + return ret + complex((0.5 * c) * y * (sum2 + sum3), + (0.5 * c) * copysign(sum5 - sum4, real(z))) +end + +# erfcx(x) = exp(x^2) erfc(x) function, for real x, written by +# Steven G. Johnson, October 2012. + +# This function combines a few different ideas. + +# First, for x > 50, it uses a continued-fraction expansion (same as +# for the Faddeeva function, but with algebraic simplifications for z=i*x). + +# Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations, +# but with two twists: + +# a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation, +# inspired by a similar transformation in the octave-forge/specfun +# erfcx by Soren Hauberg, results in much faster Chebyshev convergence +# than other simple transformations I have examined. + +# b) Instead of using a single Chebyshev polynomial for the entire +# [0,1] y interval, we break the interval up into 100 equal +# subintervals, with a switch/lookup table, and use much lower +# degree Chebyshev polynomials in each subinterval. This greatly +# improves performance in my tests. + +# For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x), +# with the usual checks for overflow etcetera. + +# Performance-wise, it seems to be substantially faster than either +# the SLATEC DERFC function [or an erfcx function derived therefrom] +# or Cody's CALERF function (from netlib.org/specfun), while +# retaining near machine precision in accuracy. */ + +# Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x). + +# Uses a look-up table of 100 different Chebyshev polynomials +# for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated +# with the help of Maple and a little shell script. This allows +# the Chebyshev polynomials to be of significantly lower degree (about 1/4) +# compared to fitting the whole [0,1] interval with a single polynomial. +const caselist_erfcx_y100 = [0.70878032454106438663e-3 0.71234091047026302958e-3 0.35779077297597742384e-5 0.17403143962587937815e-7 0.81710660047307788845e-10 0.36885022360434957634e-12 0.15917038551111111111e-14 ; +0.21479143208285144230e-2 0.72686402367379996033e-3 0.36843175430938995552e-5 0.18071841272149201685e-7 0.85496449296040325555e-10 0.38852037518534291510e-12 0.16868473576888888889e-14 ; +0.36165255935630175090e-2 0.74182092323555510862e-3 0.37948319957528242260e-5 0.18771627021793087350e-7 0.89484715122415089123e-10 0.40935858517772440862e-12 0.17872061464888888889e-14 ; +0.51154983860031979264e-2 0.75722840734791660540e-3 0.39096425726735703941e-5 0.19504168704300468210e-7 0.93687503063178993915e-10 0.43143925959079664747e-12 0.18939926435555555556e-14 ; +0.66457513172673049824e-2 0.77310406054447454920e-3 0.40289510589399439385e-5 0.20271233238288381092e-7 0.98117631321709100264e-10 0.45484207406017752971e-12 0.20076352213333333333e-14 ; +0.82082389970241207883e-2 0.78946629611881710721e-3 0.41529701552622656574e-5 0.21074693344544655714e-7 0.10278874108587317989e-9 0.47965201390613339638e-12 0.21285907413333333333e-14 ; +0.98039537275352193165e-2 0.80633440108342840956e-3 0.42819241329736982942e-5 0.21916534346907168612e-7 0.10771535136565470914e-9 0.50595972623692822410e-12 0.22573462684444444444e-14 ; +0.11433927298290302370e-1 0.82372858383196561209e-3 0.44160495311765438816e-5 0.22798861426211986056e-7 0.11291291745879239736e-9 0.53386189365816880454e-12 0.23944209546666666667e-14 ; +0.13099232878814653979e-1 0.84167002467906968214e-3 0.45555958988457506002e-5 0.23723907357214175198e-7 0.11839789326602695603e-9 0.56346163067550237877e-12 0.25403679644444444444e-14 ; +0.14800987015587535621e-1 0.86018092946345943214e-3 0.47008265848816866105e-5 0.24694040760197315333e-7 0.12418779768752299093e-9 0.59486890370320261949e-12 0.26957764568888888889e-14 ; +0.16540351739394069380e-1 0.87928458641241463952e-3 0.48520195793001753903e-5 0.25711774900881709176e-7 0.13030128534230822419e-9 0.62820097586874779402e-12 0.28612737351111111111e-14 ; +0.18318536789842392647e-1 0.89900542647891721692e-3 0.50094684089553365810e-5 0.26779777074218070482e-7 0.13675822186304615566e-9 0.66358287745352705725e-12 0.30375273884444444444e-14 ; +0.20136801964214276775e-1 0.91936908737673676012e-3 0.51734830914104276820e-5 0.27900878609710432673e-7 0.14357976402809042257e-9 0.70114790311043728387e-12 0.32252476000000000000e-14 ; +0.21996459598282740954e-1 0.94040248155366777784e-3 0.53443911508041164739e-5 0.29078085538049374673e-7 0.15078844500329731137e-9 0.74103813647499204269e-12 0.34251892320000000000e-14 ; +0.23898877187226319502e-1 0.96213386835900177540e-3 0.55225386998049012752e-5 0.30314589961047687059e-7 0.15840826497296335264e-9 0.78340500472414454395e-12 0.36381553564444444445e-14 ; +0.25845480155298518485e-1 0.98459293067820123389e-3 0.57082915920051843672e-5 0.31613782169164830118e-7 0.16646478745529630813e-9 0.82840985928785407942e-12 0.38649975768888888890e-14 ; +0.27837754783474696598e-1 0.10078108563256892757e-2 0.59020366493792212221e-5 0.32979263553246520417e-7 0.17498524159268458073e-9 0.87622459124842525110e-12 0.41066206488888888890e-14 ; +0.29877251304899307550e-1 0.10318204245057349310e-2 0.61041829697162055093e-5 0.34414860359542720579e-7 0.18399863072934089607e-9 0.92703227366365046533e-12 0.43639844053333333334e-14 ; +0.31965587178596443475e-1 0.10566560976716574401e-2 0.63151633192414586770e-5 0.35924638339521924242e-7 0.19353584758781174038e-9 0.98102783859889264382e-12 0.46381060817777777779e-14 ; +0.34104450552588334840e-1 0.10823541191350532574e-2 0.65354356159553934436e-5 0.37512918348533521149e-7 0.20362979635817883229e-9 0.10384187833037282363e-11 0.49300625262222222221e-14 ; +0.36295603928292425716e-1 0.11089526167995268200e-2 0.67654845095518363577e-5 0.39184292949913591646e-7 0.21431552202133775150e-9 0.10994259106646731797e-11 0.52409949102222222221e-14 ; +0.38540888038840509795e-1 0.11364917134175420009e-2 0.70058230641246312003e-5 0.40943644083718586939e-7 0.22563034723692881631e-9 0.11642841011361992885e-11 0.55721092871111111110e-14 ; +0.40842225954785960651e-1 0.11650136437945673891e-2 0.72569945502343006619e-5 0.42796161861855042273e-7 0.23761401711005024162e-9 0.12332431172381557035e-11 0.59246802364444444445e-14 ; +0.43201627431540222422e-1 0.11945628793917272199e-2 0.75195743532849206263e-5 0.44747364553960993492e-7 0.25030885216472953674e-9 0.13065684400300476484e-11 0.63000532853333333334e-14 ; +0.45621193513810471438e-1 0.12251862608067529503e-2 0.77941720055551920319e-5 0.46803119830954460212e-7 0.26375990983978426273e-9 0.13845421370977119765e-11 0.66996477404444444445e-14 ; +0.48103121413299865517e-1 0.12569331386432195113e-2 0.80814333496367673980e-5 0.48969667335682018324e-7 0.27801515481905748484e-9 0.14674637611609884208e-11 0.71249589351111111110e-14 ; +0.50649709676983338501e-1 0.12898555233099055810e-2 0.83820428414568799654e-5 0.51253642652551838659e-7 0.29312563849675507232e-9 0.15556512782814827846e-11 0.75775607822222222221e-14 ; +0.53263363664388864181e-1 0.13240082443256975769e-2 0.86967260015007658418e-5 0.53662102750396795566e-7 0.30914568786634796807e-9 0.16494420240828493176e-11 0.80591079644444444445e-14 ; +0.55946601353500013794e-1 0.13594491197408190706e-2 0.90262520233016380987e-5 0.56202552975056695376e-7 0.32613310410503135996e-9 0.17491936862246367398e-11 0.85713381688888888890e-14 ; +0.58702059496154081813e-1 0.13962391363223647892e-2 0.93714365487312784270e-5 0.58882975670265286526e-7 0.34414937110591753387e-9 0.18552853109751857859e-11 0.91160736711111111110e-14 ; +0.61532500145144778048e-1 0.14344426411912015247e-2 0.97331446201016809696e-5 0.61711860507347175097e-7 0.36325987418295300221e-9 0.19681183310134518232e-11 0.96952238400000000000e-14 ; +0.64440817576653297993e-1 0.14741275456383131151e-2 0.10112293819576437838e-4 0.64698236605933246196e-7 0.38353412915303665586e-9 0.20881176114385120186e-11 0.10310784480000000000e-13 ; +0.67430045633130393282e-1 0.15153655418916540370e-2 0.10509857606888328667e-4 0.67851706529363332855e-7 0.40504602194811140006e-9 0.22157325110542534469e-11 0.10964842115555555556e-13 ; +0.70503365513338850709e-1 0.15582323336495709827e-2 0.10926868866865231089e-4 0.71182482239613507542e-7 0.42787405890153386710e-9 0.23514379522274416437e-11 0.11659571751111111111e-13 ; +0.73664114037944596353e-1 0.16028078812438820413e-2 0.11364423678778207991e-4 0.74701423097423182009e-7 0.45210162777476488324e-9 0.24957355004088569134e-11 0.12397238257777777778e-13 ; +0.76915792420819562379e-1 0.16491766623447889354e-2 0.11823685320041302169e-4 0.78420075993781544386e-7 0.47781726956916478925e-9 0.26491544403815724749e-11 0.13180196462222222222e-13 ; +0.80262075578094612819e-1 0.16974279491709504117e-2 0.12305888517309891674e-4 0.82350717698979042290e-7 0.50511496109857113929e-9 0.28122528497626897696e-11 0.14010889635555555556e-13 ; +0.83706822008980357446e-1 0.17476561032212656962e-2 0.12812343958540763368e-4 0.86506399515036435592e-7 0.53409440823869467453e-9 0.29856186620887555043e-11 0.14891851591111111111e-13 ; +0.87254084284461718231e-1 0.17999608886001962327e-2 0.13344443080089492218e-4 0.90900994316429008631e-7 0.56486134972616465316e-9 0.31698707080033956934e-11 0.15825697795555555556e-13 ; +0.90908120182172748487e-1 0.18544478050657699758e-2 0.13903663143426120077e-4 0.95549246062549906177e-7 0.59752787125242054315e-9 0.33656597366099099413e-11 0.16815130613333333333e-13 ; +0.94673404508075481121e-1 0.19112284419887303347e-2 0.14491572616545004930e-4 0.10046682186333613697e-6 0.63221272959791000515e-9 0.35736693975589130818e-11 0.17862931591111111111e-13 ; +0.98554641648004456555e-1 0.19704208544725622126e-2 0.15109836875625443935e-4 0.10567036667675984067e-6 0.66904168640019354565e-9 0.37946171850824333014e-11 0.18971959040000000000e-13 ; +0.10255677889470089531e0 0.20321499629472857418e-2 0.15760224242962179564e-4 0.11117756071353507391e-6 0.70814785110097658502e-9 0.40292553276632563925e-11 0.20145143075555555556e-13 ; +0.10668502059865093318e0 0.20965479776148731610e-2 0.16444612377624983565e-4 0.11700717962026152749e-6 0.74967203250938418991e-9 0.42783716186085922176e-11 0.21385479360000000000e-13 ; +0.11094484319386444474e0 0.21637548491908170841e-2 0.17164995035719657111e-4 0.12317915750735938089e-6 0.79376309831499633734e-9 0.45427901763106353914e-11 0.22696025653333333333e-13 ; +0.11534201115268804714e0 0.22339187474546420375e-2 0.17923489217504226813e-4 0.12971465288245997681e-6 0.84057834180389073587e-9 0.48233721206418027227e-11 0.24079890062222222222e-13 ; +0.11988259392684094740e0 0.23071965691918689601e-2 0.18722342718958935446e-4 0.13663611754337957520e-6 0.89028385488493287005e-9 0.51210161569225846701e-11 0.25540227111111111111e-13 ; +0.12457298393509812907e0 0.23837544771809575380e-2 0.19563942105711612475e-4 0.14396736847739470782e-6 0.94305490646459247016e-9 0.54366590583134218096e-11 0.27080225920000000000e-13 ; +0.12941991566142438816e0 0.24637684719508859484e-2 0.20450821127475879816e-4 0.15173366280523906622e-6 0.99907632506389027739e-9 0.57712760311351625221e-11 0.28703099555555555556e-13 ; +0.13443048593088696613e0 0.25474249981080823877e-2 0.21385669591362915223e-4 0.15996177579900443030e-6 0.10585428844575134013e-8 0.61258809536787882989e-11 0.30412080142222222222e-13 ; +0.13961217543434561353e0 0.26349215871051761416e-2 0.22371342712572567744e-4 0.16868008199296822247e-6 0.11216596910444996246e-8 0.65015264753090890662e-11 0.32210394506666666666e-13 ; +0.14497287157673800690e0 0.27264675383982439814e-2 0.23410870961050950197e-4 0.17791863939526376477e-6 0.11886425714330958106e-8 0.68993039665054288034e-11 0.34101266222222222221e-13 ; +0.15052089272774618151e0 0.28222846410136238008e-2 0.24507470422713397006e-4 0.18770927679626136909e-6 0.12597184587583370712e-8 0.73203433049229821618e-11 0.36087889048888888890e-13 ; +0.15626501395774612325e0 0.29226079376196624949e-2 0.25664553693768450545e-4 0.19808568415654461964e-6 0.13351257759815557897e-8 0.77658124891046760667e-11 0.38173420035555555555e-13 ; +0.16221449434620737567e0 0.30276865332726475672e-2 0.26885741326534564336e-4 0.20908350604346384143e-6 0.14151148144240728728e-8 0.82369170665974313027e-11 0.40360957457777777779e-13 ; +0.16837910595412130659e0 0.31377844510793082301e-2 0.28174873844911175026e-4 0.22074043807045782387e-6 0.14999481055996090039e-8 0.87348993661930809254e-11 0.42653528977777777779e-13 ; +0.17476916455659369953e0 0.32531815370903068316e-2 0.29536024347344364074e-4 0.23309632627767074202e-6 0.15899007843582444846e-8 0.92610375235427359475e-11 0.45054073102222222221e-13 ; +0.18139556223643701364e0 0.33741744168096996041e-2 0.30973511714709500836e-4 0.24619326937592290996e-6 0.16852609412267750744e-8 0.98166442942854895573e-11 0.47565418097777777779e-13 ; +0.18826980194443664549e0 0.35010775057740317997e-2 0.32491914440014267480e-4 0.26007572375886319028e-6 0.17863299617388376116e-8 0.10403065638343878679e-10 0.50190265831111111110e-13 ; +0.19540403413693967350e0 0.36342240767211326315e-2 0.34096085096200907289e-4 0.27479061117017637474e-6 0.18934228504790032826e-8 0.11021679075323598664e-10 0.52931171733333333334e-13 ; +0.20281109560651886959e0 0.37739673859323597060e-2 0.35791165457592409054e-4 0.29038742889416172404e-6 0.20068685374849001770e-8 0.11673891799578381999e-10 0.55790523093333333334e-13 ; +0.21050455062669334978e0 0.39206818613925652425e-2 0.37582602289680101704e-4 0.30691836231886877385e-6 0.21270101645763677824e-8 0.12361138551062899455e-10 0.58770520160000000000e-13 ; +0.21849873453703332479e0 0.40747643554689586041e-2 0.39476163820986711501e-4 0.32443839970139918836e-6 0.22542053491518680200e-8 0.13084879235290858490e-10 0.61873153262222222221e-13 ; +0.22680879990043229327e0 0.42366354648628516935e-2 0.41477956909656896779e-4 0.34300544894502810002e-6 0.23888264229264067658e-8 0.13846596292818514601e-10 0.65100183751111111110e-13 ; +0.23545076536988703937e0 0.44067409206365170888e-2 0.43594444916224700881e-4 0.36268045617760415178e-6 0.25312606430853202748e-8 0.14647791812837903061e-10 0.68453122631111111110e-13 ; +0.24444156740777432838e0 0.45855530511605787178e-2 0.45832466292683085475e-4 0.38352752590033030472e-6 0.26819103733055603460e-8 0.15489984390884756993e-10 0.71933206364444444445e-13 ; +0.25379911500634264643e0 0.47735723208650032167e-2 0.48199253896534185372e-4 0.40561404245564732314e-6 0.28411932320871165585e-8 0.16374705736458320149e-10 0.75541379822222222221e-13 ; +0.26354234756393613032e0 0.49713289477083781266e-2 0.50702455036930367504e-4 0.42901079254268185722e-6 0.30095422058900481753e-8 0.17303497025347342498e-10 0.79278273368888888890e-13 ; +0.27369129607732343398e0 0.51793846023052643767e-2 0.53350152258326602629e-4 0.45379208848865015485e-6 0.31874057245814381257e-8 0.18277905010245111046e-10 0.83144182364444444445e-13 ; +0.28426714781640316172e0 0.53983341916695141966e-2 0.56150884865255810638e-4 0.48003589196494734238e-6 0.33752476967570796349e-8 0.19299477888083469086e-10 0.87139049137777777779e-13 ; +0.29529231465348519920e0 0.56288077305420795663e-2 0.59113671189913307427e-4 0.50782393781744840482e-6 0.35735475025851713168e-8 0.20369760937017070382e-10 0.91262442613333333334e-13 ; +0.30679050522528838613e0 0.58714723032745403331e-2 0.62248031602197686791e-4 0.53724185766200945789e-6 0.37827999418960232678e-8 0.21490291930444538307e-10 0.95513539182222222221e-13 ; +0.31878680111173319425e0 0.61270341192339103514e-2 0.65564012259707640976e-4 0.56837930287837738996e-6 0.40035151353392378882e-8 0.22662596341239294792e-10 0.99891109760000000000e-13 ; +0.33130773722152622027e0 0.63962406646798080903e-2 0.69072209592942396666e-4 0.60133006661885941812e-6 0.42362183765883466691e-8 0.23888182347073698382e-10 0.10439349811555555556e-12 ; +0.34438138658041336523e0 0.66798829540414007258e-2 0.72783795518603561144e-4 0.63619220443228800680e-6 0.44814499336514453364e-8 0.25168535651285475274e-10 0.10901861383111111111e-12 ; +0.35803744972380175583e0 0.69787978834882685031e-2 0.76710543371454822497e-4 0.67306815308917386747e-6 0.47397647975845228205e-8 0.26505114141143050509e-10 0.11376390933333333333e-12 ; +0.37230734890119724188e0 0.72938706896461381003e-2 0.80864854542670714092e-4 0.71206484718062688779e-6 0.50117323769745883805e-8 0.27899342394100074165e-10 0.11862637614222222222e-12 ; +0.38722432730555448223e0 0.76260375162549802745e-2 0.85259785810004603848e-4 0.75329383305171327677e-6 0.52979361368388119355e-8 0.29352606054164086709e-10 0.12360253370666666667e-12 ; +0.40282355354616940667e0 0.79762880915029728079e-2 0.89909077342438246452e-4 0.79687137961956194579e-6 0.55989731807360403195e-8 0.30866246101464869050e-10 0.12868841946666666667e-12 ; +0.41914223158913787649e0 0.83456685186950463538e-2 0.94827181359250161335e-4 0.84291858561783141014e-6 0.59154537751083485684e-8 0.32441553034347469291e-10 0.13387957943111111111e-12 ; +0.43621971639463786896e0 0.87352841828289495773e-2 0.10002929142066799966e-3 0.89156148280219880024e-6 0.62480008150788597147e-8 0.34079760983458878910e-10 0.13917107176888888889e-12 ; +0.45409763548534330981e0 0.91463027755548240654e-2 0.10553137232446167258e-3 0.94293113464638623798e-6 0.65972492312219959885e-8 0.35782041795476563662e-10 0.14455745872000000000e-12 ; +0.47282001668512331468e0 0.95799574408860463394e-2 0.11135019058000067469e-3 0.99716373005509038080e-6 0.69638453369956970347e-8 0.37549499088161345850e-10 0.15003280712888888889e-12 ; +0.49243342227179841649e0 0.10037550043909497071e-1 0.11750334542845234952e-3 0.10544006716188967172e-5 0.73484461168242224872e-8 0.39383162326435752965e-10 0.15559069118222222222e-12 ; +0.51298708979209258326e0 0.10520454564612427224e-1 0.12400930037494996655e-3 0.11147886579371265246e-5 0.77517184550568711454e-8 0.41283980931872622611e-10 0.16122419680000000000e-12 ; +0.53453307979101369843e0 0.11030120618800726938e-1 0.13088741519572269581e-3 0.11784797595374515432e-5 0.81743383063044825400e-8 0.43252818449517081051e-10 0.16692592640000000000e-12 ; +0.55712643071169299478e0 0.11568077107929735233e-1 0.13815797838036651289e-3 0.12456314879260904558e-5 0.86169898078969313597e-8 0.45290446811539652525e-10 0.17268801084444444444e-12 ; +0.58082532122519320968e0 0.12135935999503877077e-1 0.14584223996665838559e-3 0.13164068573095710742e-5 0.90803643355106020163e-8 0.47397540713124619155e-10 0.17850211608888888889e-12 ; +0.60569124025293375554e0 0.12735396239525550361e-1 0.15396244472258863344e-3 0.13909744385382818253e-5 0.95651595032306228245e-8 0.49574672127669041550e-10 0.18435945564444444444e-12 ; +0.63178916494715716894e0 0.13368247798287030927e-1 0.16254186562762076141e-3 0.14695084048334056083e-5 0.10072078109604152350e-7 0.51822304995680707483e-10 0.19025081422222222222e-12 ; +0.65918774689725319200e0 0.14036375850601992063e-1 0.17160483760259706354e-3 0.15521885688723188371e-5 0.10601827031535280590e-7 0.54140790105837520499e-10 0.19616655146666666667e-12 ; +0.68795950683174433822e0 0.14741765091365869084e-1 0.18117679143520433835e-3 0.16392004108230585213e-5 0.11155116068018043001e-7 0.56530360194925690374e-10 0.20209663662222222222e-12 ; +0.71818103808729967036e0 0.15486504187117112279e-1 0.19128428784550923217e-3 0.17307350969359975848e-5 0.11732656736113607751e-7 0.58991125287563833603e-10 0.20803065333333333333e-12 ; +0.74993321911726254661e0 0.16272790364044783382e-1 0.20195505163377912645e-3 0.18269894883203346953e-5 0.12335161021630225535e-7 0.61523068312169087227e-10 0.21395783431111111111e-12 ; +0.78330143531283492729e0 0.17102934132652429240e-1 0.21321800585063327041e-3 0.19281661395543913713e-5 0.12963340087354341574e-7 0.64126040998066348872e-10 0.21986708942222222222e-12 ; +0.81837581041023811832e0 0.17979364149044223802e-1 0.22510330592753129006e-3 0.20344732868018175389e-5 0.13617902941839949718e-7 0.66799760083972474642e-10 0.22574701262222222222e-12 ; +0.85525144775685126237e0 0.18904632212547561026e-1 0.23764237370371255638e-3 0.21461248251306387979e-5 0.14299555071870523786e-7 0.69543803864694171934e-10 0.23158593688888888889e-12 ; +0.89402868170849933734e0 0.19881418399127202569e-1 0.25086793128395995798e-3 0.22633402747585233180e-5 0.15008997042116532283e-7 0.72357609075043941261e-10 0.23737194737777777778e-12 ; +0.93481333942870796363e0 0.20912536329780368893e-1 0.26481403465998477969e-3 0.23863447359754921676e-5 0.15746923065472184451e-7 0.75240468141720143653e-10 0.24309291271111111111e-12 ; +0.97771701335885035464e0 0.22000938572830479551e-1 0.27951610702682383001e-3 0.25153688325245314530e-5 0.16514019547822821453e-7 0.78191526829368231251e-10 0.24873652355555555556e-12 ; +] + + +""" + erfcx_y100(y100) + +Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x). + +Uses a look-up table of 100 different Chebyshev polynomials +for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated +with the help of Maple and a little shell script. This allows +the Chebyshev polynomials to be of significantly lower degree (about 1/4) +compared to fitting the whole [0,1] interval with a single polynomial. +""" +function erfcx_y100(y100::Float64) + if isnan(y100) + return NaN + end + theCase = Int64(floor(y100)) + 1 + if theCase <= 100 + t = 2 * y100 - (2 * theCase - 1) + return @horner(t, + caselist_erfcx_y100[theCase,1], + caselist_erfcx_y100[theCase,2], + caselist_erfcx_y100[theCase,3], + caselist_erfcx_y100[theCase,4], + caselist_erfcx_y100[theCase,5], + caselist_erfcx_y100[theCase,6], + caselist_erfcx_y100[theCase,7], + ) + else + return 1.0 + end +end + + +""" + erfcx(x) + +Compute the scaled complementary error function erfcx(x::Float64). +""" +function erfcx(x::Float64) + if x >= 0 + if x > 50 + ispi = 0.56418958354775628694807945156 + if x > 5e7 + return ispi / x + end + return ispi * ((x * x) * (x * x + 4.5) + 2) / (x * ((x * x) * (x * x + 5) + 3.75)) + end + return erfcx_y100(400 / (4 + x)) + else + return x < -26.7 ? Inf : (x < -6.1 ? 2 * exp(x * x) : 2 * exp(x * x) - erfcx_y100(400 / (4 - x))) + end +end + + + + +# Steven G. Johnson, October 2012. + +# Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x). + +# Uses a look-up table of 100 different Chebyshev polynomials +# for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated +# with the help of Maple and a little shell script. This allows +# the Chebyshev polynomials to be of significantly lower degree (about 1/30) +# compared to fitting the whole [0,1] interval with a single polynomial. +const caselist_w_im_y100 = [[ 0.28351593328822191546e-2 0.28494783221378400759e-2 0.14427470563276734183e-4 0.10939723080231588129e-6 0.92474307943275042045e-9 0.89128907666450075245e-11 0.92974121935111111110e-13 ], +[ 0.85927161243940350562e-2 0.29085312941641339862e-2 0.15106783707725582090e-4 0.11716709978531327367e-6 0.10197387816021040024e-8 0.10122678863073360769e-10 0.10917479678400000000e-12 ], +[ 0.14471159831187703054e-1 0.29703978970263836210e-2 0.15835096760173030976e-4 0.12574803383199211596e-6 0.11278672159518415848e-8 0.11547462300333495797e-10 0.12894535335111111111e-12 ], +[ 0.20476320420324610618e-1 0.30352843012898665856e-2 0.16617609387003727409e-4 0.13525429711163116103e-6 0.12515095552507169013e-8 0.13235687543603382345e-10 0.15326595042666666667e-12 ], +[ 0.26614461952489004566e-1 0.31034189276234947088e-2 0.17460268109986214274e-4 0.14582130824485709573e-6 0.13935959083809746345e-8 0.15249438072998932900e-10 0.18344741882133333333e-12 ], +[ 0.32892330248093586215e-1 0.31750557067975068584e-2 0.18369907582308672632e-4 0.15761063702089457882e-6 0.15577638230480894382e-8 0.17663868462699097951e-10 0.22126732680711111111e-12 0.30273474177737853668e-14 ], +[ 0.39317207681134336024e-1 0.32504779701937539333e-2 0.19354426046513400534e-4 0.17081646971321290539e-6 0.17485733959327106250e-8 0.20593687304921961410e-10 0.26917401949155555556e-12 0.38562123837725712270e-14 ], +[ 0.45896976511367738235e-1 0.33300031273110976165e-2 0.20423005398039037313e-4 0.18567412470376467303e-6 0.19718038363586588213e-8 0.24175006536781219807e-10 0.33059982791466666666e-12 0.49756574284439426165e-14 ], +[ 0.52640192524848962855e-1 0.34139883358846720806e-2 0.21586390240603337337e-4 0.20247136501568904646e-6 0.22348696948197102935e-8 0.28597516301950162548e-10 0.41045502119111111110e-12 0.65151614515238361946e-14 ], +[ 0.59556171228656770456e-1 0.35028374386648914444e-2 0.22857246150998562824e-4 0.22156372146525190679e-6 0.25474171590893813583e-8 0.34122390890697400584e-10 0.51593189879111111110e-12 0.86775076853908006938e-14 ], +[ 0.66655089485108212551e-1 0.35970095381271285568e-2 0.24250626164318672928e-4 0.24339561521785040536e-6 0.29221990406518411415e-8 0.41117013527967776467e-10 0.65786450716444444445e-12 0.11791885745450623331e-13 ], +[ 0.73948106345519174661e-1 0.36970297216569341748e-2 0.25784588137312868792e-4 0.26853012002366752770e-6 0.33763958861206729592e-8 0.50111549981376976397e-10 0.85313857496888888890e-12 0.16417079927706899860e-13 ], +[ 0.81447508065002963203e-1 0.38035026606492705117e-2 0.27481027572231851896e-4 0.29769200731832331364e-6 0.39336816287457655076e-8 0.61895471132038157624e-10 0.11292303213511111111e-11 0.23558532213703884304e-13 ], +[ 0.89166884027582716628e-1 0.39171301322438946014e-2 0.29366827260422311668e-4 0.33183204390350724895e-6 0.46276006281647330524e-8 0.77692631378169813324e-10 0.15335153258844444444e-11 0.35183103415916026911e-13 ], +[ 0.97121342888032322019e-1 0.40387340353207909514e-2 0.31475490395950776930e-4 0.37222714227125135042e-6 0.55074373178613809996e-8 0.99509175283990337944e-10 0.21552645758222222222e-11 0.55728651431872687605e-13 ], +[ 0.10532778218603311137e0 0.41692873614065380607e-2 0.33849549774889456984e-4 0.42064596193692630143e-6 0.66494579697622432987e-8 0.13094103581931802337e-9 0.31896187409777777778e-11 0.97271974184476560742e-13 ], +[ 0.11380523107427108222e0 0.43099572287871821013e-2 0.36544324341565929930e-4 0.47965044028581857764e-6 0.81819034238463698796e-8 0.17934133239549647357e-9 0.50956666166186293627e-11 0.18850487318190638010e-12 0.79697813173519853340e-14 ], +[ 0.12257529703447467345e0 0.44621675710026986366e-2 0.39634304721292440285e-4 0.55321553769873381819e-6 0.10343619428848520870e-7 0.26033830170470368088e-9 0.87743837749108025357e-11 0.34427092430230063401e-12 0.10205506615709843189e-13 ], +[ 0.13166276955656699478e0 0.46276970481783001803e-2 0.43225026380496399310e-4 0.64799164020016902656e-6 0.13580082794704641782e-7 0.39839800853954313927e-9 0.14431142411840000000e-10 0.42193457308830027541e-12 ], +[ 0.14109647869803356475e0 0.48088424418545347758e-2 0.47474504753352150205e-4 0.77509866468724360352e-6 0.18536851570794291724e-7 0.60146623257887570439e-9 0.18533978397305276318e-10 0.41033845938901048380e-13 -0.46160680279304825485e-13 ], +[ 0.15091057940548936603e0 0.50086864672004685703e-2 0.52622482832192230762e-4 0.95034664722040355212e-6 0.25614261331144718769e-7 0.80183196716888606252e-9 0.12282524750534352272e-10 -0.10531774117332273617e-11 -0.86157181395039646412e-13 ], +[ 0.16114648116017010770e0 0.52314661581655369795e-2 0.59005534545908331315e-4 0.11885518333915387760e-5 0.33975801443239949256e-7 0.82111547144080388610e-9 -0.12357674017312854138e-10 -0.24355112256914479176e-11 -0.75155506863572930844e-13 ], +[ 0.17185551279680451144e0 0.54829002967599420860e-2 0.67013226658738082118e-4 0.14897400671425088807e-5 0.40690283917126153701e-7 0.44060872913473778318e-9 -0.52641873433280000000e-10 -0.30940587864543343124e-11 ], +[ 0.18310194559815257381e0 0.57701559375966953174e-2 0.76948789401735193483e-4 0.18227569842290822512e-5 0.41092208344387212276e-7 -0.44009499965694442143e-9 -0.92195414685628803451e-10 -0.22657389705721753299e-11 0.10004784908106839254e-12 ], +[ 0.19496527191546630345e0 0.61010853144364724856e-2 0.88812881056342004864e-4 0.21180686746360261031e-5 0.30652145555130049203e-7 -0.16841328574105890409e-8 -0.11008129460612823934e-9 -0.12180794204544515779e-12 0.15703325634590334097e-12 ], +[ 0.20754006813966575720e0 0.64825787724922073908e-2 0.10209599627522311893e-3 0.22785233392557600468e-5 0.73495224449907568402e-8 -0.29442705974150112783e-8 -0.94082603434315016546e-10 0.23609990400179321267e-11 0.14141908654269023788e-12 ], +[ 0.22093185554845172146e0 0.69182878150187964499e-2 0.11568723331156335712e-3 0.22060577946323627739e-5 -0.26929730679360840096e-7 -0.38176506152362058013e-8 -0.47399503861054459243e-10 0.40953700187172127264e-11 0.69157730376118511127e-13 ], +[ 0.23524827304057813918e0 0.74063350762008734520e-2 0.12796333874615790348e-3 0.18327267316171054273e-5 -0.66742910737957100098e-7 -0.40204740975496797870e-8 0.14515984139495745330e-10 0.44921608954536047975e-11 -0.18583341338983776219e-13 ], +[ 0.25058626331812744775e0 0.79377285151602061328e-2 0.13704268650417478346e-3 0.11427511739544695861e-5 -0.10485442447768377485e-6 -0.34850364756499369763e-8 0.72656453829502179208e-10 0.36195460197779299406e-11 -0.84882136022200714710e-13 ], +[ 0.26701724900280689785e0 0.84959936119625864274e-2 0.14112359443938883232e-3 0.17800427288596909634e-6 -0.13443492107643109071e-6 -0.23512456315677680293e-8 0.11245846264695936769e-9 0.19850501334649565404e-11 -0.11284666134635050832e-12 ], +[ 0.28457293586253654144e0 0.90581563892650431899e-2 0.13880520331140646738e-3 -0.97262302362522896157e-6 -0.15077100040254187366e-6 -0.88574317464577116689e-9 0.12760311125637474581e-9 0.20155151018282695055e-12 -0.10514169375181734921e-12 ], +[ 0.30323425595617385705e0 0.95968346790597422934e-2 0.12931067776725883939e-3 -0.21938741702795543986e-5 -0.15202888584907373963e-6 0.61788350541116331411e-9 0.11957835742791248256e-9 -0.12598179834007710908e-11 -0.75151817129574614194e-13 ], +[ 0.32292521181517384379e0 0.10082957727001199408e-1 0.11257589426154962226e-3 -0.33670890319327881129e-5 -0.13910529040004008158e-6 0.19170714373047512945e-8 0.94840222377720494290e-10 -0.21650018351795353201e-11 -0.37875211678024922689e-13 ], +[ 0.34351233557911753862e0 0.10488575435572745309e-1 0.89209444197248726614e-4 -0.43893459576483345364e-5 -0.11488595830450424419e-6 0.28599494117122464806e-8 0.61537542799857777779e-10 -0.24935749227658002212e-11 ], +[ 0.36480946642143669093e0 0.10789304203431861366e-1 0.60357993745283076834e-4 -0.51855862174130669389e-5 -0.83291664087289801313e-7 0.33898011178582671546e-8 0.27082948188277716482e-10 -0.23603379397408694974e-11 0.19328087692252869842e-13 ], +[ 0.38658679935694939199e0 0.10966119158288804999e-1 0.27521612041849561426e-4 -0.57132774537670953638e-5 -0.48404772799207914899e-7 0.35268354132474570493e-8 -0.32383477652514618094e-11 -0.19334202915190442501e-11 0.32333189861286460270e-13 ], +[ 0.40858275583808707870e0 0.11006378016848466550e-1 -0.76396376685213286033e-5 -0.59609835484245791439e-5 -0.13834610033859313213e-7 0.33406952974861448790e-8 -0.26474915974296612559e-10 -0.13750229270354351983e-11 0.36169366979417390637e-13 ], +[ 0.43051714914006682977e0 0.10904106549500816155e-1 -0.43477527256787216909e-4 -0.59429739547798343948e-5 0.17639200194091885949e-7 0.29235991689639918688e-8 -0.41718791216277812879e-10 -0.81023337739508049606e-12 0.33618915934461994428e-13 ], +[ 0.45210428135559607406e0 0.10659670756384400554e-1 -0.78488639913256978087e-4 -0.56919860886214735936e-5 0.44181850467477733407e-7 0.23694306174312688151e-8 -0.49492621596685443247e-10 -0.31827275712126287222e-12 0.27494438742721623654e-13 ], +[ 0.47306491195005224077e0 0.10279006119745977570e-1 -0.11140268171830478306e-3 -0.52518035247451432069e-5 0.64846898158889479518e-7 0.17603624837787337662e-8 -0.51129481592926104316e-10 0.62674584974141049511e-13 0.20055478560829935356e-13 ], +[ 0.49313638965719857647e0 0.97725799114772017662e-2 -0.14122854267291533334e-3 -0.46707252568834951907e-5 0.79421347979319449524e-7 0.11603027184324708643e-8 -0.48269605844397175946e-10 0.32477251431748571219e-12 0.12831052634143527985e-13 ], +[ 0.51208057433416004042e0 0.91542422354009224951e-2 -0.16726530230228647275e-3 -0.39964621752527649409e-5 0.88232252903213171454e-7 0.61343113364949928501e-9 -0.42516755603130443051e-10 0.47910437172240209262e-12 0.66784341874437478953e-14 ], +[ 0.52968945458607484524e0 0.84400880445116786088e-2 -0.18908729783854258774e-3 -0.32725905467782951931e-5 0.91956190588652090659e-7 0.14593989152420122909e-9 -0.35239490687644444445e-10 0.54613829888448694898e-12 ], +[ 0.54578857454330070965e0 0.76474155195880295311e-2 -0.20651230590808213884e-3 -0.25364339140543131706e-5 0.91455367999510681979e-7 -0.23061359005297528898e-9 -0.27512928625244444444e-10 0.54895806008493285579e-12 ], +[ 0.56023851910298493910e0 0.67938321739997196804e-2 -0.21956066613331411760e-3 -0.18181127670443266395e-5 0.87650335075416845987e-7 -0.51548062050366615977e-9 -0.20068462174044444444e-10 0.50912654909758187264e-12 ], +[ 0.57293478057455721150e0 0.58965321010394044087e-2 -0.22841145229276575597e-3 -0.11404605562013443659e-5 0.81430290992322326296e-7 -0.71512447242755357629e-9 -0.13372664928000000000e-10 0.44461498336689298148e-12 ], +[ 0.58380635448407827360e0 0.49717469530842831182e-2 -0.23336001540009645365e-3 -0.51952064448608850822e-6 0.73596577815411080511e-7 -0.84020916763091566035e-9 -0.76700972702222222221e-11 0.36914462807972467044e-12 ], +[ 0.59281340237769489597e0 0.40343592069379730568e-2 -0.23477963738658326185e-3 0.34615944987790224234e-7 0.64832803248395814574e-7 -0.90329163587627007971e-9 -0.30421940400000000000e-11 0.29237386653743536669e-12 ], +[ 0.59994428743114271918e0 0.30976579788271744329e-2 -0.23308875765700082835e-3 0.51681681023846925160e-6 0.55694594264948268169e-7 -0.91719117313243464652e-9 0.53982743680000000000e-12 0.22050829296187771142e-12 ], +[ 0.60521224471819875444e0 0.21732138012345456060e-2 -0.22872428969625997456e-3 0.92588959922653404233e-6 0.46612665806531930684e-7 -0.89393722514414153351e-9 0.31718550353777777778e-11 0.15705458816080549117e-12 ], +[ 0.60865189969791123620e0 0.12708480848877451719e-2 -0.22212090111534847166e-3 0.12636236031532793467e-5 0.37904037100232937574e-7 -0.84417089968101223519e-9 0.49843180828444444445e-11 0.10355439441049048273e-12 ], +[ 0.61031580103499200191e0 0.39867436055861038223e-3 -0.21369573439579869291e-3 0.15339402129026183670e-5 0.29787479206646594442e-7 -0.77687792914228632974e-9 0.61192452741333333334e-11 0.60216691829459295780e-13 ], +[ 0.61027109047879835868e0 -0.43680904508059878254e-3 -0.20383783788303894442e-3 0.17421743090883439959e-5 0.22400425572175715576e-7 -0.69934719320045128997e-9 0.67152759655111111110e-11 0.26419960042578359995e-13 ], +[ 0.60859639489217430521e0 -0.12305921390962936873e-2 -0.19290150253894682629e-3 0.18944904654478310128e-5 0.15815530398618149110e-7 -0.61726850580964876070e-9 0.68987888999111111110e-11 ], +[ 0.60537899426486075181e0 -0.19790062241395705751e-2 -0.18120271393047062253e-3 0.19974264162313241405e-5 0.10055795094298172492e-7 -0.53491997919318263593e-9 0.67794550295111111110e-11 -0.17059208095741511603e-13 ], +[ 0.60071229457904110537e0 -0.26795676776166354354e-2 -0.16901799553627508781e-3 0.20575498324332621581e-5 0.51077165074461745053e-8 -0.45536079828057221858e-9 0.64488005516444444445e-11 -0.29311677573152766338e-13 ], +[ 0.59469361520112714738e0 -0.33308208190600993470e-2 -0.15658501295912405679e-3 0.20812116912895417272e-5 0.93227468760614182021e-9 -0.38066673740116080415e-9 0.59806790359111111110e-11 -0.36887077278950440597e-13 ], +[ 0.58742228631775388268e0 -0.39321858196059227251e-2 -0.14410441141450122535e-3 0.20743790018404020716e-5 -0.25261903811221913762e-8 -0.31212416519526924318e-9 0.54328422462222222221e-11 -0.40864152484979815972e-13 ], +[ 0.57899804200033018447e0 -0.44838157005618913447e-2 -0.13174245966501437965e-3 0.20425306888294362674e-5 -0.53330296023875447782e-8 -0.25041289435539821014e-9 0.48490437205333333334e-11 -0.42162206939169045177e-13 ], +[ 0.56951968796931245974e0 -0.49864649488074868952e-2 -0.11963416583477567125e-3 0.19906021780991036425e-5 -0.75580140299436494248e-8 -0.19576060961919820491e-9 0.42613011928888888890e-11 -0.41539443304115604377e-13 ], +[ 0.55908401930063918964e0 -0.54413711036826877753e-2 -0.10788661102511914628e-3 0.19229663322982839331e-5 -0.92714731195118129616e-8 -0.14807038677197394186e-9 0.36920870298666666666e-11 -0.39603726688419162617e-13 ], +[ 0.54778496152925675315e0 -0.58501497933213396670e-2 -0.96582314317855227421e-4 0.18434405235069270228e-5 -0.10541580254317078711e-7 -0.10702303407788943498e-9 0.31563175582222222222e-11 -0.36829748079110481422e-13 ], +[ 0.53571290831682823999e0 -0.62147030670760791791e-2 -0.85782497917111760790e-4 0.17553116363443470478e-5 -0.11432547349815541084e-7 -0.72157091369041330520e-10 0.26630811607111111111e-11 -0.33578660425893164084e-13 ], +[ 0.52295422962048434978e0 -0.65371404367776320720e-2 -0.75530164941473343780e-4 0.16613725797181276790e-5 -0.12003521296598910761e-7 -0.42929753689181106171e-10 0.22170894940444444444e-11 -0.30117697501065110505e-13 ], +[ 0.50959092577577886140e0 -0.68197117603118591766e-2 -0.65852936198953623307e-4 0.15639654113906716939e-5 -0.12308007991056524902e-7 -0.18761997536910939570e-10 0.18198628922666666667e-11 -0.26638355362285200932e-13 ], +[ 0.49570040481823167970e0 -0.70647509397614398066e-2 -0.56765617728962588218e-4 0.14650274449141448497e-5 -0.12393681471984051132e-7 0.92904351801168955424e-12 0.14706755960177777778e-11 -0.23272455351266325318e-13 ], +[ 0.48135536250935238066e0 -0.72746293327402359783e-2 -0.48272489495730030780e-4 0.13661377309113939689e-5 -0.12302464447599382189e-7 0.16707760028737074907e-10 0.11672928324444444444e-11 -0.20105801424709924499e-13 ], +[ 0.46662374675511439448e0 -0.74517177649528487002e-2 -0.40369318744279128718e-4 0.12685621118898535407e-5 -0.12070791463315156250e-7 0.29105507892605823871e-10 0.90653314645333333334e-12 -0.17189503312102982646e-13 ], +[ 0.45156879030168268778e0 -0.75983560650033817497e-2 -0.33045110380705139759e-4 0.11732956732035040896e-5 -0.11729986947158201869e-7 0.38611905704166441308e-10 0.68468768305777777779e-12 -0.14549134330396754575e-13 ], +[ 0.43624909769330896904e0 -0.77168291040309554679e-2 -0.26283612321339907756e-4 0.10811018836893550820e-5 -0.11306707563739851552e-7 0.45670446788529607380e-10 0.49782492549333333334e-12 -0.12191983967561779442e-13 ], +[ 0.42071877443548481181e0 -0.78093484015052730097e-2 -0.20064596897224934705e-4 0.99254806680671890766e-6 -0.10823412088884741451e-7 0.50677203326904716247e-10 0.34200547594666666666e-12 -0.10112698698356194618e-13 ], +[ 0.40502758809710844280e0 -0.78780384460872937555e-2 -0.14364940764532853112e-4 0.90803709228265217384e-6 -0.10298832847014466907e-7 0.53981671221969478551e-10 0.21342751381333333333e-12 -0.82975901848387729274e-14 ], +[ 0.38922115269731446690e0 -0.79249269708242064120e-2 -0.91595258799106970453e-5 0.82783535102217576495e-6 -0.97484311059617744437e-8 0.55889029041660225629e-10 0.10851981336888888889e-12 -0.67278553237853459757e-14 ], +[ 0.37334112915460307335e0 -0.79519385109223148791e-2 -0.44219833548840469752e-5 0.75209719038240314732e-6 -0.91848251458553190451e-8 0.56663266668051433844e-10 0.23995894257777777778e-13 -0.53819475285389344313e-14 ], +[ 0.35742543583374223085e0 -0.79608906571527956177e-2 -0.12530071050975781198e-6 0.68088605744900552505e-6 -0.86181844090844164075e-8 0.56530784203816176153e-10 -0.43120012248888888890e-13 -0.42372603392496813810e-14 ], +[ 0.34150846431979618536e0 -0.79534924968773806029e-2 0.37576885610891515813e-5 0.61419263633090524326e-6 -0.80565865409945960125e-8 0.55684175248749269411e-10 -0.95486860764444444445e-13 -0.32712946432984510595e-14 ], +[ 0.32562129649136346824e0 -0.79313448067948884309e-2 0.72539159933545300034e-5 0.55195028297415503083e-6 -0.75063365335570475258e-8 0.54281686749699595941e-10 -0.13545424295111111111e-12 ], +[ 0.30979191977078391864e0 -0.78959416264207333695e-2 0.10389774377677210794e-4 0.49404804463196316464e-6 -0.69722488229411164685e-8 0.52469254655951393842e-10 -0.16507860650666666667e-12 ], +[ 0.29404543811214459904e0 -0.78486728990364155356e-2 0.13190885683106990459e-4 0.44034158861387909694e-6 -0.64578942561562616481e-8 0.50354306498006928984e-10 -0.18614473550222222222e-12 ], +[ 0.27840427686253660515e0 -0.77908279176252742013e-2 0.15681928798708548349e-4 0.39066226205099807573e-6 -0.59658144820660420814e-8 0.48030086420373141763e-10 -0.20018995173333333333e-12 ], +[ 0.26288838011163800908e0 -0.77235993576119469018e-2 0.17886516796198660969e-4 0.34482457073472497720e-6 -0.54977066551955420066e-8 0.45572749379147269213e-10 -0.20852924954666666667e-12 ], +[ 0.24751539954181029717e0 -0.76480877165290370975e-2 0.19827114835033977049e-4 0.30263228619976332110e-6 -0.50545814570120129947e-8 0.43043879374212005966e-10 -0.21228012028444444444e-12 ], +[ 0.23230087411688914593e0 -0.75653060136384041587e-2 0.21524991113020016415e-4 0.26388338542539382413e-6 -0.46368974069671446622e-8 0.40492715758206515307e-10 -0.21238627815111111111e-12 ], +[ 0.21725840021297341931e0 -0.74761846305979730439e-2 0.23000194404129495243e-4 0.22837400135642906796e-6 -0.42446743058417541277e-8 0.37958104071765923728e-10 -0.20963978568888888889e-12 ], +[ 0.20239979200788191491e0 -0.73815761980493466516e-2 0.24271552727631854013e-4 0.19590154043390012843e-6 -0.38775884642456551753e-8 0.35470192372162901168e-10 -0.20470131678222222222e-12 ], +[ 0.18773523211558098962e0 -0.72822604530339834448e-2 0.25356688567841293697e-4 0.16626710297744290016e-6 -0.35350521468015310830e-8 0.33051896213898864306e-10 -0.19811844544000000000e-12 ], +[ 0.17327341258479649442e0 -0.71789490089142761950e-2 0.26272046822383820476e-4 0.13927732375657362345e-6 -0.32162794266956859603e-8 0.30720156036105652035e-10 -0.19034196304000000000e-12 ], +[ 0.15902166648328672043e0 -0.70722899934245504034e-2 0.27032932310132226025e-4 0.11474573347816568279e-6 -0.29203404091754665063e-8 0.28487010262547971859e-10 -0.18174029063111111111e-12 ], +[ 0.14498609036610283865e0 -0.69628725220045029273e-2 0.27653554229160596221e-4 0.92493727167393036470e-7 -0.26462055548683583849e-8 0.26360506250989943739e-10 -0.17261211260444444444e-12 ], +[ 0.13117165798208050667e0 -0.68512309830281084723e-2 0.28147075431133863774e-4 0.72351212437979583441e-7 -0.23927816200314358570e-8 0.24345469651209833155e-10 -0.16319736960000000000e-12 ], +[ 0.11758232561160626306e0 -0.67378491192463392927e-2 0.28525664781722907847e-4 0.54156999310046790024e-7 -0.21589405340123827823e-8 0.22444150951727334619e-10 -0.15368675584000000000e-12 ], +[ 0.10422112945361673560e0 -0.66231638959845581564e-2 0.28800551216363918088e-4 0.37758983397952149613e-7 -0.19435423557038933431e-8 0.20656766125421362458e-10 -0.14422990012444444444e-12 ], +[ 0.91090275493541084785e-1 -0.65075691516115160062e-2 0.28982078385527224867e-4 0.23014165807643012781e-7 -0.17454532910249875958e-8 0.18981946442680092373e-10 -0.13494234691555555556e-12 ], +[ 0.78191222288771379358e-1 -0.63914190297303976434e-2 0.29079759021299682675e-4 0.97885458059415717014e-8 -0.15635596116134296819e-8 0.17417110744051331974e-10 -0.12591151763555555556e-12 ], +[ 0.65524757106147402224e-1 -0.62750311956082444159e-2 0.29102328354323449795e-4 -0.20430838882727954582e-8 -0.13967781903855367270e-8 0.15958771833747057569e-10 -0.11720175765333333333e-12 ], +[ 0.53091065838453612773e-1 -0.61586898417077043662e-2 0.29057796072960100710e-4 -0.12597414620517987536e-7 -0.12440642607426861943e-8 0.14602787128447932137e-10 -0.10885859114666666667e-12 ], +[ 0.40889797115352738582e-1 -0.60426484889413678200e-2 0.28953496450191694606e-4 -0.21982952021823718400e-7 -0.11044169117553026211e-8 0.13344562332430552171e-10 -0.10091231402844444444e-12 ] +] + +# the length of each row in caselist_w_im_y10 +const size_caselist_w_im_y100 = [7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 9 9 8 9 9 9 8 9 9 9 9 9 9 9 9 9 9 8 9 9 9 9 9 9 9 9 8 8 8 8 8 8 8 8 8 8 8 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ] + +""" + w_im_y100(y100,x) + +Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x). + +Uses a look-up table of 100 different Chebyshev polynomials +for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated +with the help of Maple and a little shell script. This allows +the Chebyshev polynomials to be of significantly lower degree (about 1/30) +compared to fitting the whole [0,1] interval with a single polynomial. +""" +function w_im_y100(y100::Float64, x::Float64) + if isnan(y100) || isnan(x) + return NaN + end + theCase = Int64(floor(y100)) + 1 + if theCase <= 97 + t = 2 * y100 - (2 * theCase - 1) + if size_caselist_w_im_y100[theCase] == 7 + return @horner(t, + caselist_w_im_y100[theCase][1], + caselist_w_im_y100[theCase][2], + caselist_w_im_y100[theCase][3], + caselist_w_im_y100[theCase][4], + caselist_w_im_y100[theCase][5], + caselist_w_im_y100[theCase][6], + caselist_w_im_y100[theCase][7] + ) + elseif size_caselist_w_im_y100[theCase] == 8 + return @horner(t, + caselist_w_im_y100[theCase][1], + caselist_w_im_y100[theCase][2], + caselist_w_im_y100[theCase][3], + caselist_w_im_y100[theCase][4], + caselist_w_im_y100[theCase][5], + caselist_w_im_y100[theCase][6], + caselist_w_im_y100[theCase][7], + caselist_w_im_y100[theCase][8] + ) + elseif size_caselist_w_im_y100[theCase] == 9 + return @horner(t, + caselist_w_im_y100[theCase][1], + caselist_w_im_y100[theCase][2], + caselist_w_im_y100[theCase][3], + caselist_w_im_y100[theCase][4], + caselist_w_im_y100[theCase][5], + caselist_w_im_y100[theCase][6], + caselist_w_im_y100[theCase][7], + caselist_w_im_y100[theCase][8], + caselist_w_im_y100[theCase][9] + ) + end + elseif 98 <= theCase <= 101 + x2 = x^2 + return x * (1.1283791670955125739 + - x2 * (0.75225277806367504925 + - x2 * (0.30090111122547001970 + - x2 * (0.085971746064420005629 + - x2 * 0.016931216931216931217)))) + else + return NaN + end +end + +""" + w_im(x) + +Compute a scaled Dawson integral + - w_im(x) = 2*Dawson(x)/sqrt(pi) +equivalent to the imaginary part w(x) for real x. + +Uses methods similar to the erfcx calculation above: continued fractions +for large |x|, a lookup table of Chebyshev polynomials for smaller |x|, +and finally a Taylor expansion for |x|<0.01. +""" +function w_im(x::Float64) + const ispi = 0.56418958354775628694807945156 + if x >= 0 + if x > 45 + if x > 5e7 + return ispi / x + end + return ispi * ((x * x) * (x * x - 4.5) + 2) / (x * ((x * x) * (x * x - 5) + 3.75)) + end + return w_im_y100(100 / (1 + x), x) + else + if x < -45 + if x < -5e7 + return ispi / x + end + return ispi * ((x * x) * (x * x - 4.5) + 2) / (x * ((x * x) * (x * x - 5) + 3.75)) + end + return -w_im_y100(100 / (1 - x), -x) + end +end