Skip to content

Commit fb1224c

Browse files
committed
Merge pull request #1 from lbenet/master
Scattered corrections and improvements, and tests for python and julia versions
2 parents c1cf5f3 + 5eacc6a commit fb1224c

File tree

6 files changed

+403
-289
lines changed

6 files changed

+403
-289
lines changed

src/julia/Intervals.jl

Lines changed: 39 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,17 @@
11
## Intervals.jl
22
##
3-
## Last modification: 2014.04.10
3+
## Last modification: 2014.06.09
44
## Luis Benet and David P. Sanders, UNAM
55
##
66
## Julia module for handling Interval arithmetics
77
##
88

99
module Intervals
1010

11-
import
12-
Base: in, zero, one, show,
13-
sqrt, exp, log, sin, cos, tan,
14-
union, intersect, isempty
11+
import Base: in, zero, one, show,
12+
sqrt, exp, log, sin, cos, tan, inv,
13+
union, intersect, isempty,
14+
convert, promote_rule
1515

1616
export
1717
@roundingDown, @roundingUp, @roundingDirect, @rounded_interval,
@@ -40,7 +40,7 @@ macro roundingUp(expr)
4040
end
4141

4242
## Interval constructor
43-
immutable Interval <: Real
43+
immutable Interval <: Number
4444
lo :: Real
4545
hi :: Real
4646

@@ -58,7 +58,6 @@ immutable Interval <: Real
5858
elseif T == Rational{Int64} || T == Rational{BigInt} || T == Rational{Int32}
5959
lo = @roundingDown( BigFloat(a) )
6060
else
61-
#lo = BigFloat(string(a), RoundDown)
6261
lo = @roundingDown( BigFloat(string(a)) )
6362
end
6463

@@ -68,40 +67,34 @@ immutable Interval <: Real
6867
elseif T == Rational{Int64} || T == Rational{BigInt} || T == Rational{Int32}
6968
hi = @roundingUp( BigFloat(b) )
7069
else
71-
#hi = BigFloat(string(b), RoundUp)
7270
hi = @roundingUp( BigFloat(string(b)) )
7371
end
7472

7573
new(lo, hi)
7674
end
7775
end
7876
Interval(a::Interval) = a
79-
Interval(a::Tuple) = Interval(a[1], a[2])
80-
#?Interval(a::BigFloat) = Interval( @roundingDown("$a"), @roundingUp("$a") )
81-
Interval(b::Bool) = Interval(int(b)) ## Included because a warning: Bool <: Real --> true
77+
Interval(a::Tuple) = Interval(a...)
8278
Interval(a::Real) = Interval(a, a)
8379

80+
# Convertion and promotion
81+
convert(::Type{Interval},x::Real) = Interval(x)
82+
promote_rule{A<:Real}(::Type{Interval}, ::Type{A}) = Interval
83+
8484
## Equalities and neg-equalities
8585
==(a::Interval, b::Interval) = a.lo == b.lo && a.hi == b.hi
86-
==(a::Interval, x::Real) = ==(a, Interval(x))
87-
==(x::Real, a::Interval) = ==(a, Interval(x))
8886
!=(a::Interval, b::Interval) = a.lo != b.lo || a.hi != b.hi
89-
!=(a::Interval, x::Real) = !=(a, Interval(x))
90-
!=(x::Real, a::Interval) = !=(a, Interval(x))
9187

9288
## Inclusion/containment functions
9389
in(a::Interval, b::Interval) = b.lo <= a.lo && a.hi <= b.hi
94-
in(x::Rational, a::Interval) = a.lo <= BigFloat(x) <= a.hi
95-
in(x::BigFloat, a::Interval) = a.lo <= x <= a.hi
96-
in(x::Real, a::Interval) = a.lo <= BigFloat(string(x)) <= a.hi
90+
in(x::Real, a::Interval) = in(promote(x,a)...)
91+
9792
isinside(a::Interval, b::Interval) = b.lo < a.lo && a.hi < b.hi
98-
isinside(x::Rational, a::Interval) = a.lo < BigFloat(x) < a.hi
99-
isinside(x::BigFloat, a::Interval) = a.lo < x < a.hi
100-
isinside(x::Real, a::Interval) = a.lo < BigFloat(string(x)) < a.hi
93+
isinside(x::Real, a::Interval) = isinside(promote(x,a)...)
10194

10295
## zero and one functions
103-
zero(a::Interval) = Interval(BigFloat("0"))
104-
one(a::Interval) = Interval(BigFloat("1"))
96+
zero(a::Interval) = Interval(zero(BigFloat))
97+
one(a::Interval) = Interval(one(BigFloat))
10598

10699
## Addition
107100
function +(a::Interval, b::Interval)
@@ -112,7 +105,6 @@ function +(a::Interval, b::Interval)
112105
set_rounding(BigFloat, RoundNearest)
113106
Interval( lo, hi )
114107
end
115-
#?+(a::Interval, b::Interval) = @rounded_interval( a.lo+b.lo, a.hi+b.hi )
116108
+(a::Interval) = a
117109

118110
## Substraction
@@ -124,7 +116,6 @@ function -(a::Interval, b::Interval)
124116
set_rounding(BigFloat, RoundNearest)
125117
Interval( lo, hi )
126118
end
127-
#?-(a::Interval, b::Interval) = @rounded_interval( a.lo-b.hi, a.hi-b.lo )
128119
-(a::Interval) = Interval(-a.hi,-a.lo)
129120

130121
## Multiplication
@@ -136,15 +127,6 @@ function *(a::Interval, b::Interval)
136127
set_rounding(BigFloat, RoundNearest)
137128
Interval( lo, hi )
138129
end
139-
# The following is needed because to avoid a silly warning
140-
*(a::Interval, b::Bool) = b ? a : zero(a)
141-
*(b::Bool, a::Interval) = b ? a : zero(a)
142-
143-
## Operations mixing Interval and Real
144-
for op in (:+, :-, :*)
145-
@eval $(op)(a::Interval, x::Real) = $(op)(a, Interval(x))
146-
@eval $(op)(x::Real, a::Interval) = $(op)(Interval(x), a)
147-
end
148130

149131
## Division
150132
function reciprocal(a::Interval)
@@ -161,19 +143,17 @@ function reciprocal(a::Interval)
161143
set_rounding(BigFloat, RoundNearest)
162144
Interval( lo, hi )
163145
end
146+
inv(a::Interval) = reciprocal(a)
164147
/(a::Interval, b::Interval) = a*reciprocal(b)
165-
/(a::Interval, x::Real) = a*reciprocal( Interval(x) )
166-
/(x::Real, a::Interval) = Interval(x)*reciprocal(a)
167148

168-
## Some scalar functions on intervals
149+
## Some scalar functions on intervals; no direct rounding used
169150
diam(a::Interval) = a.hi - a.lo
170151
mid(a::Interval) = BigFloat(0.5) * (a.hi + a.lo)
171152
mag(a::Interval) = max( abs(a.lo), abs(a.hi) )
172153
mig(a::Interval) = in(BigFloat(0.0),a) ? BigFloat(0.0) : min( abs(a.lo), abs(a.hi) )
173154

174155
## Intersection
175156
isempty(a::Interval, b::Interval) = a.hi < b.lo || b.hi < a.lo
176-
#?isempty(a::Interval, x::Real) = isempty(a, Interval(x))
177157
function intersect(a::Interval, b::Interval)
178158
if isempty(a,b)
179159
warn("Intersection is empty")
@@ -198,18 +178,19 @@ function hull(a::Interval, b::Interval)
198178
set_rounding(BigFloat, RoundNearest)
199179
Interval( lo, hi )
200180
end
201-
hull(a::Interval, x::Real) = hull(a, Interval(x))
202-
hull(x::Real, a::Interval) = hull(a, Interval(x))
203181

204182
## union
205183
function union(a::Interval, b::Interval)
206184
isempty(a,b) && warn("Empty intersection; union is computed as hull")
207185
hull(a,b)
208186
end
209-
union(a::Interval, x::Real) = union(a, Interval(x))
210-
union(x::Real, a::Interval) = union(Interval(x), a)
211187

212-
#----- From here on, NEEDS TESTING ------
188+
## Extending operators that mix Interval and Real
189+
for fn in (:intersect, :hull, :union)
190+
@eval $(fn)(a::Interval, x::Real) = $(fn)(promote(a,x)...)
191+
@eval $(fn)(x::Real, a::Interval) = $(fn)(promote(x,a)...)
192+
end
193+
213194
## Int power
214195
function ^(a::Interval, n::Integer)
215196
n < zero(n) && return reciprocal( a^(-n) )
@@ -270,26 +251,28 @@ function ^(a::Interval, x::Interval)
270251
# Is x a thin interval?
271252
diam( x ) < eps( mid(x) ) && return a^(x.lo)
272253
z = zero(BigFloat)
273-
z > a.hi && error("Undefined operation; Interval is strictly negative and power is not an integer")
254+
z > a.hi && error("Undefined operation;\n",
255+
"Interval is strictly negative and power is not an integer")
274256
#
275257
domainPow = Interval(z, inf(BigFloat))
276258
aRestricted = intersect(a, domainPow)
277259
set_rounding(BigFloat, RoundDown)
278-
lolo = aRestricted.lo^x.lo
279-
lohi = aRestricted.lo^x.hi
260+
lolo = aRestricted.lo^(x.lo)
261+
lohi = aRestricted.lo^(x.hi)
280262
set_rounding(BigFloat, RoundUp)
281-
hilo = aRestricted.hi^x.lo
282-
hihi = aRestricted.hi^x.hi
263+
hilo = aRestricted.hi^(x.lo)
264+
hihi = aRestricted.hi^(x.hi)
283265
set_rounding(BigFloat, RoundNearest)
284-
lo = min( lolo, lohi)
285-
hi = min( hilo, hihi)
266+
lo = min( lolo, lohi )
267+
hi = min( hilo, hihi )
286268
Interval( lo, hi )
287269
end
288270

289271
## sqrt
290272
function sqrt(a::Interval)
291273
z = zero(BigFloat)
292-
z > a.hi && error("Undefined operation; Interval is strictly negative and power is not an integer")
274+
z > a.hi && error("Undefined operation;\n",
275+
"Interval is strictly negative and power is not an integer")
293276
#
294277
domainSqrt = Interval(z, inf(BigFloat))
295278
aRestricted = intersect(a, domainSqrt)
@@ -326,6 +309,7 @@ function log(a::Interval)
326309
Interval( lo, hi )
327310
end
328311

312+
#----- From here on, NEEDS TESTING ------
329313
## sin
330314
function sin(a::Interval)
331315
piHalf = pi*BigFloat("0.5")
@@ -469,13 +453,13 @@ function tan(a::Interval)
469453
hi = tan( a.hi )
470454
set_rounding(BigFloat, RoundNearest)
471455

472-
if (loHalf > hiHalf) || ( loHalf == hiHalf && lo_modpi < hi_modpi)
456+
if (loHalf > hiHalf) || ( loHalf == hiHalf && loModpi <= hiModpi)
473457
return Interval( lo, hi )
474458
end
475459

476-
disjoint2 = Interval( tan_xlo, BigFloat(Inf) )
477-
disjoint1 = Interval( BigFloat(-Inf), tan_xhi )
478-
info(string("\n The resulting interval is disjoint:\n", disjoint1, disjoint2,
460+
disjoint2 = Interval( lo, BigFloat(Inf) )
461+
disjoint1 = Interval( BigFloat(-Inf), hi )
462+
info(string("The resulting interval is disjoint:\n", disjoint1, "\n", disjoint2,
479463
"\n The hull of the disjoint subintervals is considered:\n", domainTan))
480464
return domainTan
481465
end

src/julia/Readme.md

Lines changed: 54 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,6 @@
33
Interval arithmetics on the real line with julia
44

55
#### Related Work ####
6-
76
- [MPFI.jl][1]: Wrapper of [MPFI library][http://perso.ens-lyon.fr/nathalie.revol/software.html] (multiple precision interval arithmetic library based on MPFR) for julia.
87

98
#### References ####
@@ -14,11 +13,11 @@ Interval arithmetics on the real line with julia
1413
- Luis Benet, Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México (UNAM)
1514
- David P. Sanders, Facultad de Ciencias, Universidad Nacional Autónoma de México (UNAM)
1615

17-
## Design and implementation ##
1816

17+
## Design and implementation ##
1918
`Intervals.jl` is a julia module to perform *validated numerics*, i.e. *rigorous* computations with finite-precision floating-point arithmetic.
2019

21-
The basic constructor is `Interval`, which takes pairs of `Real`-type values as arguments, or alternatively, single values. The idea of including the latter is to produce *thin intervals*, i.e., intervals which consist of a single real value. Yet, if the number is not exactly representable in base 2, rounding-off errors appear due to the finite precision floating-point arithmetic. Incorporating *directed rounding* may thus create intervals whose *diameter* or *width* is strictly positive. Directed rounding is incorporated for any `Real`-type values *except* for `BigFloat`, since this type includes rounding. Direct access to the `RoundDown` and `RoundUp` modes is obtained through the macros `@roundingDown` and `@roundingUp`.
20+
The basic constructor is `Interval`, which takes pairs of `Real`-type values as arguments, or alternatively, single values for *thin* intervals. Yet, if the number is not exactly representable in base 2, rounding-off errors appear due to the finite precision floating-point arithmetic. Incorporating *directed rounding* may thus create intervals whose *diameter* or *width* is non-zero, despite of *a priori* being thin intervals. Directed rounding is incorporated for any `Real`-type values *except* for `BigFloat`, since this type implements already rounding. Direct access to `RoundDown` and `RoundUp` modes is obtained through the macros `@roundingDown` and `@roundingUp`.
2221

2322
The following provides some examples.
2423
```julia
@@ -38,19 +37,21 @@ Interval (constructor with 5 methods)
3837
julia> Interval(1//10)
3938
[9.9999999999999992e-02, 1.0000000000000001e-01] with 53 bits of precision
4039

40+
julia> Interval(0.1,0.1)
41+
[9.9999999999999992e-02, 1.0000000000000001e-01] with 53 bits of precision
42+
4143
julia> Interval(BigFloat(0.1))
4244
[1.0000000000000001e-01, 1.0000000000000001e-01] with 53 bits of precision
4345

4446
julia> b1 = Interval(BigFloat("0.1"))
4547
[1.0000000000000001e-01, 1.0000000000000001e-01] with 53 bits of precision
4648

47-
julia> b2 = Interval( @roundingDown(BigFloat("0.1")), @roundingDown(BigFloat("0.1")) )
48-
[9.9999999999999992e-02, 9.9999999999999992e-02] with 53 bits of precision
49+
julia> b2 = Interval( @roundingDown(BigFloat("0.1")), @roundingUp(BigFloat("0.1")) )
50+
[9.9999999999999992e-02, 1.0000000000000001e-01] with 53 bits of precision
4951

5052
```
5153

52-
Note that all possibilities above yield the same results *except* `Interval(BigFloat(0.1))`, which yields a thin interval, that noticeably *does not* contain the real value 0.1; a way to obtain a non-thin interval using `BigFloat` of strings is illustrated in the last line. It is worth emphasizing that the behavior implemented by the definition of `Interval` above is *not the same* as that implemented by [MPFI][1], which in all cases above
53-
would yield thin intervals; this is the main design difference.
54+
Note that all possibilities above yield the same results *except* `Interval(BigFloat(0.1))` and the definition of `b1`, which yield thin intervals that noticeably *do not* contain the real value 0.1; a way to obtain a non-thin interval using `BigFloat` of strings is illustrated by the definition of `b2`. It is worth emphasizing that the behavior implemented by the definition of `Interval` above is *not the same* as that implemented by [MPFI][1], which in all cases above would yield thin intervals. (Using `MPFI`, the behaviour above is obtained with `Interval("0.1")`, which is not defined in our implementation.) This is the main design difference.
5455

5556
The limits of the interval are accessed with the fields `lo` and `hi`. These fields are of `BigFloat` type. The diameter is obtained using `diam(a)`; note that for non-exactly representable numbers in base 2 (i.e., whose *binary* expansion is infinite) the diameter corresponds to the local `eps` value.
5657
```julia
@@ -60,7 +61,7 @@ julia> a.lo
6061
julia> typeof(a.lo)
6162
BigFloat (constructor with 14 methods)
6263

63-
julia> a.hi.prec ## this displays the precision of the BigFloat number
64+
julia> a.hi.prec ## this displays the precision of the BigFloat number
6465
53
6566

6667
julia> diam(a) == eps(0.1)
@@ -74,10 +75,53 @@ true
7475

7576
```
7677

77-
The basic operations among intervals (`+`, `-`, `*`, `/`, `^`) are included following the usual definitions, and include consisting *directed rounding*. In particular, in the case `^` we have distinguished for non-integer powers, if they define a thin interval or not. Some simple functions (`exp`, `log`, `sin`, `cos`, `tan`) have been incorporated and more will come with time (and patience).
78+
## Operations and functions ##
79+
The basic operations among intervals (`+`, `-`, `*`, `/`, `^`) have been implemented following the usual definitions; they are implemented with consisting *directed rounding* mode. That is, the resulting lower (right) edge is systematically rounded down, while the upper (left) one is rounded up. In particular, the implementation for `^` distinguishes whether non-integer powers define a thin interval or not.
80+
81+
```julia
82+
julia> a = Interval(1.1)
83+
[1.0999999999999999e+00, 1.1000000000000001e+00] with 53 bits of precision
84+
85+
julia> b = Interval(1,2)
86+
[1e+00, 2e+00] with 53 bits of precision
87+
88+
julia> a + b
89+
[2.0999999999999996e+00, 3.1000000000000001e+00] with 53 bits of precision
90+
91+
julia> a - b
92+
[-9.0000000000000013e-01, 1.0000000000000009e-01] with 53 bits of precision
93+
94+
julia> a * b
95+
[1.0999999999999999e+00, 2.2000000000000002e+00] with 53 bits of precision
96+
97+
julia> a / b
98+
[5.4999999999999993e-01, 1.1000000000000001e+00] with 53 bits of precision
99+
100+
julia> b/Interval(-1,1)
101+
WARNING:
102+
Interval in denominator contains 0.
103+
[-inf, inf] with 53 bits of precision
104+
105+
julia> b/Interval(0,2)
106+
[5e-01, inf] with 53 bits of precision
107+
108+
julia> c = Interval(-0.1, 0.2)
109+
[-1.0000000000000001e-01, 2.0000000000000001e-01] with 53 bits of precision
110+
111+
julia> c*c
112+
[-2.0000000000000004e-02, 4.0000000000000008e-02] with 53 bits of precision
113+
114+
julia> c^2
115+
[0e+00, 4.0000000000000008e-02] with 53 bits of precision
116+
117+
```
118+
119+
The last results shows that `^` is implemented independently from the product; this is related to the *dependency problem*.
120+
121+
Some simple functions (`exp`, `log`, `sin`, `cos`, `tan`) have also been implemented, and more will come with time (and patience).
78122

79123

80-
### Acknowledgements ###
124+
## Acknowledgements ##
81125
This project was developed in a masters' course in the postgraduate programs in Physics and in Mathematics at UNAM during the second half of 2013. We thank the participants of the course for putting up with the half-baked material and contributing energy and ideas.
82126

83127
Financial support is acknowledged from DGAPA-UNAM PAPIME grants PE-105911 and PE-107114. LB acknowledges support through a *Cátedra Moshinsky* (2013).

0 commit comments

Comments
 (0)