From 9a7ae873cb54442ef6bde8a49bfeacbdef176293 Mon Sep 17 00:00:00 2001 From: OlivierHnt <38465572+OlivierHnt@users.noreply.github.com> Date: Fri, 20 Feb 2026 19:31:34 +0800 Subject: [PATCH 1/2] Add `:ulp`rounding mode via CORE-MATH --- Project.toml | 8 +- src/IntervalArithmetic.jl | 4 +- src/intervals/rounding.jl | 226 ++++++++++++++++---------------------- 3 files changed, 99 insertions(+), 139 deletions(-) diff --git a/Project.toml b/Project.toml index 3dfb05e63..1647abe6a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,15 @@ name = "IntervalArithmetic" uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" -repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git" version = "1.0.2" +repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git" [deps] CRlibm = "96374032-68de-5a5b-8d9e-752f78720389" +CoreMath = "b7a15901-be09-4a0e-87d2-2e66b0e09b5a" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705" [weakdeps] @@ -32,14 +33,15 @@ IntervalArithmeticSparseArraysExt = "SparseArrays" [compat] Arblib = "1.3.0" CRlibm = "1.0.2" +CoreMath = "0.1.0" DiffRules = "1" ForwardDiff = "1" IntervalSets = "0.7" LinearAlgebra = "1.10" MacroTools = "0.5" OpenBLASConsistentFPCSR_jll = "0.3.29" -Random = "1.10" Printf = "1.10" +Random = "1.10" RecipesBase = "1" RoundingEmulator = "0.2" SparseArrays = "1.10" diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index d12f81414..d837bec23 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -44,7 +44,7 @@ with decorations and up to 6 significant digits. """ module IntervalArithmetic -import RoundingEmulator, CRlibm, Base.MPFR +import RoundingEmulator, CRlibm, CoreMath, Base.MPFR using MacroTools: MacroTools, prewalk, postwalk, @capture @@ -102,7 +102,7 @@ function configure_flavor(flavor::Symbol) end function configure_rounding(rounding::Symbol) - rounding ∈ (:correct, :none) || return throw(ArgumentError("only the rounding mode `:correct` and `:none` are available")) + rounding ∈ (:correct, :ulp, :none) || return throw(ArgumentError("only the rounding mode `:correct`, `:ulp` and `:none` are available")) @eval default_rounding() = IntervalRounding{$(QuoteNode(rounding))}() return rounding end diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl index 4b8bf295e..17acb0660 100644 --- a/src/intervals/rounding.jl +++ b/src/intervals/rounding.jl @@ -4,8 +4,8 @@ Interval rounding type. Available rounding types: -- `:correct`: rounding via [RoundingEmulator.jl](https://github.com/JuliaIntervals/RoundingEmulator.jl) -and [CRlibm.jl](https://github.com/JuliaInterval/IntervalArithmetic.jl). +- `:correct`: rounding via [RoundingEmulator.jl](https://github.com/JuliaIntervals/RoundingEmulator.jl) and [CRlibm.jl](https://github.com/JuliaInterval/IntervalArithmetic.jl); fallback to MPFR. +- `:ulp`: rounding via [CoreMath.jl](https://github.com/JuliaMath/CoreMath.jl) (default rounding to nearest) with `prevfloat` and `nextfloat`; fallback to MPFR. - `:none`: no rounding (non-rigorous numerics). """ struct IntervalRounding{T} end @@ -15,22 +15,21 @@ struct IntervalRounding{T} end _fround(f::Function, x, y, r) = _fround(f, default_rounding(), x, y, r) _fround(f::Function, x, r) = _fround(f, default_rounding(), x, r) -# +# +, -, *, /, inv, sqrt for (f, fname) ∈ ((:+, :add), (:-, :sub), (:*, :mul), (:/, :div)) + + mpfr_f = Symbol(:mpfr_, fname) + @eval begin _fround(::typeof($f), ::IntervalRounding, x::T, y::T, ::RoundingMode) where {T<:Rational} = $f(x, y) # exact operation - _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} = - RoundingEmulator.$(Symbol(fname, :_down))(x, y) - _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} = - RoundingEmulator.$(Symbol(fname, :_up))(x, y) - function _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} + function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} prec = max(precision(x), precision(y)) bigx = BigFloat(x; precision = prec) bigy = BigFloat(y; precision = prec) bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.$(Symbol(:mpfr_, fname))( + @ccall Base.MPFR.libmpfr.$mpfr_f( bigz::Ref{BigFloat}, bigx::Ref{BigFloat}, bigy::Ref{BigFloat}, @@ -39,179 +38,138 @@ for (f, fname) ∈ ((:+, :add), (:-, :sub), (:*, :mul), (:/, :div)) return bigz end - _fround(::typeof($f), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y) - end -end - -# - -_fround(::typeof(^), ::IntervalRounding, x::T, y::T, ::RoundingMode) where {T<:Rational} = ^(x, y) # exact operation -_fround(::typeof(^), ::IntervalRounding, x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation + _fround(::typeof($f), ir::IntervalRounding{:correct}, x::Float16, y::Float16, r::RoundingMode) = + Float16(_fround($f, ir, Float64(x), Float64(y), r), r) + _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} = + RoundingEmulator.$(Symbol(fname, :_down))(x, y) + _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} = + RoundingEmulator.$(Symbol(fname, :_up))(x, y) -_fround(::typeof(^), ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) = - _fround(^, ir, promote(x, n)..., r) + _fround(::typeof($f), ir::IntervalRounding{:ulp}, x::T, y::T, r::RoundingMode) where {T<:Union{Float16,Float32}} = + T(_fround($f, ir, Float64(x), Float64(y), r), r) + _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Float64, y::Float64, ::RoundingMode{:Down}) = prevfloat($f(x, y)) + _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Float64, y::Float64, ::RoundingMode{:Up}) = nextfloat($f(x, y)) -function _fround(::typeof(^), ::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} - prec = max(precision(x), precision(y)) - bigx = BigFloat(x; precision = prec) - bigy = BigFloat(y; precision = prec) - bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.mpfr_pow( - bigz::Ref{BigFloat}, - bigx::Ref{BigFloat}, - bigy::Ref{BigFloat}, - r::Base.MPFR.MPFRRoundingMode - )::Int32 - return bigz + _fround(::typeof($f), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y) + end end -_fround(::typeof(^), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = ^(x, y) +_fround(::typeof(inv), ::IntervalRounding, x::Rational, ::RoundingMode) = inv(x) # exact operation -# +_fround(::typeof(inv), ::IntervalRounding, x::AbstractFloat, r::RoundingMode) = _fround(/, one(x), x, r) -_fround(::typeof(inv), ::IntervalRounding, x::Rational, ::RoundingMode) = inv(x) # exact operation +_fround(::typeof(inv), ::IntervalRounding{:none}, x::T, ::RoundingMode) where {T<:AbstractFloat} = inv(x) -_fround(::typeof(inv), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) = - RoundingEmulator.div_down(one(x), x) -_fround(::typeof(inv), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) = - RoundingEmulator.div_up(one(x), x) -function _fround(::typeof(inv), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode) +function _fround(::typeof(sqrt), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::T, r::RoundingMode) where {T<:AbstractFloat} prec = precision(x) bigx = BigFloat(x; precision = prec) bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.mpfr_div( + @ccall Base.MPFR.libmpfr.mpfr_sqrt( bigz::Ref{BigFloat}, - one(bigx)::Ref{BigFloat}, bigx::Ref{BigFloat}, r::Base.MPFR.MPFRRoundingMode )::Int32 return bigz end -_fround(::typeof(inv), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = inv(x) - -# - -_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) = +_fround(::typeof(sqrt), ir::IntervalRounding{:correct}, x::Float16, r::RoundingMode) = + Float16(_fround(sqrt, ir, Float64(x), r), r) +_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} = RoundingEmulator.sqrt_down(x) -_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) = +_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} = RoundingEmulator.sqrt_up(x) -function _fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode) - prec = precision(x) - bigx = BigFloat(x; precision = prec) - bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.mpfr_sqrt( - bigz::Ref{BigFloat}, - bigx::Ref{BigFloat}, - r::Base.MPFR.MPFRRoundingMode - )::Int32 - return bigz -end + +_fround(::typeof(sqrt), ir::IntervalRounding{:ulp}, x::T, r::RoundingMode) where {T<:Union{Float16,Float32}} = + T(_fround(sqrt, ir, Float64(x), r), r) +_fround(::typeof(sqrt), ::IntervalRounding{:ulp}, x::Float64, ::RoundingMode{:Down}) = prevfloat(sqrt(x)) +_fround(::typeof(sqrt), ::IntervalRounding{:ulp}, x::Float64, ::RoundingMode{:Up}) = nextfloat(sqrt(x)) _fround(::typeof(sqrt), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = sqrt(x) -# +# 1-argument functions + +for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1, # exponential + :log, :log2, :log10, :log1p, # logarithmic + :sin, :sinpi, :cos, :cospi, :tan, :cot, :sec, :csc, :asin, :acos, :atan, :acot, # trigonometric + :sinh, :tanh, :asinh, :cosh, :coth, :sech, :csch, :acosh, :atanh, :acoth) # hyperbolic + + mpfr_f = Symbol(:mpfr_, f) + coremath_f = Symbol(:cr_, f) + + @eval function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode) + prec = precision(x) + bigx = BigFloat(x; precision = prec) + bigz = BigFloat(; precision = prec) + @ccall Base.MPFR.libmpfr.$mpfr_f( + bigz::Ref{BigFloat}, + bigx::Ref{BigFloat}, + r::Base.MPFR.MPFRRoundingMode + )::Int32 + return bigz + end -function rootn end # defined in arithmetic/power.jl + if f ∈ CRlibm.functions + @eval _fround(::typeof($f), ::IntervalRounding{:correct}, x::Float16, r::RoundingMode) = Float16(_fround($f, Float64(x), r), r) + @eval _fround(::typeof($f), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, r::RoundingMode) = CRlibm.$f(x, r) + end -function _fround(::typeof(rootn), ::IntervalRounding{:correct}, x::AbstractFloat, n::Integer, r::RoundingMode) - prec = precision(x) - bigx = BigFloat(x; precision = prec) - bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.mpfr_rootn_ui( - bigz::Ref{BigFloat}, - bigx::Ref{BigFloat}, - n::Culong, - r::Base.MPFR.MPFRRoundingMode - )::Int32 - return bigz + @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Down}) = prevfloat(CoreMath.$coremath_f(x)) + @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Up}) = nextfloat(CoreMath.$coremath_f(x)) + + @eval _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x) end -_fround(::typeof(rootn), ::IntervalRounding{:none}, x::AbstractFloat, n::Integer, ::RoundingMode) = x^(1//n) +# 2-argument functions -# +_fround(::typeof(^), ::IntervalRounding, x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation -function _fround(::typeof(atan), ::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} - prec = max(precision(x), precision(y)) - bigx = BigFloat(x; precision = prec) - bigy = BigFloat(y; precision = prec) - bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.mpfr_atan2( - bigz::Ref{BigFloat}, - bigx::Ref{BigFloat}, - bigy::Ref{BigFloat}, - r::Base.MPFR.MPFRRoundingMode - )::Int32 - return bigz -end +_fround(::typeof(^), ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) = + _fround(^, ir, promote(x, n)..., r) -_fround(::typeof(atan), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = atan(x, y) +for (f, fname) ∈ ((:^, :pow), (:atan, :atan2)) -# + mpfr_f = Symbol(:mpfr_, fname) + coremath_f = Symbol(:cr_, fname) -for f ∈ (:cbrt, :exp2, :exp10, :cot, :sec, :csc, :tanh, :coth, :sech, :csch, :asinh, :acosh, :atanh) @eval begin - function _fround(::typeof($f), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode) - prec = precision(x) + function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} + prec = max(precision(x), precision(y)) bigx = BigFloat(x; precision = prec) + bigy = BigFloat(y; precision = prec) bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.$(Symbol(:mpfr_, f))( + @ccall Base.MPFR.libmpfr.$mpfr_f( bigz::Ref{BigFloat}, bigx::Ref{BigFloat}, + bigy::Ref{BigFloat}, r::Base.MPFR.MPFRRoundingMode )::Int32 return bigz end - _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x) - end -end - -for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh)) - @eval begin - function _fround(::typeof($f), ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Down}) - prec = precision(x) - bigx = BigFloat(x; precision = prec + 10) - bigz = BigFloat(; precision = prec + 10) - @ccall Base.MPFR.libmpfr.mpfr_div( - bigz::Ref{BigFloat}, - one(bigx)::Ref{BigFloat}, - bigx::Ref{BigFloat}, - RoundUp::Base.MPFR.MPFRRoundingMode - )::Int32 - bigw = _fround($g, ir, bigz, r) - return BigFloat(bigw, r; precision = prec) - end - function _fround(::typeof($f), ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Up}) - prec = precision(x) - bigx = BigFloat(x; precision = prec + 32) - bigz = BigFloat(; precision = prec + 32) - @ccall Base.MPFR.libmpfr.mpfr_div( - bigz::Ref{BigFloat}, - one(bigx)::Ref{BigFloat}, - bigx::Ref{BigFloat}, - RoundDown::Base.MPFR.MPFRRoundingMode - )::Int32 - bigw = _fround($g, ir, bigz, r) - return BigFloat(bigw, r; precision = prec) - end + _fround(::typeof($f), ::IntervalRounding{:ulp}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float16,Float32,Float64}} = prevfloat(CoreMath.$coremath_f(x, y)) + _fround(::typeof($f), ::IntervalRounding{:ulp}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float16,Float32,Float64}} = nextfloat(CoreMath.$coremath_f(x, y)) - _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x) + _fround(::typeof($f), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y) end end -# CRlibm functions - -for f ∈ CRlibm.functions - if f !== :atanpi # currently not defined in IntervalArithmetic - @eval begin - _fround(::typeof($f), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode) = CRlibm.$f(x, r) +function rootn end # defined in arithmetic/power.jl - _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x) - end - end +function _fround(::typeof(rootn), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, n::Integer, r::RoundingMode) + prec = precision(x) + bigx = BigFloat(x; precision = prec) + bigz = BigFloat(; precision = prec) + @ccall Base.MPFR.libmpfr.mpfr_rootn_ui( + bigz::Ref{BigFloat}, + bigx::Ref{BigFloat}, + n::Culong, + r::Base.MPFR.MPFRRoundingMode + )::Int32 + return bigz end +_fround(::typeof(rootn), ::IntervalRounding{:none}, x::AbstractFloat, n::Integer, ::RoundingMode) = x^(1//n) + From 2aa7e5463d5e95cd4a465c0729b64f1c119fd6dd Mon Sep 17 00:00:00 2001 From: OlivierHnt <38465572+OlivierHnt@users.noreply.github.com> Date: Sat, 21 Feb 2026 03:28:48 +0800 Subject: [PATCH 2/2] Specialize `_fround` implementations for `acot` and `acoth` --- src/intervals/rounding.jl | 50 +++++++++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 10 deletions(-) diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl index 17acb0660..e84142e2a 100644 --- a/src/intervals/rounding.jl +++ b/src/intervals/rounding.jl @@ -96,16 +96,46 @@ for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1, # exponential mpfr_f = Symbol(:mpfr_, f) coremath_f = Symbol(:cr_, f) - @eval function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode) - prec = precision(x) - bigx = BigFloat(x; precision = prec) - bigz = BigFloat(; precision = prec) - @ccall Base.MPFR.libmpfr.$mpfr_f( - bigz::Ref{BigFloat}, - bigx::Ref{BigFloat}, - r::Base.MPFR.MPFRRoundingMode - )::Int32 - return bigz + if f ∈ (:acot, :acoth) + g = f === :acot ? :atan : :atanh + @eval function _fround(::typeof($f), ir::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode{:Down}) + prec = precision(x) + bigx = BigFloat(x; precision = prec + 10) + bigz = BigFloat(; precision = prec + 10) + @ccall Base.MPFR.libmpfr.mpfr_div( + bigz::Ref{BigFloat}, + one(bigx)::Ref{BigFloat}, + bigx::Ref{BigFloat}, + RoundUp::Base.MPFR.MPFRRoundingMode + )::Int32 + bigw = _fround($g, ir, bigz, r) + return BigFloat(bigw, r; precision = prec) + end + @eval function _fround(::typeof($f), ir::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode{:Up}) + prec = precision(x) + bigx = BigFloat(x; precision = prec + 32) + bigz = BigFloat(; precision = prec + 32) + @ccall Base.MPFR.libmpfr.mpfr_div( + bigz::Ref{BigFloat}, + one(bigx)::Ref{BigFloat}, + bigx::Ref{BigFloat}, + RoundDown::Base.MPFR.MPFRRoundingMode + )::Int32 + bigw = _fround($g, ir, bigz, r) + return BigFloat(bigw, r; precision = prec) + end + else + @eval function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode) + prec = precision(x) + bigx = BigFloat(x; precision = prec) + bigz = BigFloat(; precision = prec) + @ccall Base.MPFR.libmpfr.$mpfr_f( + bigz::Ref{BigFloat}, + bigx::Ref{BigFloat}, + r::Base.MPFR.MPFRRoundingMode + )::Int32 + return bigz + end end if f ∈ CRlibm.functions