Skip to content

Commit dce986d

Browse files
KolaruOlivierHnt
andauthored
New macro for guaranteed operations (#627)
Co-authored-by: OlivierHnt <38465572+OlivierHnt@users.noreply.github.com>
1 parent 118e0d0 commit dce986d

File tree

8 files changed

+247
-10
lines changed

8 files changed

+247
-10
lines changed

Project.toml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ version = "0.22.9"
55

66
[deps]
77
CRlibm_jll = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
8+
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
89
RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
910

1011
[weakdeps]

docs/src/manual/api.md

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,10 @@
11
# API
22

3+
```@meta
4+
CollapsedDocStrings = true
5+
```
6+
37
```@autodocs
48
Modules = [IntervalArithmetic]
5-
Order = [:type, :function, :macro, :constant]
9+
Private = false
610
```

docs/src/manual/construction.md

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -54,9 +54,9 @@ x = interval(0.5, 3)
5454
sqrt(x)
5555
```
5656

57-
Both input `x` and output `sqrt(x)` are common intervals since they are closed, bounded, non-empty and that ``\sqrt`` is continuous over ``[1/2, 3]``.
57+
Both input `x` and output `sqrt(x)` are common intervals since they are closed, bounded, non-empty and the square root is continuous over ``[1/2, 3]``.
5858

59-
Observe that these decorations, together with the fact that any element of the interval `sqrt(x)` is also in the interval `x`, imply that the [Schauder Fixed-Point Theorem](https://en.wikipedia.org/wiki/Schauder_fixed-point_theorem) is satisfied. More precisely, this computation proves the existence of a fixed-point of ``\sqrt`` in ``[1/2, 3]`` (in this simple example, ``\sqrt(1) = 1``).
59+
Observe that these decorations, together with the fact that any element of the interval `sqrt(x)` is also in the interval `x`, imply that the [Schauder Fixed-Point Theorem](https://en.wikipedia.org/wiki/Schauder_fixed-point_theorem) is satisfied. More precisely, this computation proves the existence of a fixed-point of the square root in ``[1/2, 3]`` (in this simple example, ``\sqrt(1) = 1``).
6060

6161
##### Defined and continuous
6262

@@ -90,7 +90,7 @@ x = interval(-3.5, 4)
9090
sqrt(x)
9191
```
9292

93-
The negative part of `x` is discarded before evaluating the ``\sqrt`` function since its domain is ``[0, \infty)``. The process of discarding parts of an interval that are not in the domain of a function is called *loose evaluation*. This event has been recorded by degrading the decoration of the resulting interval to `trv`, indicating that nothing is known about the relationship between `x` and `sqrt(x)`.
93+
The negative part of `x` is discarded before evaluating the square root since its domain is ``[0, \infty)``. The process of discarding parts of an interval that are not in the domain of a function is called *loose evaluation*. This event has been recorded by degrading the decoration of the resulting interval to `trv`, indicating that nothing is known about the relationship between `x` and `sqrt(x)`.
9494

9595
In this case, we know why the decoration was reduced to `trv`. Generally, if this were just a single step in a longer calculation, a resulting decoration `trv` shows only that something like this occured at some step.
9696

@@ -133,7 +133,7 @@ These are all examples of ill-formed intervals, resulting in the decoration `ill
133133

134134
### Guarantee
135135

136-
A *guarantee* is yet another label, independent of decorations, and not described by the IEEE Standard 1788-2015 specifications. Its purpose is to accomodate for Julia's extensive conversion and promotion system, while retaining reliability in computations. Specifically, an interval `x` constructed via [`interval`](@ref) satisfies `isguaranteed(x) == true`. However, if a call to `convert(::Type{<:Interval}, ::Real)` occurs, then the resulting interval `x` satisfies `isguaranteed(x) == false`, receiving the 'NG' (not guaranteed) label. For instance, consider the following examples:
136+
A *guarantee* is yet another label, independent of decorations, and not described by the IEEE Standard 1788-2015 specifications. Its purpose is to accomodate for Julia's extensive conversion and promotion system, while retaining reliability in computations. Specifically, an interval `x` constructed via [`interval`](@ref) satisfies `isguaranteed(x) == true`. However, if a call to `convert(::Type{<:Interval}, ::Real)` occurs, then the resulting interval `x` satisfies `isguaranteed(x) == false`, receiving the "NG" (not guaranteed) label. For instance, consider the following examples:
137137

138138
```@repl construction
139139
convert(Interval{Float64}, 1.) # considered "not guaranteed" as this call can be done implicitly

src/IntervalArithmetic.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ module IntervalArithmetic
33
import CRlibm_jll
44
import RoundingEmulator
55
import Base.MPFR
6+
using MacroTools: MacroTools, prewalk, postwalk, @capture
67

78
#
89

src/intervals/construction.jl

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -402,17 +402,19 @@ Interval{BigFloat}(1.0, 3.141592653589793238462643383279502884197169399375105820
402402
```
403403
"""
404404
function interval(::Type{T}, a, b, d::Decoration = com; format::Symbol = :infsup) where {T}
405-
format === :infsup && return _interval_infsup(T, a, b, d)
406-
format === :midpoint && return _interval_midpoint(T, a, b, d)
405+
format === :infsup && return _interval_infsup(T, _value(a), _value(b), d)
406+
format === :midpoint && return _interval_midpoint(T, _value(a), _value(b), d)
407407
return throw(ArgumentError("`format` must be `:infsup` or `:midpoint`"))
408408
end
409-
interval(a, b, d::Decoration = com; format::Symbol = :infsup) = interval(promote_numtype(numtype(a), numtype(b)), a, b, d; format = format)
409+
interval(a, b, d::Decoration = com; format::Symbol = :infsup) = interval(promote_numtype(numtype(_value(a)), numtype(_value(b))), _value(a), _value(b), d; format = format)
410410

411411
function interval(::Type{T}, a, d::Decoration = com; format::Symbol = :infsup) where {T}
412-
(format === :infsup) | (format === :midpoint) && return _interval_infsup(T, a, a, d)
412+
(format === :infsup) | (format === :midpoint) && return _interval_infsup(T, _value(a), _value(a), d)
413413
return throw(ArgumentError("`format` must be `:infsup` or `:midpoint`"))
414414
end
415-
interval(a, d::Decoration = com; format::Symbol = :infsup) = interval(promote_numtype(numtype(a), numtype(a)), a, d; format = format)
415+
interval(a, d::Decoration = com; format::Symbol = :infsup) = interval(promote_numtype(numtype(_value(a)), numtype(_value(a))), _value(a), d; format = format)
416+
417+
_value(a) = a # convenient hook, used in exact_literals.jl
416418

417419
# some useful extra constructor
418420
interval(::Type{T}, a::Tuple, d::Decoration = com; format::Symbol = :infsup) where {T} = interval(T, a..., d; format = format)

src/intervals/exact_literals.jl

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,202 @@
1+
"""
2+
ExactReal{T<:Real} <: Real
3+
4+
Real numbers with the assurance that they precisely correspond to the number
5+
described by their binary form. The purpose is to guarantee that a non interval
6+
number is exact, so that `ExactReal` can be used with `Interval` without
7+
producing the "NG" flag.
8+
9+
!!! danger
10+
By using `ExactReal`, users acknowledge the responsibility of ensuring that
11+
the number they input corresponds to their intended value.
12+
For example, `ExactReal(0.1)` implies that the user knows that ``0.1`` can
13+
not be represented exactly as a binary number, and that they are using a
14+
slightly different number than ``0.1``.
15+
To help identify the binary number, `ExactReal` is displayed without any
16+
rounding.
17+
18+
```julia
19+
julia> ExactReal(0.1)
20+
ExactReal{Float64}(0.1000000000000000055511151231257827021181583404541015625)
21+
```
22+
23+
In case of doubt, [`has_exact_display`](@ref) can be use to check if the
24+
string representation of a `Real` is equal to its binary value.
25+
26+
# Examples
27+
28+
```jldoctest
29+
julia> setdisplay(:full);
30+
31+
julia> 0.5 * interval(1)
32+
Interval{Float64}(0.5, 0.5, com, NG)
33+
34+
julia> ExactReal(0.5) * interval(1)
35+
Interval{Float64}(0.5, 0.5, com)
36+
37+
julia> [1, interval(2)]
38+
2-element Vector{Interval{Float64}}:
39+
Interval{Float64}(1.0, 1.0, com, NG)
40+
Interval{Float64}(2.0, 2.0, com)
41+
42+
julia> [ExactReal(1), interval(2)]
43+
2-element Vector{Interval{Float64}}:
44+
Interval{Float64}(1.0, 1.0, com)
45+
Interval{Float64}(2.0, 2.0, com)
46+
```
47+
48+
See also: [`@exact`](@ref).
49+
"""
50+
struct ExactReal{T<:Real} <: Real
51+
value :: T
52+
53+
ExactReal(value::T) where {T<:Real} = new{T}(value)
54+
end
55+
56+
_value(x::ExactReal) = x.value # hook for interval constructor
57+
58+
# conversion and promotion
59+
60+
Base.convert(::Type{ExactReal{T}}, x::ExactReal{T}) where {T<:Real} = x
61+
62+
Base.promote_rule(::Type{ExactReal{T}}, ::Type{ExactReal{S}}) where {T<:Real,S<:Real} =
63+
throw(ArgumentError("promotion between `ExactReal` is not allowed"))
64+
65+
# to Interval
66+
67+
Base.convert(::Type{Interval{T}}, x::ExactReal) where {T<:NumTypes} =
68+
interval(T, x.value)
69+
70+
Base.promote_rule(::Type{Interval{T}}, ::Type{ExactReal{S}}) where {T<:NumTypes,S<:Real} =
71+
Interval{promote_numtype(T, S)}
72+
Base.promote_rule(::Type{ExactReal{T}}, ::Type{Interval{S}}) where {T<:Real,S<:NumTypes} =
73+
Interval{promote_numtype(T, S)}
74+
75+
# to Real
76+
77+
Base.convert(::Type{T}, x::ExactReal) where {T<:Real} =
78+
convert(T, x.value)
79+
80+
Base.promote_rule(::Type{T}, ::Type{ExactReal{S}}) where {T<:Real,S<:Real} =
81+
promote_type(T, S)
82+
Base.promote_rule(::Type{ExactReal{T}}, ::Type{S}) where {T<:Real,S<:Real} =
83+
promote_type(T, S)
84+
# needed to resolve method ambiguities
85+
Base.promote_rule(::Type{ExactReal{T}}, ::Type{Bool}) where {T<:Real} =
86+
promote_type(T, Bool)
87+
Base.promote_rule(::Type{Bool}, ::Type{ExactReal{T}}) where {T<:Real} =
88+
promote_type(Bool, T)
89+
90+
# display
91+
92+
Base.string(x::ExactReal{T}) where {T<:AbstractFloat} =
93+
Base.Ryu.writefixed(x.value, 2000, false, false, false, UInt8('.'), true)
94+
95+
Base.string(x::ExactReal) = string(x.value)
96+
97+
Base.show(io::IO, ::MIME"text/plain", x::ExactReal{T}) where {T<:AbstractFloat} =
98+
print(io, "ExactReal{$T}($(string(x)))")
99+
100+
# always exact
101+
102+
Base.:-(x::ExactReal) = ExactReal(-x.value)
103+
104+
function Base.:+(x::ExactReal, y::Complex{<:ExactReal})
105+
iszero(real(y).value) && return complex(x, imag(y))
106+
return complex(x + real(y), imag(y))
107+
end
108+
109+
"""
110+
has_exact_display(x::Real)
111+
112+
Determine if the display of `x` up to 2000 decimals is equal to the bitwise
113+
value of `x`. This is famously not true for the float displayed as `0.1`.
114+
"""
115+
has_exact_display(x::Real) = string(x) == string(ExactReal(x))
116+
117+
#
118+
119+
struct ExactIm end
120+
121+
Base.:*(x::ExactReal, ::ExactIm) = complex(ExactReal(zero(x.value)), x)
122+
Base.:*(::ExactIm, x::ExactReal) = complex(ExactReal(zero(x.value)), x)
123+
124+
Base.:*(x::Real, ::ExactIm) = complex(zero(x), x)
125+
Base.:*(::ExactIm, x::Real) = complex(zero(x), x)
126+
127+
Base.:*(x::Complex, ::ExactIm) = complex(-imag(x), real(x))
128+
Base.:*(::ExactIm, x::Complex) = complex(-imag(x), real(x))
129+
130+
#
131+
132+
"""
133+
@exact
134+
135+
Wrap every literal numbers of the expression in an [`ExactReal`](@ref). This
136+
macro allows defining generic functions, seamlessly accepting both `Number` and
137+
[`Interval`](@ref) arguments, without producing the "NG" flag.
138+
139+
!!! danger
140+
By using [`ExactReal`](@ref), users acknowledge the responsibility of
141+
ensuring that the number they input corresponds to their intended value.
142+
For example, `ExactReal(0.1)` implies that the user knows that ``0.1`` can
143+
not be represented exactly as a binary number, and that they are using a
144+
slightly different number than ``0.1``.
145+
To help identify the binary number, `ExactReal` is displayed without any
146+
rounding.
147+
148+
```julia
149+
julia> ExactReal(0.1)
150+
ExactReal{Float64}(0.1000000000000000055511151231257827021181583404541015625)
151+
```
152+
153+
In case of doubt, [`has_exact_display`](@ref) can be use to check if the
154+
string representation of a `Real` is equal to its binary value.
155+
156+
# Examples
157+
158+
```jldoctest
159+
julia> setdisplay(:full);
160+
161+
julia> f(x) = 1.2*x + 0.1
162+
f (generic function with 1 method)
163+
164+
julia> f(interval(1, 2))
165+
Interval{Float64}(1.2999999999999998, 2.5, com, NG)
166+
167+
julia> @exact g(x) = 1.2*x + 0.1
168+
g (generic function with 1 method)
169+
170+
julia> g(interval(1, 2))
171+
Interval{Float64}(1.2999999999999998, 2.5, com)
172+
173+
julia> g(1.4)
174+
1.78
175+
```
176+
177+
See also: [`ExactReal`](@ref).
178+
"""
179+
macro exact(expr)
180+
exact_expr = postwalk(expr) do x
181+
x isa Real && return ExactReal(x)
182+
return x
183+
end
184+
185+
exact_expr = prewalk(exact_expr) do x
186+
if @capture(x, b_ * im) || @capture(x, im * b_)
187+
if b isa ExactReal
188+
return :(complex(ExactReal(zero($b.value)), $b))
189+
end
190+
end
191+
192+
if @capture(x, a_ + b_ * im) || @capture(x, a_ + im * b_) || @capture(x, b_ * im + a_) || @capture(x, im * b_ + a_)
193+
if a isa ExactReal && b isa ExactReal
194+
return :(complex($a, $b))
195+
end
196+
end
197+
198+
return x
199+
end
200+
201+
return esc(exact_expr)
202+
end

src/intervals/intervals.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,8 @@ include("construction.jl")
44
Interval, interval, isguaranteed, @interval, @I_str
55
include("parsing.jl")
66
include("real_interface.jl")
7+
include("exact_literals.jl")
8+
export ExactReal, @exact, has_exact_display
79

810
# Rounding
911
include("rounding.jl")

test/interval_tests/exact_literals.jl

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
@testset "Exact literals" begin
2+
x = @exact 0.5
3+
@test (2 * x) isa Float64
4+
Y = interval(2) * x
5+
@test isguaranteed(Y)
6+
@test in_interval(1, Y)
7+
8+
@exact function f(x)
9+
return x^2 - 2x + 1
10+
end
11+
12+
Z = f(interval(1))
13+
@test f(1.22) isa Real
14+
@test isguaranteed(Z)
15+
@test in_interval(0, Z)
16+
17+
@test_throws MethodError convert(ExactReal{Float64}, 2)
18+
19+
@test has_exact_display(0.5)
20+
@test !has_exact_display(0.1)
21+
22+
@test (@exact 2im) isa Complex{<:ExactReal}
23+
@test (@exact 1.2 + 3.4im) isa Complex{<:ExactReal}
24+
@test_throws ArgumentError (@exact 1.2 + 3im)
25+
end

0 commit comments

Comments
 (0)