Skip to content

Commit af8be29

Browse files
OlivierHntOlivierHnt
committed
Use CRlibm.jl and introduce configure function for rounding
1 parent a4dfb0d commit af8be29

File tree

4 files changed

+93
-184
lines changed

4 files changed

+93
-184
lines changed

‎Project.toml

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git"
44
version = "0.22.29"
55

66
[deps]
7-
CRlibm_jll = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
7+
CRlibm = "96374032-68de-5a5b-8d9e-752f78720389"
88
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
99
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
1010
OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af"
@@ -23,7 +23,7 @@ IntervalArithmeticIntervalSetsExt = "IntervalSets"
2323
IntervalArithmeticRecipesBaseExt = "RecipesBase"
2424

2525
[compat]
26-
CRlibm_jll = "1"
26+
CRlibm = "1.0.2"
2727
DiffRules = "1"
2828
ForwardDiff = "0.10, 1"
2929
IntervalSets = "0.7"

‎src/IntervalArithmetic.jl

Lines changed: 47 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,8 @@ Learn more: https://github.com/JuliaIntervals/IntervalArithmetic.jl
77
"""
88
module IntervalArithmetic
99

10-
import CRlibm_jll
11-
import RoundingEmulator
12-
import Base.MPFR
10+
import RoundingEmulator, CRlibm, Base.MPFR
11+
1312
using MacroTools: MacroTools, prewalk, postwalk, @capture
1413

1514
#
@@ -49,4 +48,49 @@ end
4948

5049
include("matmul.jl")
5150

51+
#
52+
53+
function configure_rounding(rounding::Symbol)
54+
rounding ∈ (:correct, :none) || return throw(ArgumentError("only the rounding mode `:correct` and `:none` are available"))
55+
56+
for f ∈ (:add, :sub, :mul, :div)
57+
f_round = Symbol(:_, f, :_round)
58+
@eval $f_round(x, y, r) = $f_round(IntervalRounding{$(QuoteNode(rounding))}(), x, y, r)
59+
end
60+
61+
@eval _pow_round(x, y, r) = _pow_round(IntervalRounding{$(QuoteNode(rounding))}(), x, y, r)
62+
63+
@eval _inv_round(x, r) = _inv_round(IntervalRounding{$(QuoteNode(rounding))}(), x, r)
64+
65+
@eval _sqrt_round(x, r) = _sqrt_round(IntervalRounding{$(QuoteNode(rounding))}(), x, r)
66+
67+
@eval _rootn_round(x, n, r) = _rootn_round(IntervalRounding{$(QuoteNode(rounding))}(), x, n, r)
68+
69+
@eval _atan_round(x, y, r) = _atan_round(IntervalRounding{$(QuoteNode(rounding))}(), x, y, r)
70+
71+
for f ∈ [:cbrt, :exp2, :exp10, :cot, :sec, :csc, :tanh, :coth, :sech, :csch, :asinh, :acosh, :atanh]
72+
f_round = Symbol(:_, f, :_round)
73+
@eval $f_round(x, r) = $f_round(IntervalRounding{$(QuoteNode(rounding))}(), x, r)
74+
end
75+
76+
for f ∈ (:acot, :acoth)
77+
f_round = Symbol(:_, f, :_round)
78+
@eval $f_round(x, r) = $f_round(IntervalRounding{$(QuoteNode(rounding))}(), x, r)
79+
end
80+
81+
for f ∈ CRlibm.functions
82+
f_round = Symbol(:_, f, :_round)
83+
@eval $f_round(x, r) = $f_round(IntervalRounding{$(QuoteNode(rounding))}(), x, r)
84+
end
85+
86+
return rounding
87+
end
88+
89+
function configure(; rounding::Symbol=:correct)
90+
configure_rounding(rounding)
91+
return rounding
92+
end
93+
94+
configure()
95+
5296
end

‎src/intervals/rounding.jl

Lines changed: 29 additions & 164 deletions
Original file line numberDiff line numberDiff line change
@@ -4,41 +4,26 @@
44
Interval rounding type.
55
66
Available rounding types:
7-
- `:fast`: rounding via `prevfloat` and `nextfloat`. Currently only the functions
8-
`+`, `-`, `*`, `/`, `inv` and `sqrt` are supported.
9-
- `:tight`: rounding via [RoundingEmulator.jl](https://github.com/matsueushi/RoundingEmulator.jl).
10-
- `:slow`: rounding via `BigFloat`. The other rounding modes default to `:slow`
11-
whenever a faster algorithm is not available.
7+
- `:correct`: rounding via [RoundingEmulator.jl](https://github.com/JuliaIntervals/RoundingEmulator.jl)
8+
and [CRlibm.jl](https://github.com/JuliaInterval/IntervalArithmetic.jl).
129
- `:none`: no rounding (non-rigorous numerics).
1310
"""
1411
struct IntervalRounding{T} end
1512

16-
interval_rounding() = IntervalRounding{:tight}()
17-
1813
#
1914

2015
for (f, fname) ∈ ((:+, :add), (:-, :sub), (:*, :mul), (:/, :div))
2116
g = Symbol(:_, fname, :_round)
2217
mpfr_f = Symbol(:mpfr_, fname)
2318

2419
@eval begin
25-
$g(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} = $g(interval_rounding(), x, y, r)
26-
$g(x::T, y::T, ::RoundingMode) where {T<:Rational} = $f(x, y) # exact operation
27-
28-
$g(::IntervalRounding, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} =
29-
$g(IntervalRounding{:slow}(), x, y, r)
20+
$g(::IntervalRounding, x::T, y::T, ::RoundingMode) where {T<:Rational} = $f(x, y) # exact operation
3021

31-
$g(::IntervalRounding{:fast}, x::T, y::T, ::RoundingMode{:Down}) where {T<:AbstractFloat} =
32-
prevfloat($f(x, y))
33-
$g(::IntervalRounding{:fast}, x::T, y::T, ::RoundingMode{:Up}) where {T<:AbstractFloat} =
34-
nextfloat($f(x, y))
35-
36-
$g(::IntervalRounding{:tight}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} =
22+
$g(::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} =
3723
RoundingEmulator.$(Symbol(fname, :_down))(x, y)
38-
$g(::IntervalRounding{:tight}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} =
24+
$g(::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} =
3925
RoundingEmulator.$(Symbol(fname, :_up))(x, y)
40-
41-
function $g(::IntervalRounding{:slow}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
26+
function $g(::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
4227
prec = max(precision(x), precision(y))
4328
bigx = BigFloat(x; precision = prec)
4429
bigy = BigFloat(y; precision = prec)
@@ -58,20 +43,13 @@ end
5843

5944
#
6045

61-
_pow_round(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} = _pow_round(interval_rounding(), x, y, r)
62-
_pow_round(x::AbstractFloat, n::Integer, r::RoundingMode) = _pow_round(promote(x, n)..., r)
63-
_pow_round(x::T, y::T, ::RoundingMode) where {T<:Rational} = ^(x, y) # exact operation
64-
_pow_round(x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation
46+
_pow_round(ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) =
47+
_pow_round(ir, promote(x, n)..., r)
6548

66-
_pow_round(::IntervalRounding, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} =
67-
_pow_round(IntervalRounding{:slow}(), x, y, r)
49+
_pow_round(::IntervalRounding, x::T, y::T, ::RoundingMode) where {T<:Rational} = ^(x, y) # exact operation
50+
_pow_round(::IntervalRounding, x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation
6851

69-
# _pow_round(::IntervalRounding{:fast}, x::T, y::T, ::RoundingMode{:Down}) where {T<:AbstractFloat} =
70-
# prevfloat(^(x, y))
71-
# _pow_round(::IntervalRounding{:fast}, x::T, y::T, ::RoundingMode{:Up}) where {T<:AbstractFloat} =
72-
# nextfloat(^(x, y))
73-
74-
function _pow_round(::IntervalRounding{:slow}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
52+
function _pow_round(::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
7553
prec = max(precision(x), precision(y))
7654
bigx = BigFloat(x; precision = prec)
7755
bigy = BigFloat(y; precision = prec)
@@ -89,23 +67,13 @@ _pow_round(::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:Abst
8967

9068
#
9169

92-
_inv_round(x::AbstractFloat, r::RoundingMode) = _inv_round(interval_rounding(), x, r)
93-
_inv_round(x::Rational, ::RoundingMode) = inv(x) # exact operation
94-
95-
_inv_round(::IntervalRounding, x::AbstractFloat, r::RoundingMode) =
96-
_inv_round(IntervalRounding{:slow}(), x, r)
70+
_inv_round(::IntervalRounding, x::Rational, ::RoundingMode) = inv(x) # exact operation
9771

98-
_inv_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Down}) =
99-
prevfloat(inv(x))
100-
_inv_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Up}) =
101-
nextfloat(inv(x))
102-
103-
_inv_round(::IntervalRounding{:tight}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
72+
_inv_round(::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
10473
RoundingEmulator.div_down(one(x), x)
105-
_inv_round(::IntervalRounding{:tight}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
74+
_inv_round(::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
10675
RoundingEmulator.div_up(one(x), x)
107-
108-
function _inv_round(::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode)
76+
function _inv_round(::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode)
10977
prec = precision(x)
11078
bigx = BigFloat(x; precision = prec)
11179
bigz = BigFloat(; precision = prec)
@@ -122,22 +90,11 @@ _inv_round(::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = inv(x)
12290

12391
#
12492

125-
_sqrt_round(x::AbstractFloat, r::RoundingMode) = _sqrt_round(interval_rounding(), x, r)
126-
127-
_sqrt_round(::IntervalRounding, x::AbstractFloat, r::RoundingMode) =
128-
_sqrt_round(IntervalRounding{:slow}(), x, r)
129-
130-
_sqrt_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Down}) =
131-
prevfloat(sqrt(x))
132-
_sqrt_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Up}) =
133-
nextfloat(sqrt(x))
134-
135-
_sqrt_round(::IntervalRounding{:tight}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
93+
_sqrt_round(::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
13694
RoundingEmulator.sqrt_down(x)
137-
_sqrt_round(::IntervalRounding{:tight}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
95+
_sqrt_round(::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
13896
RoundingEmulator.sqrt_up(x)
139-
140-
function _sqrt_round(::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode)
97+
function _sqrt_round(::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode)
14198
prec = precision(x)
14299
bigx = BigFloat(x; precision = prec)
143100
bigz = BigFloat(; precision = prec)
@@ -153,17 +110,7 @@ _sqrt_round(::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = sqrt(
153110

154111
#
155112

156-
_rootn_round(x::AbstractFloat, n::Integer, r::RoundingMode) = _rootn_round(interval_rounding(), x, n, r)
157-
158-
_rootn_round(::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) =
159-
_rootn_round(IntervalRounding{:slow}(), x, n, r)
160-
161-
# _rootn_round(::IntervalRounding{:fast}, x::AbstractFloat, n::Integer, ::RoundingMode{:Down}) =
162-
# prevfloat(x^(1//n))
163-
# _rootn_round(::IntervalRounding{:fast}, x::AbstractFloat, n::Integer, ::RoundingMode{:Up}) =
164-
# nextfloat(x^(1//n))
165-
166-
function _rootn_round(::IntervalRounding{:slow}, x::AbstractFloat, n::Integer, r::RoundingMode)
113+
function _rootn_round(::IntervalRounding{:correct}, x::AbstractFloat, n::Integer, r::RoundingMode)
167114
prec = precision(x)
168115
bigx = BigFloat(x; precision = prec)
169116
bigz = BigFloat(; precision = prec)
@@ -180,17 +127,7 @@ _rootn_round(::IntervalRounding{:none}, x::AbstractFloat, n::Integer, ::Rounding
180127

181128
#
182129

183-
_atan_round(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} = _atan_round(interval_rounding(), x, y, r)
184-
185-
_atan_round(::IntervalRounding, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} =
186-
_atan_round(IntervalRounding{:slow}(), x, y, r)
187-
188-
# _atan_round(::IntervalRounding{:fast}, x::T, y::T, ::RoundingMode{:Down}) where {T<:AbstractFloat} =
189-
# prevfloat(atan(x, y))
190-
# _atan_round(::IntervalRounding{:fast}, x::T, y::T, ::RoundingMode{:Up}) where {T<:AbstractFloat} =
191-
# nextfloat(atan(x, y))
192-
193-
function _atan_round(::IntervalRounding{:slow}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
130+
function _atan_round(::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
194131
prec = max(precision(x), precision(y))
195132
bigx = BigFloat(x; precision = prec)
196133
bigy = BigFloat(y; precision = prec)
@@ -213,16 +150,7 @@ for f ∈ [:cbrt, :exp2, :exp10, :cot, :sec, :csc, :tanh, :coth, :sech, :csch, :
213150
mpfr_f = Symbol(:mpfr_, f)
214151

215152
@eval begin
216-
$f_round(x::AbstractFloat, r::RoundingMode) = $f_round(interval_rounding(), x, r)
217-
218-
$f_round(::IntervalRounding, x::AbstractFloat, r::RoundingMode) = $f_round(IntervalRounding{:slow}(), x, r)
219-
220-
# $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Down}) =
221-
# prevfloat($f(x))
222-
# $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Up}) =
223-
# nextfloat($f(x))
224-
225-
function $f_round(::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode)
153+
function $f_round(::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode)
226154
prec = precision(x)
227155
bigx = BigFloat(x; precision = prec)
228156
bigz = BigFloat(; precision = prec)
@@ -243,16 +171,7 @@ for (f, g) ∈ [(:acot, :atan), (:acoth, :atanh)]
243171
g_round = Symbol(:_, g, :_round)
244172

245173
@eval begin
246-
$f_round(x::AbstractFloat, r::RoundingMode) = $f_round(interval_rounding(), x, r)
247-
248-
$f_round(::IntervalRounding, x::AbstractFloat, r::RoundingMode) = $f_round(IntervalRounding{:slow}(), x, r)
249-
250-
# $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Down}) =
251-
# prevfloat($f(x))
252-
# $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Up}) =
253-
# nextfloat($f(x))
254-
255-
function $f_round(ir::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode{:Down})
174+
function $f_round(ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Down})
256175
prec = precision(x)
257176
bigx = BigFloat(x; precision = prec + 10)
258177
bigz = BigFloat(; precision = prec + 10)
@@ -265,7 +184,7 @@ for (f, g) ∈ [(:acot, :atan), (:acoth, :atanh)]
265184
bigw = $g_round(ir, bigz, r)
266185
return BigFloat(bigw, r; precision = prec)
267186
end
268-
function $f_round(ir::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode{:Up})
187+
function $f_round(ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Up})
269188
prec = precision(x)
270189
bigx = BigFloat(x; precision = prec + 32)
271190
bigz = BigFloat(; precision = prec + 32)
@@ -285,67 +204,13 @@ end
285204

286205
# CRlibm functions
287206

288-
for f ∈ [:exp, :expm1, :log, :log1p, :log2, :log10, :sin, :cos, :tan, :asin, :acos, :atan, :sinh, :cosh, :sinpi, :cospi]
289-
if isdefined(Base, f)
290-
f_round = Symbol(:_, f, :_round)
291-
crlibm_f_d = string(f, "_rd")
292-
crlibm_f_u = string(f, "_ru")
293-
mpfr_f = Symbol(:mpfr_, f)
294-
295-
if Int == Int32 # issues with CRlibm for 32 bit systems, use MPFR (only available since Julia v1.10)
296-
if VERSION ≥ v"1.10" || f ∉ (:sinpi, :cospi)
297-
@eval $f_round(x::AbstractFloat, r::RoundingMode) = $f_round(interval_rounding(), x, r)
298-
299-
@eval $f_round(::IntervalRounding, x::AbstractFloat, r::RoundingMode) = $f_round(IntervalRounding{:slow}(), x, r)
300-
301-
# @eval $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Down}) =
302-
# prevfloat($f(x))
303-
# @eval $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Up}) =
304-
# nextfloat($f(x))
305-
306-
@eval function $f_round(::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode)
307-
prec = precision(x)
308-
bigx = BigFloat(x; precision = prec)
309-
bigz = BigFloat(; precision = prec)
310-
@ccall Base.MPFR.libmpfr.$mpfr_f(
311-
bigz::Ref{BigFloat},
312-
bigx::Ref{BigFloat},
313-
r::Base.MPFR.MPFRRoundingMode
314-
)::Int32
315-
return bigz
316-
end
317-
318-
@eval $f_round(::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
319-
end
320-
else
321-
@eval $f_round(x::AbstractFloat, r::RoundingMode) = $f_round(interval_rounding(), x, r)
322-
323-
@eval $f_round(::IntervalRounding, x::AbstractFloat, r::RoundingMode) = $f_round(IntervalRounding{:slow}(), x, r)
324-
325-
# @eval $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Down}) =
326-
# prevfloat($f(x))
327-
# @eval $f_round(::IntervalRounding{:fast}, x::AbstractFloat, ::RoundingMode{:Up}) =
328-
# nextfloat($f(x))
329-
330-
@eval $f_round(::IntervalRounding{:tight}, x::Float16, r::RoundingMode) = Float16($f_round(Float64(x), r), r)
331-
@eval $f_round(::IntervalRounding{:tight}, x::Float32, r::RoundingMode) = Float32($f_round(Float64(x), r), r)
332-
@eval $f_round(::IntervalRounding{:tight}, x::Float64, r::RoundingMode{:Down}) = ccall(($crlibm_f_d, CRlibm_jll.libcrlibm), Float64, (Float64,), x)
333-
@eval $f_round(::IntervalRounding{:tight}, x::Float64, r::RoundingMode{:Up}) = ccall(($crlibm_f_u, CRlibm_jll.libcrlibm), Float64, (Float64,), x)
334-
335-
@eval function $f_round(::IntervalRounding{:slow}, x::AbstractFloat, r::RoundingMode)
336-
prec = precision(x)
337-
bigx = BigFloat(x; precision = prec)
338-
bigz = BigFloat(; precision = prec)
339-
@ccall Base.MPFR.libmpfr.$mpfr_f(
340-
bigz::Ref{BigFloat},
341-
bigx::Ref{BigFloat},
342-
r::Base.MPFR.MPFRRoundingMode
343-
)::Int32
344-
return bigz
345-
end
207+
for f ∈ CRlibm.functions
208+
f_round = Symbol(:_, f, :_round)
346209

347-
@eval $f_round(::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
348-
end
210+
@eval begin
211+
$f_round(::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode) = CRlibm.$f(x, r)
212+
213+
$f_round(::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x)
349214
end
350215
end
351216

0 commit comments

Comments
 (0)