src/Provers/Arith/fast_lin_arith.ML
changeset 14372 51ddf8963c95
parent 14360 e654599b114e
child 14386 ad1ffcc90162
     1.1 --- a/src/Provers/Arith/fast_lin_arith.ML	Mon Feb 02 12:23:46 2004 +0100
     1.2 +++ b/src/Provers/Arith/fast_lin_arith.ML	Tue Feb 03 10:19:21 2004 +0100
     1.3 @@ -153,36 +153,40 @@
     1.4  (* Finding a (counter) example from the trace of a failed elimination        *)
     1.5  (* ------------------------------------------------------------------------- *)
     1.6  (* Examples are represented as rational numbers,                             *)
     1.7 -(* although at the moment all examples are rounded to integers -             *)
     1.8 -(* thus it does not yet work for type real.                                  *)
     1.9  (* Dont blame John Harrison for this code - it is entirely mine. TN          *)
    1.10  
    1.11  exception NoEx;
    1.12 -exception NotYetImpl;
    1.13  
    1.14 -fun elim_eqns(ineqs,Lineq(i,Le,cs,_)) = (i,cs)::ineqs
    1.15 -  | elim_eqns(ineqs,Lineq(i,Eq,cs,_)) = (i,cs)::(~i,map ~ cs)::ineqs
    1.16 -  | elim_eqns(ineqs,Lineq(i,Lt,cs,_)) = raise NotYetImpl;
    1.17 +(* Coding: (i,true,cs) means i <= cs and (i,false,cs) means i < cs.
    1.18 +   In general, true means the bound is included, false means it is excluded.
    1.19 +   Need to know if it is a lower or upper bound for unambiguous interpretation!
    1.20 +*)
    1.21 +
    1.22 +fun elim_eqns(ineqs,Lineq(i,Le,cs,_)) = (i,true,cs)::ineqs
    1.23 +  | elim_eqns(ineqs,Lineq(i,Eq,cs,_)) = (i,true,cs)::(~i,true,map ~ cs)::ineqs
    1.24 +  | elim_eqns(ineqs,Lineq(i,Lt,cs,_)) = (i,false,cs)::ineqs;
    1.25  
    1.26  val rat0 = rat_of_int 0;
    1.27  
    1.28  (* PRE: ex[v] must be 0! *)
    1.29 -fun eval (ex:rat list) v (a:int,cs:int list) =
    1.30 +fun eval (ex:rat list) v (a:int,le,cs:int list) =
    1.31    let val rs = map rat_of_int cs
    1.32        val rsum = foldl ratadd (rat0,map ratmul (rs ~~ ex))
    1.33 -  in ratmul(ratadd(rat_of_int a,ratneg rsum), ratinv(el v rs)) end;
    1.34 +  in (ratmul(ratadd(rat_of_int a,ratneg rsum), ratinv(el v rs)), le) end;
    1.35 +(* If el v rs < 0, le should be negated.
    1.36 +   Instead this swap is taken into account in ratrelmin2.
    1.37 +*)
    1.38  
    1.39 -(*
    1.40 -fun ratge0(Rat(a,p,q)) = (p = 0 orelse a)
    1.41 -*)
    1.42  fun ratge0 r = fst(rep_rat r) >= 0;
    1.43 -fun ratle(r,s) = ratge0(ratadd(s,ratneg r))
    1.44 +fun ratle(r,s) = ratge0(ratadd(s,ratneg r));
    1.45  
    1.46 -fun ratmin2(r,s) = if ratle(r,s) then r else s;
    1.47 -fun ratmax2(r,s) = if ratle(r,s) then s else r;
    1.48 +fun ratrelmin2(x as (r,ler),y as (s,les)) =
    1.49 +  if r=s then (r, (not ler) andalso (not les)) else if ratle(r,s) then x else y;
    1.50 +fun ratrelmax2(x as (r,ler),y as (s,les)) =
    1.51 +  if r=s then (r,ler andalso les) else if ratle(r,s) then y else x;
    1.52  
    1.53 -val ratmin = foldr1 ratmin2;
    1.54 -val ratmax = foldr1 ratmax2;
    1.55 +val ratrelmin = foldr1 ratrelmin2;
    1.56 +val ratrelmax = foldr1 ratrelmax2;
    1.57  
    1.58  fun ratroundup r = let val (p,q) = rep_rat r
    1.59                     in if q=1 then r else rat_of_int((p div q) + 1) end
    1.60 @@ -190,9 +194,23 @@
    1.61  fun ratrounddown r = let val (p,q) = rep_rat r
    1.62                       in if q=1 then r else rat_of_int((p div q) - 1) end
    1.63  
    1.64 -fun choose2 d (lb,ub) =
    1.65 -  if ratle(lb,rat0) andalso ratle(rat0,ub) then rat0 else
    1.66 -  if not d then (if ratge0 lb then lb else ub) else
    1.67 +fun ratexact up (r,exact) =
    1.68 +  if exact then r else
    1.69 +  let val (p,q) = rep_rat r
    1.70 +      val nth = ratinv(rat_of_int q)
    1.71 +  in ratadd(r,if up then nth else ratneg nth) end;
    1.72 +
    1.73 +fun ratmiddle(r,s) = ratmul(ratadd(r,s),ratinv(rat_of_int 2));
    1.74 +
    1.75 +fun choose2 d ((lb,exactl),(ub,exactu)) =
    1.76 +  if ratle(lb,rat0) andalso (lb <> rat0 orelse exactl) andalso
    1.77 +     ratle(rat0,ub) andalso (ub <> rat0 orelse exactu)
    1.78 +  then rat0 else
    1.79 +  if not d
    1.80 +  then (if ratge0 lb
    1.81 +        then if exactl then lb else ratmiddle(lb,ub)
    1.82 +        else if exactu then ub else ratmiddle(lb,ub))
    1.83 +  else (* discrete domain, both bounds must be exact *)
    1.84    if ratge0 lb then let val lb' = ratroundup lb
    1.85                      in if ratle(lb',ub) then lb' else raise NoEx end
    1.86                 else let val ub' = ratrounddown ub
    1.87 @@ -201,40 +219,45 @@
    1.88  fun findex1 discr (ex,(v,lineqs)) =
    1.89    let val nz = filter (fn (Lineq(_,_,cs,_)) => el v cs <> 0) lineqs;
    1.90        val ineqs = foldl elim_eqns ([],nz)
    1.91 -      val (ge,le) = partition (fn (_,cs) => el v cs > 0) ineqs
    1.92 -      val lb = ratmax(map (eval ex v) ge)
    1.93 -      val ub = ratmin(map (eval ex v) le)
    1.94 +      val (ge,le) = partition (fn (_,_,cs) => el v cs > 0) ineqs
    1.95 +      val lb = ratrelmax(map (eval ex v) ge)
    1.96 +      val ub = ratrelmin(map (eval ex v) le)
    1.97    in nth_update (choose2 (nth_elem(v,discr)) (lb,ub)) (v,ex) end;
    1.98  
    1.99  fun findex discr = foldl (findex1 discr);
   1.100  
   1.101  fun elim1 v x =
   1.102 -  map (fn (a,bs) => (ratadd(a,ratneg(ratmul(el v bs,x))),
   1.103 -                     nth_update rat0 (v,bs)));
   1.104 +  map (fn (a,le,bs) => (ratadd(a,ratneg(ratmul(el v bs,x))), le,
   1.105 +                        nth_update rat0 (v,bs)));
   1.106  
   1.107 -fun single_var v cs = (filter_out (equal rat0) cs = [el v cs]);
   1.108 +fun single_var v (_,_,cs) = (filter_out (equal rat0) cs = [el v cs]);
   1.109  
   1.110  (* The base case:
   1.111     all variables occur only with positive or only with negative coefficients *)
   1.112  fun pick_vars discr (ineqs,ex) =
   1.113 -  let val nz = filter_out (forall (equal rat0) o snd) ineqs
   1.114 -  in if null nz then ex
   1.115 -     else let val v = find_index (not o equal rat0) (snd(hd nz))
   1.116 -              val d = nth_elem(v,discr)
   1.117 -              val sv = filter (single_var v o snd) nz
   1.118 -              val minmax = if ratge0(el v (snd(hd nz)))
   1.119 -                           then if d then ratroundup o ratmax else ratmax
   1.120 -                           else if d then ratrounddown o ratmin else ratmin
   1.121 -              val bnds = map (fn (a,bs) => ratmul(a,ratinv(el v bs))) sv
   1.122 -              val x = minmax(rat0::bnds)
   1.123 -              val ineqs' = elim1 v x nz
   1.124 -              val ex' = nth_update x (v,ex)
   1.125 -          in pick_vars discr (ineqs',ex') end
   1.126 +  let val nz = filter_out (fn (_,_,cs) => forall (equal rat0) cs) ineqs
   1.127 +  in case nz of [] => ex
   1.128 +     | (_,_,cs) :: _ =>
   1.129 +       let val v = find_index (not o equal rat0) cs
   1.130 +           val d = nth_elem(v,discr)
   1.131 +           val pos = ratge0(el v cs)
   1.132 +           val sv = filter (single_var v) nz
   1.133 +           val minmax =
   1.134 +             if pos then if d then ratroundup o fst o ratrelmax
   1.135 +                         else ratexact true o ratrelmax
   1.136 +                    else if d then ratrounddown o fst o ratrelmin
   1.137 +                         else ratexact false o ratrelmin
   1.138 +           val bnds = map (fn (a,le,bs) => (ratmul(a,ratinv(el v bs)),le)) sv
   1.139 +           val x = minmax((rat0,if pos then true else false)::bnds)
   1.140 +           val ineqs' = elim1 v x nz
   1.141 +           val ex' = nth_update x (v,ex)
   1.142 +       in pick_vars discr (ineqs',ex') end
   1.143    end;
   1.144  
   1.145  fun findex0 discr n lineqs =
   1.146    let val ineqs = foldl elim_eqns ([],lineqs)
   1.147 -      val rineqs = map (fn (a,cs) => (rat_of_int a, map rat_of_int cs)) ineqs
   1.148 +      val rineqs = map (fn (a,le,cs) => (rat_of_int a, le, map rat_of_int cs))
   1.149 +                       ineqs
   1.150    in pick_vars discr (rineqs,replicate n rat0) end;
   1.151  
   1.152  (* ------------------------------------------------------------------------- *)