3 (* port of interval.hl,
5 This file gives a simple implementation of interval arithmetic,
6 together with the basic arithmetic operations on intervals.
8 It has been incompletely implemented.
10 For now, I am not implementing directed roundings.
11 However, McLaughlin implemented directed rounding several years ago:
12 See http://perso.ens-lyon.fr/nathalie.revol/mpfi.html
13 ~/Library/McLaughlinOCAML/ocaml/src/extensions/ocaml-mpfi/
17 module Interval =struct
20 let mk_interval (a,b) = { lo = a; hi = b; };;
22 let string_of_interval x = Printf.sprintf "[%f;%f]" x.lo x.hi;;
24 (* let izero = mk_interval(0.0,0.0);; *)
25 let zero = mk_interval(0.0,0.0);;
26 let one = mk_interval(1.0,1.0);;
27 let two = mk_interval(2.0,2.0);;
28 let four = mk_interval(4.0,4.0);;
30 let is_zero x =(x.lo=0.0)&&(x.hi=0.0);;
32 let pos x = if (x.lo >= 0.0) then x else
33 mk_interval(0.0, if (x.hi < 0.0) then 0.0 else x.hi );;
35 let imax (x,y) = let t=max x.hi y.hi in mk_interval(t,t);;
37 let imin (x,y) = let t = min x.lo y.lo in mk_interval(t,t);;
39 let imin3(x,y,z) = imin(x,imin(y,z));;
41 let imax3(x,y,z) = imax(x,imax(y,z));;
43 let imax4(w,x,y,z) = imax(imax(w,x),imax(y,z));;
49 let iabs x = max x.hi (~-. (x.lo));;
51 let ilt x y = (x.hi < y.lo);;
53 let igt x y = (x.lo > y.hi);;
55 let ieq x y = (x.lo = y.lo && x.hi = y.hi);;
57 (* need rounding modes -- BUG *)
60 (* start of bug section *)
62 let up() = ( (* bug *) );;
63 let down() = ( (* bug *) );;
64 let nearest() = ( (* bug *) );;
65 let upadd x y = ( x +. y);; (* bug *)
66 let upmul x y = (x *. y);;
67 let updiv x y = (x /. y);;
68 let upsub x y = (x -. y);;
69 let downadd x y = (x +. y);;
70 let downmul x y = (x *. y);;
71 let downdiv x y = (x /. y);;
72 let downsub x y = (x -. y);;
74 (* end of bug section *)
76 let interval_of_string =
77 let dbl_min =1.0e-300 in
79 let ( - ) = (down(); downsub) in
80 let lo = float_of_string s1 - dbl_min in
81 let ( + ) = (up(); upadd) in
82 let hi = float_of_string s2 + dbl_min in
85 let interval_of_single s = interval_of_string (s,s);;
87 let ineg x = mk_interval(~-. (x.hi), ~-. (x.lo));;
89 let iadd x y = mk_interval((down(); downadd x.lo y.lo), (up(); upadd x.hi y.hi));;
93 (if (y.lo >= 0.0) then (x.lo,y.lo,x.hi,y.hi)
94 else if (y.hi <= 0.0) then (x.hi,y.lo,x.lo,y.hi ) else (x.hi,y.lo,x.hi,y.hi))
95 else if (x.hi <= 0.0) then
96 (if (y.hi <= 0.0) then (x.hi,y.hi,x.lo,y.lo)
97 else if (y.lo >= 0.0) then (x.lo,y.hi,x.hi,y.lo) else (x.lo,y.hi,x.lo,y.lo))
99 (if (y.lo >=0.0) then (x.lo,y.hi,x.hi,y.hi)
100 else if (y.hi <=0.0) then (x.hi,y.lo,x.lo,y.lo)
101 else (let lo = (down(); min (downmul x.hi y.lo) (downmul x.lo y.hi)) in
102 let hi = (up(); max (upmul x.hi y.hi) (upmul x.lo y.lo)) in (lo,1.0,hi,1.0)));;
105 let (xlo,ylo,xhi,yhi) = slowcases x y in
106 mk_interval((down(); downmul xlo ylo),(up(); upmul xhi yhi));;
109 let test_slowmul x y =
110 let all = [x.lo *. y.lo; x.hi *. y.lo; x.lo *. y.hi; x.hi *. y.hi] in
111 let m = end_itlist min all in
112 let M = end_itlist max all in
113 ( mk_interval(m,M) = slowmul x y) in
114 let xs = map mk_interval [(~-. 7.0, ~-. 5.0);(~-. 3.0,9.0);(11.0,13.0)] in
115 let ys = map mk_interval [(~-. 16.0, ~-. 14.0);(~-. 10.0,12.0); (18.0,22.0)] in
116 let test i j = test_slowmul (List.nth xs i) (List.nth ys j) or
117 failwith (Printf.sprintf "%d %d" i j) in
120 let _ = test i j in ();
123 let imul x y = if (x.lo > 0.0 && y.lo > 0.0) then
124 mk_interval((down(); downmul x.lo y.lo, (up(); upmul x.hi y.hi))) else slowmul x y;;
126 let isub x y = mk_interval((down();downsub x.lo y.hi),(up(); upsub x.hi y.lo));;
129 let sqrt = Pervasives.sqrt in
130 fun x -> mk_interval(
131 (if (x.lo <= 0.0) then 0.0 else (down(); sqrt(x.lo))),
132 (if (x.hi <= 0.0) then 0.0 else (up(); sqrt(x.hi))));;
136 mk_interval((down(); atan x.lo),(up(); atan x.hi));;
140 mk_interval((down(); acos x.hi),(up(); acos x.lo));;
142 let combine x y = mk_interval(inf(imin(x,y)),sup(imax(x,y)));;
145 let random_int_seed = 81757 in
146 let _ = Random.init(random_int_seed) in
147 fun _ -> Random.float(1.0);;
149 let bounded_from_zero =
150 let dbl_epsilon = 1.0e-8 in
151 fun x-> (x.hi < ~-. dbl_epsilon or x.lo > dbl_epsilon);;
153 let idiv x y = if (bounded_from_zero y) then
154 imul x (mk_interval((down(); downdiv 1.0 y.hi),(up(); updiv 1.0 y.lo)))
155 else raise Unstable;;
157 (* overload arithmetic ops *)