From b8bd2960708106d585d113b039ecef527288dba1 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 1 Oct 2021 01:14:02 +0200 Subject: [PATCH 1/2] Avoid StackOverflowError with erf* functions --- src/erf.jl | 86 ++++++++++++++++++++++++++++------------------------- test/erf.jl | 24 +++++++++++---- 2 files changed, 64 insertions(+), 46 deletions(-) diff --git a/src/erf.jl b/src/erf.jl index 4cac6bc9..f86aee61 100644 --- a/src/erf.jl +++ b/src/erf.jl @@ -4,43 +4,49 @@ using Base.Math: @horner using Base.MPFR: ROUNDING_MODE for f in (:erf, :erfc) + internalf = Symbol(:_, f) + libopenlibmf = QuoteNode(f) + libopenlibmf0 = QuoteNode(Symbol(f, :f)) + openspecfunf = QuoteNode(Symbol(:Faddeeva_, f)) + mpfrf = QuoteNode(Symbol(:mpfr_, f)) @eval begin - ($f)(x::Float64) = ccall(($(string(f)),libopenlibm), Float64, (Float64,), x) - ($f)(x::Float32) = ccall(($(string(f,"f")),libopenlibm), Float32, (Float32,), x) - ($f)(x::Real) = ($f)(float(x)) - ($f)(a::Float16) = Float16($f(Float32(a))) - ($f)(a::Complex{Float16}) = Complex{Float16}($f(Complex{Float32}(a))) - function ($f)(x::BigFloat) + $f(x::Number) = $internalf(float(x)) + + $internalf(x::Float64) = ccall(($libopenlibmf, libopenlibm), Float64, (Float64,), x) + $internalf(x::Float32) = ccall(($libopenlibmf0, libopenlibm), Float32, (Float32,), x) + $internalf(x::Float16) = Float16($internalf(Float32(x))) + + $internalf(z::Complex{Float64}) = Complex{Float64}(ccall(($openspecfunf, libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) + $internalf(z::Complex{Float32}) = Complex{Float32}(ccall(($openspecfunf, libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex{Float64}(z), Float64(eps(Float32)))) + $internalf(z::Complex{Float16}) = Complex{Float16}($internalf(Complex{Float32}(z))) + + function $internalf(x::BigFloat) z = BigFloat() - ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[]) + ccall(($mpfrf, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[]) return z end - ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x)) end end -for f in (:erf, :erfc, :erfcx, :erfi, :Dawson) - fname = (f === :Dawson) ? :dawson : f +for f in (:erfcx, :erfi, :dawson) + internalf = Symbol(:_, f) + openspecfunfsym = Symbol(:Faddeeva_, f === :dawson ? :Dawson : f) + openspecfunfF64 = QuoteNode(Symbol(openspecfunfsym, :_re)) + openspecfunfCF64 = QuoteNode(openspecfunfsym) @eval begin - ($fname)(z::Complex{Float64}) = Complex{Float64}(ccall(($(string("Faddeeva_",f)),libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) - ($fname)(z::Complex{Float32}) = Complex{Float32}(ccall(($(string("Faddeeva_",f)),libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex{Float64}(z), Float64(eps(Float32)))) + $f(x::Number) = $internalf(float(x)) - ($fname)(z::Complex) = ($fname)(float(z)) - ($fname)(z::Complex{<:AbstractFloat}) = throw(MethodError($fname,(z,))) - end -end - -for f in (:erfcx, :erfi, :Dawson) - fname = (f === :Dawson) ? :dawson : f - @eval begin - ($fname)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),libopenspecfun), Float64, (Float64,), x) - ($fname)(x::Float32) = Float32(ccall(($(string("Faddeeva_",f,"_re")),libopenspecfun), Float64, (Float64,), Float64(x))) + $internalf(x::Float64) = ccall(($openspecfunfF64, libopenspecfun), Float64, (Float64,), x) + $internalf(x::Float32) = Float32($internalf(Float64(x))) + $internalf(x::Float16) = Float16($internalf(Float64(x))) - ($fname)(x::Real) = ($fname)(float(x)) - ($fname)(x::AbstractFloat) = throw(MethodError($fname,(x,))) + $internalf(z::Complex{Float64}) = Complex{Float64}(ccall(($openspecfunfCF64, libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) + $internalf(z::Complex{Float32}) = Complex{Float32}(ccall(($openspecfunfCF64, libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex{Float64}(z), Float64(eps(Float32)))) + $internalf(z::Complex{Float16}) = Complex{Float16}($internalf(Complex{Float32}(z))) end end + @doc raw""" erf(x) @@ -204,7 +210,9 @@ Using the rational approximants tabulated in: > combined with Newton iterations for `BigFloat`. """ -function erfinv(x::Float64) +erfinv(x::Real) = _erfinv(float(x)) + +function _erfinv(x::Float64) a = abs(x) if a >= 1.0 if x == 1.0 @@ -272,7 +280,7 @@ function erfinv(x::Float64) end end -function erfinv(x::Float32) +function _erfinv(x::Float32) a = abs(x) if a >= 1.0f0 if x == 1.0f0 @@ -315,7 +323,6 @@ function erfinv(x::Float32) end end -erfinv(x::Union{Integer,Rational}) = erfinv(float(x)) @doc raw""" erfcinv(x) @@ -341,7 +348,9 @@ Using the rational approximants tabulated in: > combined with Newton iterations for `BigFloat`. """ -function erfcinv(y::Float64) +erfcinv(x::Real) = _erfcinv(float(x)) + +function _erfcinv(y::Float64) if y > 0.0625 return erfinv(1.0 - y) elseif y <= 0.0 @@ -393,7 +402,7 @@ function erfcinv(y::Float64) end end -function erfcinv(y::Float32) +function _erfcinv(y::Float32) if y > 0.0625f0 return erfinv(1.0f0 - y) elseif y <= 0.0f0 @@ -415,7 +424,7 @@ function erfcinv(y::Float32) end end -function erfinv(y::BigFloat) +function _erfinv(y::BigFloat) xfloat = erfinv(Float64(y)) if isfinite(xfloat) x = BigFloat(xfloat) @@ -435,7 +444,7 @@ function erfinv(y::BigFloat) return x end -function erfcinv(y::BigFloat) +function _erfcinv(y::BigFloat) yfloat = Float64(y) xfloat = erfcinv(yfloat) if isfinite(xfloat) @@ -461,7 +470,6 @@ function erfcinv(y::BigFloat) return x end -erfcinv(x::Union{Integer,Rational}) = 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 @@ -511,7 +519,9 @@ See also: [`erfcx(x)`](@ref erfcx). Based on the [`erfc(x)`](@ref erfc) and [`erfcx(x)`](@ref erfcx) functions. Currently only implemented for `Float32`, `Float64`, and `BigFloat`. """ -function logerfc(x::Union{Float32, Float64, BigFloat}) +logerfc(x::Real) = _logerfc(float(x)) + +function _logerfc(x::Union{Float32, Float64, BigFloat}) # Don't include Float16 in the Union, otherwise logerfc would currently work for x <= 0.0, but not x > 0.0 if x > 0.0 return log(erfcx(x)) - x^2 @@ -520,9 +530,6 @@ function logerfc(x::Union{Float32, Float64, BigFloat}) end end -logerfc(x::Real) = logerfc(float(x)) -logerfc(x::AbstractFloat) = throw(MethodError(logerfc, x)) - @doc raw""" logerfcx(x) @@ -543,7 +550,9 @@ See also: [`erfcx(x)`](@ref erfcx). Based on the [`erfc(x)`](@ref erfc) and [`erfcx(x)`](@ref erfcx) functions. Currently only implemented for `Float32`, `Float64`, and `BigFloat`. """ -function logerfcx(x::Union{Float32, Float64, BigFloat}) +logerfcx(x::Real) = _logerfcx(float(x)) + +function _logerfcx(x::Union{Float32, Float64, BigFloat}) # Don't include Float16 in the Union, otherwise logerfc would currently work for x <= 0.0, but not x > 0.0 if x < 0.0 return log(erfc(x)) + x^2 @@ -552,9 +561,6 @@ function logerfcx(x::Union{Float32, Float64, BigFloat}) end end -logerfcx(x::Real) = logerfcx(float(x)) -logerfcx(x::AbstractFloat) = throw(MethodError(logerfcx, x)) - @doc raw""" logerf(x, y) diff --git a/test/erf.jl b/test/erf.jl index 91f94cf4..16c77846 100644 --- a/test/erf.jl +++ b/test/erf.jl @@ -8,7 +8,7 @@ @test erfc(Float32(1)) ≈ 0.15729920705028513066 rtol=2*eps(Float32) @test erfc(Float64(1)) ≈ 0.15729920705028513066 rtol=2*eps(Float64) - @test_throws MethodError erfcx(Float16(1)) + @test erfcx(Float16(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float16) @test erfcx(Float32(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float32) @test erfcx(Float64(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float64) @@ -38,7 +38,7 @@ @test logerfcx(Float32(1000)) ≈ -7.48012072190621214066734919080 rtol=2eps(Float32) @test logerfcx(Float64(1000)) ≈ -7.48012072190621214066734919080 rtol=2eps(Float64) - @test_throws MethodError erfi(Float16(1)) + @test erfi(Float16(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float16) @test erfi(Float32(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float32) @test erfi(Float64(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float64) @@ -52,7 +52,7 @@ @test erfcinv(Float32(0.15729920705028513066)) ≈ 1 rtol=2*eps(Float32) @test erfcinv(Float64(0.15729920705028513066)) ≈ 1 rtol=2*eps(Float64) - @test_throws MethodError dawson(Float16(1)) + @test dawson(Float16(1)) ≈ 0.53807950691276841914 rtol=2*eps(Float16) @test dawson(Float32(1)) ≈ 0.53807950691276841914 rtol=2*eps(Float32) @test dawson(Float64(1)) ≈ 0.53807950691276841914 rtol=2*eps(Float64) end @@ -66,11 +66,11 @@ @test erfc(ComplexF32(1+2im)) ≈ 1.5366435657785650340+5.0491437034470346695im @test erfc(ComplexF64(1+2im)) ≈ 1.5366435657785650340+5.0491437034470346695im - @test_throws MethodError erfcx(ComplexF16(1)) + @test erfcx(ComplexF16(1+2im)) ≈ 0.14023958136627794370-0.22221344017989910261im @test erfcx(ComplexF32(1+2im)) ≈ 0.14023958136627794370-0.22221344017989910261im @test erfcx(ComplexF64(1+2im)) ≈ 0.14023958136627794370-0.22221344017989910261im - @test_throws MethodError erfi(ComplexF16(1)) + @test erfi(ComplexF16(1+2im)) ≈ -0.011259006028815025076+1.0036063427256517509im @test erfi(ComplexF32(1+2im)) ≈ -0.011259006028815025076+1.0036063427256517509im @test erfi(ComplexF64(1+2im)) ≈ -0.011259006028815025076+1.0036063427256517509im @@ -78,7 +78,7 @@ @test_throws MethodError erfcinv(Complex(1)) - @test_throws MethodError dawson(ComplexF16(1)) + @test dawson(ComplexF16(1+2im)) ≈ -13.388927316482919244-11.828715103889593303im @test dawson(ComplexF32(1+2im)) ≈ -13.388927316482919244-11.828715103889593303im @test dawson(ComplexF64(1+2im)) ≈ -13.388927316482919244-11.828715103889593303im end @@ -116,6 +116,18 @@ end end + @testset "Other float types" begin + struct NotAFloat <: AbstractFloat end + + @test_throws MethodError erf(NotAFloat()) + @test_throws MethodError erfc(NotAFloat()) + @test_throws MethodError erfcx(NotAFloat()) + @test_throws MethodError erfi(NotAFloat()) + @test_throws MethodError erfinv(NotAFloat()) + @test_throws MethodError erfcinv(NotAFloat()) + @test_throws MethodError dawson(NotAFloat()) + end + @testset "inverse" begin for elty in [Float32,Float64] for x in exp10.(range(-200, stop=-0.01, length=50)) From d0867829181fdcad7ba5466e6818df1820613ecb Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 1 Oct 2021 01:14:43 +0200 Subject: [PATCH 2/2] Rearrange BigFloat implementations --- src/erf.jl | 95 ++++++++++++++++++++++++++---------------------------- 1 file changed, 46 insertions(+), 49 deletions(-) diff --git a/src/erf.jl b/src/erf.jl index f86aee61..1eb080f6 100644 --- a/src/erf.jl +++ b/src/erf.jl @@ -46,6 +46,33 @@ for f in (:erfcx, :erfi, :dawson) end end +# 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π) + end +end @doc raw""" erf(x) @@ -323,6 +350,25 @@ function _erfinv(x::Float32) end end +function _erfinv(y::BigFloat) + xfloat = erfinv(Float64(y)) + if isfinite(xfloat) + x = BigFloat(xfloat) + else + # Float64 overflowed, use asymptotic estimate instead + # from erfc(x) ≈ exp(-x²)/x√π ≈ y ⟹ -log(yπ) ≈ x² + log(x) ≈ x² + x = copysign(sqrt(-log((1-abs(y))*sqrtπ)), y) + isfinite(x) || return x + end + sqrtπhalf = sqrtπ * big(0.5) + tol = 2eps(abs(x)) + while true # Newton iterations + Δx = sqrtπhalf * (erf(x) - y) * exp(x^2) + x -= Δx + abs(Δx) < tol && break + end + return x +end @doc raw""" erfcinv(x) @@ -424,26 +470,6 @@ function _erfcinv(y::Float32) end end -function _erfinv(y::BigFloat) - xfloat = erfinv(Float64(y)) - if isfinite(xfloat) - x = BigFloat(xfloat) - else - # Float64 overflowed, use asymptotic estimate instead - # from erfc(x) ≈ exp(-x²)/x√π ≈ y ⟹ -log(yπ) ≈ x² + log(x) ≈ x² - x = copysign(sqrt(-log((1-abs(y))*sqrtπ)), y) - isfinite(x) || return x - end - sqrtπhalf = sqrtπ * big(0.5) - tol = 2eps(abs(x)) - while true # Newton iterations - Δx = sqrtπhalf * (erf(x) - y) * exp(x^2) - x -= Δx - abs(Δx) < tol && break - end - return x -end - function _erfcinv(y::BigFloat) yfloat = Float64(y) xfloat = erfcinv(yfloat) @@ -470,35 +496,6 @@ function _erfcinv(y::BigFloat) return x end - -# 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π) - end -end - @doc raw""" logerfc(x)