formal_ineqs/verifier/interval_m/interval.ml

   1 (* =========================================================== *)
   2 (* OCaml interval arithmetic                                   *)
   3 (* Author: Thomas C. Hales                                     *)
   4 (* Date: 2011-08-21                                            *)
   5 (* =========================================================== *)
   6
   7 (* port of interval.cc,
   8
   9 This file gives a simple implementation of interval arithmetic,
  10 together with the basic arithmetic operations on intervals.
  11
  12 It has been incompletely implemented.
  13
  14 For now, I am not implementing directed roundings.
  15 However, McLaughlin implemented directed rounding several years ago:
  16 See http://perso.ens-lyon.fr/nathalie.revol/mpfi.html
  17  ~/Library/McLaughlinOCAML/ocaml/src/extensions/ocaml-mpfi/
  18
  19  *)
  20
  21 needs "verifier/interval_m/types.ml";;
  22
  23 module Interval = struct
  24
  25 open Interval_types;;
  26
  27 let mk_interval (a,b) = { lo = a; hi = b; };;
  28
  29 let string_of_interval x = Printf.sprintf "[%f;%f]" x.lo x.hi;;
  30
  31 (* let izero = mk_interval(0.0,0.0);; *)
  32 let zero = mk_interval(0.0,0.0);;
  33 let one = mk_interval(1.0,1.0);;
  34 let two = mk_interval(2.0,2.0);;
  35 let four = mk_interval(4.0,4.0);;
  36
  37 let is_zero x =(x.lo=0.0)&&(x.hi=0.0);;
  38
  39 let pos x = if (x.lo >= 0.0) then x else
  40    mk_interval(0.0,   if (x.hi < 0.0) then 0.0 else x.hi );;
  41
  42 let imax (x,y) = let t=max x.hi y.hi in mk_interval(t,t);;
  43
  44 let imin (x,y) = let t = min x.lo y.lo in mk_interval(t,t);;
  45
  46 let imin3(x,y,z) = imin(x,imin(y,z));;
  47
  48 let imax3(x,y,z) = imax(x,imax(y,z));;
  49
  50 let imax4(w,x,y,z) = imax(imax(w,x),imax(y,z));;
  51
  52 let sup x = x.hi;;
  53
  54 let inf x = x.lo;;
  55
  56 let iabs x = max x.hi (~-. (x.lo));;
  57
  58 let ilt x y = (x.hi < y.lo);;
  59
  60 let igt x y = (x.lo > y.hi);;
  61
  62 let ieq x y = (x.lo = y.lo && x.hi = y.hi);;
  63
  64 (* need rounding modes -- BUG *)
  65
  66
  67 (* start of bug section *)
  68
  69 let up() = (  (* bug *) );;
  70 let down() = ( (* bug *) );;
  71 let nearest() = ( (* bug *) );;
  72 let upadd x y = ( x +. y);;  (* bug *)
  73 let upmul x y = (x *. y);;
  74 let updiv x y = (x /. y);;
  75 let upsub x y = (x -. y);;
  76 let downadd x y = (x +. y);;
  77 let downmul x y = (x *. y);;
  78 let downdiv x y = (x /. y);;
  79 let downsub x y = (x -. y);;
  80
  81 (* end of bug section *)
  82
  83 let interval_of_string =
  84   let dbl_min =1.0e-300 in
  85     fun (s1,s2) ->
  86       let ( - ) = (down(); downsub) in
  87       let lo = float_of_string s1 - dbl_min in
  88       let ( + ) = (up(); upadd) in
  89       let hi = float_of_string s2 + dbl_min in
  90         mk_interval(lo,hi);;
  91
  92 let interval_of_single s = interval_of_string (s,s);;
  93
  94 let ineg x = mk_interval(~-. (x.hi), ~-. (x.lo));;
  95
  96 let iadd x y = mk_interval((down(); downadd x.lo  y.lo), (up(); upadd x.hi y.hi));;
  97
  98 let slowcases x y =
  99   if (x.lo >= 0.0) then
 100     (if (y.lo >= 0.0) then (x.lo,y.lo,x.hi,y.hi)
 101     else if (y.hi <= 0.0) then (x.hi,y.lo,x.lo,y.hi ) else (x.hi,y.lo,x.hi,y.hi))
 102   else if (x.hi <= 0.0) then
 103     (if (y.hi <= 0.0) then (x.hi,y.hi,x.lo,y.lo)
 104      else if (y.lo >= 0.0) then (x.lo,y.hi,x.hi,y.lo) else (x.lo,y.hi,x.lo,y.lo))
 105   else
 106     (if (y.lo >=0.0) then (x.lo,y.hi,x.hi,y.hi)
 107      else if (y.hi <=0.0) then (x.hi,y.lo,x.lo,y.lo)
 108      else (let lo = (down(); min (downmul x.hi  y.lo) (downmul x.lo  y.hi)) in
 109            let hi = (up(); max (upmul x.hi  y.hi) (upmul x.lo  y.lo)) in (lo,1.0,hi,1.0)));;
 110
 111 let slowmul x y =
 112   let (xlo,ylo,xhi,yhi) = slowcases x y in
 113   mk_interval((down(); downmul xlo  ylo),(up(); upmul xhi  yhi));;
 114
 115 let  _ =
 116   let test_slowmul x y =
 117     let all = [x.lo *. y.lo; x.hi *. y.lo; x.lo *. y.hi; x.hi *. y.hi] in
 118     let m = end_itlist min all in
 119     let M = end_itlist max all in
 120       ( mk_interval(m,M) = slowmul x  y) in
 121   let xs = map mk_interval [(~-. 7.0, ~-. 5.0);(~-. 3.0,9.0);(11.0,13.0)] in
 122   let ys = map mk_interval [(~-. 16.0, ~-. 14.0);(~-. 10.0,12.0); (18.0,22.0)] in
 123   let test i j = test_slowmul (List.nth xs i) (List.nth ys j) or
 124     failwith (Printf.sprintf "%d %d" i j) in
 125     for  i=0 to 2 do
 126       for j= 0 to 2 do
 127         let _ = test i j in ();
 128       done; done;;
 129
 130 let imul x y = if (x.lo > 0.0 && y.lo > 0.0) then
 131   mk_interval((down(); downmul x.lo  y.lo, (up(); upmul x.hi  y.hi))) else slowmul x y;;
 132
 133 let isub x y = mk_interval((down();downsub x.lo  y.hi),(up(); upsub x.hi  y.lo));;
 134
 135 let isqrt   =
 136   let sqrt = Pervasives.sqrt in
 137     fun x -> mk_interval(
 138   (if (x.lo <= 0.0) then 0.0 else (down(); sqrt(x.lo))),
 139   (if (x.hi <= 0.0) then 0.0 else (up(); sqrt(x.hi))));;
 140
 141 let iatan x =
 142   let _ = nearest() in
 143     mk_interval((down(); atan x.lo),(up(); atan x.hi));;
 144
 145 let iacos x =
 146   let _ = nearest() in
 147     mk_interval((down(); acos x.hi),(up(); acos x.lo));;
 148
 149 let combine x y = mk_interval(inf(imin(x,y)),sup(imax(x,y)));;
 150
 151 let rand01 =
 152   let random_int_seed = 81757 in
 153   let _ = Random.init(random_int_seed) in
 154     fun _ -> Random.float(1.0);;
 155
 156 let bounded_from_zero  =
 157   let  dbl_epsilon = 1.0e-8 in
 158     fun x-> (x.hi < ~-. dbl_epsilon or  x.lo > dbl_epsilon);;
 159
 160 let idiv x y = if (bounded_from_zero y) then
 161                      imul x  (mk_interval((down(); downdiv 1.0 y.hi),(up(); updiv 1.0 y.lo)))
 162                    else raise Unstable;;
 163
 164 (* overload arithmetic ops *)
 165
 166 (*
 167 let (+) = iadd;;
 168 let (-) = isub;;
 169 let (/) = idiv;;
 170 let (~-) = ineg;;
 171 *)
 172
 173 end;;