Skip to content

Add :ulp rounding mode based on prevfloat and nextfloat #744

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
OlivierHnt wants to merge 2 commits into
base: master
Choose a base branch
from
+118 −128
Open

Add :ulp rounding mode based on prevfloat and nextfloat #744

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9a7ae87
Add `:ulp`rounding mode via CORE-MATH
OlivierHnt Feb 20, 2026
2aa7e54
Specialize `_fround` implementations for `acot` and `acoth`
OlivierHnt Feb 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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]
Expand All @@ -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"
Expand Down
4 changes: 2 additions & 2 deletions src/IntervalArithmetic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
234 changes: 111 additions & 123 deletions src/intervals/rounding.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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},
Expand All @@ -39,64 +38,29 @@ 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($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(^), ::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{: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))

_fround(::typeof(^), ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) =
_fround(^, ir, promote(x, n)..., r)

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{: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)
prec = precision(x)
bigx = BigFloat(x; precision = prec)
bigz = BigFloat(; precision = prec)
@ccall Base.MPFR.libmpfr.mpfr_div(
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(inv), ::IntervalRounding, x::AbstractFloat, r::RoundingMode) = _fround(/, one(x), x, r)

#
_fround(::typeof(inv), ::IntervalRounding{:none}, x::T, ::RoundingMode) where {T<:AbstractFloat} = inv(x)

_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
RoundingEmulator.sqrt_down(x)
_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
RoundingEmulator.sqrt_up(x)
function _fround(::typeof(sqrt), ::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)
Expand All @@ -108,68 +72,33 @@ function _fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::AbstractFloat,
return bigz
end

_fround(::typeof(sqrt), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = sqrt(x)

#

function rootn end # defined in arithmetic/power.jl

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
end

_fround(::typeof(rootn), ::IntervalRounding{:none}, x::AbstractFloat, n::Integer, ::RoundingMode) = x^(1//n)

#
_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::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} =
RoundingEmulator.sqrt_up(x)

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(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(atan), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = atan(x, y)
_fround(::typeof(sqrt), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = sqrt(x)

#
# 1-argument functions

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)
bigx = BigFloat(x; precision = prec)
bigz = BigFloat(; precision = prec)
@ccall Base.MPFR.libmpfr.$(Symbol(:mpfr_, f))(
bigz::Ref{BigFloat},
bigx::Ref{BigFloat},
r::Base.MPFR.MPFRRoundingMode
)::Int32
return bigz
end
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

_fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
end
end
mpfr_f = Symbol(:mpfr_, f)
coremath_f = Symbol(:cr_, f)

for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh))
@eval begin
function _fround(::typeof($f), ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Down})
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)
Expand All @@ -182,7 +111,7 @@ for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh))
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})
@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)
Expand All @@ -195,23 +124,82 @@ for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh))
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

_fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
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

@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

# CRlibm functions
# 2-argument functions

_fround(::typeof(^), ::IntervalRounding, x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation

_fround(::typeof(^), ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) =
_fround(^, ir, promote(x, n)..., r)

for (f, fname) ∈ ((:^, :pow), (:atan, :atan2))

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)
mpfr_f = Symbol(:mpfr_, fname)
coremath_f = Symbol(:cr_, fname)

_fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x)
@eval begin
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.$mpfr_f(
bigz::Ref{BigFloat},
bigx::Ref{BigFloat},
bigy::Ref{BigFloat},
r::Base.MPFR.MPFRRoundingMode
)::Int32
return bigz
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::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y)
end
end

function rootn end # defined in arithmetic/power.jl

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)




Expand Down
Loading