erfcx(BigFloat) misses a connection segment

[putianyi889]

Plots.plot(20000:40000,erfcx.(BigFloat(20000):40000), yaxis=:log)

[stevengj]

Rather than a plot, can you just give a single point where it is not giving the expected output?

[putianyi889]

NaN begins from 27282 and ends at 32768 (only counted integer values).

julia> erfcx(BigFloat(27281))
2.068067824378598630288975599374805428229073376310733390485590903274295760840468e-05

julia> erfcx(BigFloat(27282))
NaN

julia> erfcx(BigFloat(32768))
NaN

julia> erfcx(BigFloat(32769))
1.721717425875220708758243853735950836098892035259101544766730697542870688938946e-05

[devmotion]

It seems the heuristic in

SpecialFunctions.jl/src/erf.jl

Line 52 in a1df96d

if x <= (Clong == Int32 ? 0x1p15 : 0x1p30)

is wrong and internal overflow occurs also for values less than or equal to 0x1p15 = 32768.

[kaltsoplyn]

Is there a fix planned for this?

[inkydragon]

SpecialFunctions.jl/src/erf.jl

Line 54 in a1df96d

return exp(x^2)*erfc(x)

julia> ex = big"27281"; (exp(ex^2), erfc(ex))
(1.266065248520218332393804858158913551964455781340246368117884220328223878971691e+323224954, 1.633460697855630920863644577730562244908520023549224867866220531861327684119909e-323224959)

julia> ex = big"27282"; (exp(ex^2), erfc(ex))
(Inf, 0.0)

julia> ex = big"272868"; (exp(ex^2), erfc(ex))
(Inf, 0.0)

exp(ex^2) overflow.

Related julia issue: JuliaLang/julia#40661

---

Maybe we can do something like faddeeva::erfcx(double x)
https://github.com/scipy/xsf/blob/cd55a14bcdaacaae440cc5928ccd8dd22c4cc997/include/xsf/faddeeva.h#L1287-L1294

    if (x > 50) { // continued-fraction expansion is faster
      const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
      if (x > 5e7) // 1-term expansion, important to avoid overflow
        return ispi / x;
      /* 5-term expansion (rely on compiler for CSE), simplified from:
	        ispi / (x+0.5/(x+1/(x+1.5/(x+2/x))))  */
      return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75));
    }

https://dlmf.nist.gov/7.9

julia> x = big"27281"; 1 / sqrt(BigFloat(pi)) *((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75))
2.06806782437859863028897559937480542822907341027272739490356827246976154701339e-05

julia> x = big"27282"; 1 / sqrt(BigFloat(pi)) *((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75))
2.067992020998289267517616942659877979608088875268148689261372732860471272096829e-05

julia> x = big"32768"; 1 / sqrt(BigFloat(pi)) *((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75))
1.721769968521225096297092565194187658031720000298874393245387951092613514709156e-05

julia> x = big"32769"; 1 / sqrt(BigFloat(pi)) *((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75))
1.721717425875220708758243853735950836098892039781421161259973555876921448789855e-05