src/Provers/Arith/fast_lin_arith.ML
author nipkow
Mon Jan 04 15:08:40 1999 +0100 (1999-01-04)
changeset 6056 b21813d1b701
parent 5982 aeb97860d352
child 6062 ede9fea461f5
permissions -rw-r--r--
Version 1 of linear arithmetic for nat.
nipkow@5982
     1
(*  Title:      Provers/Arith/fast_lin_arith.ML
nipkow@5982
     2
    ID:         $Id$
nipkow@5982
     3
    Author:     Tobias Nipkow
nipkow@5982
     4
    Copyright   1998  TU Munich
nipkow@5982
     5
nipkow@5982
     6
A generic linear arithmetic package. At the moment only used for nat.
nipkow@5982
     7
The two tactics provided:
nipkow@5982
     8
    lin_arith_tac:         int -> tactic
nipkow@5982
     9
cut_lin_arith_tac: thms -> int -> tactic
nipkow@5982
    10
Only take premises and conclusions into account
nipkow@5982
    11
that are already (negated) (in)equations.
nipkow@5982
    12
*)
nipkow@5982
    13
nipkow@5982
    14
(*** Data needed for setting up the linear arithmetic package ***)
nipkow@5982
    15
nipkow@5982
    16
signature LIN_ARITH_DATA =
nipkow@5982
    17
sig
nipkow@5982
    18
  val add_mono_thms:    thm list
nipkow@5982
    19
                            (* [| i rel1 j; m rel2 n |] ==> i + m rel3 j + n *)
nipkow@5982
    20
  val conjI:		thm
nipkow@5982
    21
  val ccontr:           thm (* (~ P ==> False) ==> P *)
nipkow@5982
    22
  val lessD:            thm (* m < n ==> Suc m <= n *)
nipkow@5982
    23
  val nat_neqE:         thm (* [| m ~= n; m < n ==> P; n < m ==> P |] ==> P *)
nipkow@5982
    24
  val notI:             thm (* (P ==> False) ==> ~ P *)
nipkow@5982
    25
  val not_leD:          thm (* ~(m <= n) ==> Suc n <= m *)
nipkow@5982
    26
  val not_lessD:        thm (* ~(m < n) ==> n < m *)
nipkow@5982
    27
  val sym:		thm (* x = y ==> y = x *)
nipkow@5982
    28
  val decomp: term ->
nipkow@5982
    29
             ((term * int)list * int * string * (term * int)list * int)option
nipkow@5982
    30
  val simp: thm -> thm
nipkow@6056
    31
  val is_False: thm -> bool
nipkow@6056
    32
  val is_nat: typ list * term -> bool
nipkow@6056
    33
  val mk_nat_thm: Sign.sg -> term -> thm
nipkow@5982
    34
end;
nipkow@5982
    35
(*
nipkow@5982
    36
decomp(`x Rel y') should yield (p,i,Rel,q,j)
nipkow@5982
    37
   where Rel is one of "<", "~<", "<=", "~<=" and "=" and
nipkow@5982
    38
         p/q is the decomposition of the sum terms x/y into a list
nipkow@5982
    39
         of summand * multiplicity pairs and a constant summand.
nipkow@5982
    40
nipkow@5982
    41
simp must reduce contradictory <= to False.
nipkow@5982
    42
   It should also cancel common summands to keep <= reduced;
nipkow@5982
    43
   otherwise <= can grow to massive proportions.
nipkow@6056
    44
nipkow@6056
    45
is_nat(parameter-types,t) =  t:nat
nipkow@6056
    46
mk_nat_thm(t) = "0 <= t"
nipkow@5982
    47
*)
nipkow@5982
    48
nipkow@5982
    49
functor Fast_Lin_Arith(LA_Data:LIN_ARITH_DATA) =
nipkow@5982
    50
struct
nipkow@5982
    51
nipkow@5982
    52
(*** A fast decision procedure ***)
nipkow@5982
    53
(*** Code ported from HOL Light ***)
nipkow@6056
    54
(* possible optimizations:
nipkow@6056
    55
   use (var,coeff) rep or vector rep  tp save space;
nipkow@6056
    56
   treat non-negative atoms separately rather than adding 0 <= atom
nipkow@6056
    57
*)
nipkow@5982
    58
nipkow@5982
    59
datatype lineq_type = Eq | Le | Lt;
nipkow@5982
    60
nipkow@6056
    61
datatype injust = Asm of int
nipkow@6056
    62
                | Nat of int (* index of atom *)
nipkow@5982
    63
                | Fwd of injust * thm
nipkow@5982
    64
                | Multiplied of int * injust
nipkow@5982
    65
                | Added of injust * injust;
nipkow@5982
    66
nipkow@5982
    67
datatype lineq = Lineq of int * lineq_type * int list * injust;
nipkow@5982
    68
nipkow@5982
    69
(* ------------------------------------------------------------------------- *)
nipkow@5982
    70
(* Calculate new (in)equality type after addition.                           *)
nipkow@5982
    71
(* ------------------------------------------------------------------------- *)
nipkow@5982
    72
nipkow@5982
    73
fun find_add_type(Eq,x) = x
nipkow@5982
    74
  | find_add_type(x,Eq) = x
nipkow@5982
    75
  | find_add_type(_,Lt) = Lt
nipkow@5982
    76
  | find_add_type(Lt,_) = Lt
nipkow@5982
    77
  | find_add_type(Le,Le) = Le;
nipkow@5982
    78
nipkow@5982
    79
(* ------------------------------------------------------------------------- *)
nipkow@5982
    80
(* Multiply out an (in)equation.                                             *)
nipkow@5982
    81
(* ------------------------------------------------------------------------- *)
nipkow@5982
    82
nipkow@5982
    83
fun multiply_ineq n (i as Lineq(k,ty,l,just)) =
nipkow@5982
    84
  if n = 1 then i
nipkow@5982
    85
  else if n = 0 andalso ty = Lt then sys_error "multiply_ineq"
nipkow@5982
    86
  else if n < 0 andalso (ty=Le orelse ty=Lt) then sys_error "multiply_ineq"
nipkow@5982
    87
  else Lineq(n * k,ty,map (apl(n,op * )) l,Multiplied(n,just));
nipkow@5982
    88
nipkow@5982
    89
(* ------------------------------------------------------------------------- *)
nipkow@5982
    90
(* Add together (in)equations.                                               *)
nipkow@5982
    91
(* ------------------------------------------------------------------------- *)
nipkow@5982
    92
nipkow@5982
    93
fun add_ineq (i1 as Lineq(k1,ty1,l1,just1)) (i2 as Lineq(k2,ty2,l2,just2)) =
nipkow@5982
    94
  let val l = map2 (op +) (l1,l2)
nipkow@5982
    95
  in Lineq(k1+k2,find_add_type(ty1,ty2),l,Added(just1,just2)) end;
nipkow@5982
    96
nipkow@5982
    97
(* ------------------------------------------------------------------------- *)
nipkow@5982
    98
(* Elimination of variable between a single pair of (in)equations.           *)
nipkow@5982
    99
(* If they're both inequalities, 1st coefficient must be +ve, 2nd -ve.       *)
nipkow@5982
   100
(* ------------------------------------------------------------------------- *)
nipkow@5982
   101
nipkow@5982
   102
fun gcd x y =
nipkow@5982
   103
  let fun gxd x y =
nipkow@5982
   104
    if y = 0 then x else gxd y (x mod y)
nipkow@5982
   105
  in if x < y then gxd y x else gxd x y end;
nipkow@5982
   106
nipkow@5982
   107
fun lcm x y = (x * y) div gcd x y;
nipkow@5982
   108
nipkow@5982
   109
fun el 0 (h::_) = h
nipkow@5982
   110
  | el n (_::t) = el (n - 1) t
nipkow@5982
   111
  | el _ _  = sys_error "el";
nipkow@5982
   112
nipkow@5982
   113
fun elim_var v (i1 as Lineq(k1,ty1,l1,just1)) (i2 as Lineq(k2,ty2,l2,just2)) =
nipkow@5982
   114
  let val c1 = el v l1 and c2 = el v l2
nipkow@5982
   115
      val m = lcm (abs c1) (abs c2)
nipkow@5982
   116
      val m1 = m div (abs c1) and m2 = m div (abs c2)
nipkow@5982
   117
      val (n1,n2) =
nipkow@5982
   118
        if (c1 >= 0) = (c2 >= 0)
nipkow@5982
   119
        then if ty1 = Eq then (~m1,m2)
nipkow@5982
   120
             else if ty2 = Eq then (m1,~m2)
nipkow@5982
   121
                  else sys_error "elim_var"
nipkow@5982
   122
        else (m1,m2)
nipkow@5982
   123
      val (p1,p2) = if ty1=Eq andalso ty2=Eq andalso (n1 = ~1 orelse n2 = ~1)
nipkow@5982
   124
                    then (~n1,~n2) else (n1,n2)
nipkow@5982
   125
  in add_ineq (multiply_ineq n1 i1) (multiply_ineq n2 i2) end;
nipkow@5982
   126
nipkow@5982
   127
(* ------------------------------------------------------------------------- *)
nipkow@5982
   128
(* The main refutation-finding code.                                         *)
nipkow@5982
   129
(* ------------------------------------------------------------------------- *)
nipkow@5982
   130
nipkow@5982
   131
fun is_trivial (Lineq(_,_,l,_)) = forall (fn i => i=0) l;
nipkow@5982
   132
nipkow@5982
   133
fun is_answer (ans as Lineq(k,ty,l,_)) =
nipkow@5982
   134
  case ty  of Eq => k <> 0 | Le => k > 0 | Lt => k >= 0;
nipkow@5982
   135
nipkow@5982
   136
fun calc_blowup l =
nipkow@5982
   137
  let val (p,n) = partition (apl(0,op<)) (filter (apl(0,op<>)) l)
nipkow@5982
   138
  in (length p) * (length n) end;
nipkow@5982
   139
nipkow@5982
   140
(* ------------------------------------------------------------------------- *)
nipkow@5982
   141
(* Main elimination code:                                                    *)
nipkow@5982
   142
(*                                                                           *)
nipkow@5982
   143
(* (1) Looks for immediate solutions (false assertions with no variables).   *)
nipkow@5982
   144
(*                                                                           *)
nipkow@5982
   145
(* (2) If there are any equations, picks a variable with the lowest absolute *)
nipkow@5982
   146
(* coefficient in any of them, and uses it to eliminate.                     *)
nipkow@5982
   147
(*                                                                           *)
nipkow@5982
   148
(* (3) Otherwise, chooses a variable in the inequality to minimize the       *)
nipkow@5982
   149
(* blowup (number of consequences generated) and eliminates it.              *)
nipkow@5982
   150
(* ------------------------------------------------------------------------- *)
nipkow@5982
   151
nipkow@5982
   152
fun allpairs f xs ys =
nipkow@5982
   153
  flat(map (fn x => map (fn y => f x y) ys) xs);
nipkow@5982
   154
nipkow@5982
   155
fun extract_first p =
nipkow@5982
   156
  let fun extract xs (y::ys) = if p y then (Some y,xs@ys)
nipkow@5982
   157
                               else extract (y::xs) ys
nipkow@5982
   158
        | extract xs []      = (None,xs)
nipkow@5982
   159
  in extract [] end;
nipkow@5982
   160
nipkow@6056
   161
(*
nipkow@6056
   162
fun print_ineqs ineqs =
nipkow@6056
   163
 writeln(cat_lines(""::map (fn Lineq(c,t,l,_) =>
nipkow@6056
   164
   string_of_int c ^
nipkow@6056
   165
   (case t of Eq => " =  " | Lt=> " <  " | Le => " <= ") ^
nipkow@6056
   166
   commas(map string_of_int l)) ineqs));
nipkow@6056
   167
*)
nipkow@6056
   168
nipkow@5982
   169
fun elim ineqs =
nipkow@6056
   170
  let (*val dummy = print_ineqs ineqs;*)
nipkow@6056
   171
      val (triv,nontriv) = partition is_trivial ineqs in
nipkow@5982
   172
  if not(null triv)
nipkow@5982
   173
  then case find_first is_answer triv of
nipkow@5982
   174
         None => elim nontriv | some => some
nipkow@5982
   175
  else
nipkow@5982
   176
  if null nontriv then None else
nipkow@5982
   177
  let val (eqs,noneqs) = partition (fn (Lineq(_,ty,_,_)) => ty=Eq) nontriv in
nipkow@5982
   178
  if not(null eqs) then
nipkow@5982
   179
     let val clist = foldl (fn (cs,Lineq(_,_,l,_)) => l union cs) ([],eqs)
nipkow@5982
   180
         val sclist = sort (fn (x,y) => int_ord(abs(x),abs(y)))
nipkow@5982
   181
                           (filter (fn i => i<>0) clist)
nipkow@5982
   182
         val c = hd sclist
nipkow@5982
   183
         val (Some(eq as Lineq(_,_,ceq,_)),othereqs) =
nipkow@5982
   184
               extract_first (fn Lineq(_,_,l,_) => c mem l) eqs
nipkow@5982
   185
         val v = find_index (fn k => k=c) ceq
nipkow@5982
   186
         val (ioth,roth) = partition (fn (Lineq(_,_,l,_)) => el v l = 0)
nipkow@5982
   187
                                     (othereqs @ noneqs)
nipkow@5982
   188
         val others = map (elim_var v eq) roth @ ioth
nipkow@5982
   189
     in elim others end
nipkow@5982
   190
  else
nipkow@5982
   191
  let val lists = map (fn (Lineq(_,_,l,_)) => l) noneqs
nipkow@5982
   192
      val numlist = 0 upto (length(hd lists) - 1)
nipkow@5982
   193
      val coeffs = map (fn i => map (el i) lists) numlist
nipkow@5982
   194
      val blows = map calc_blowup coeffs
nipkow@5982
   195
      val iblows = blows ~~ numlist
nipkow@5982
   196
      val nziblows = filter (fn (i,_) => i<>0) iblows
nipkow@5982
   197
  in if null nziblows then None else
nipkow@5982
   198
     let val (c,v) = hd(sort (fn (x,y) => int_ord(fst(x),fst(y))) nziblows)
nipkow@5982
   199
         val (no,yes) = partition (fn (Lineq(_,_,l,_)) => el v l = 0) ineqs
nipkow@5982
   200
         val (pos,neg) = partition(fn (Lineq(_,_,l,_)) => el v l > 0) yes
nipkow@5982
   201
     in elim (no @ allpairs (elim_var v) pos neg) end
nipkow@5982
   202
  end
nipkow@5982
   203
  end
nipkow@5982
   204
  end;
nipkow@5982
   205
nipkow@5982
   206
(* ------------------------------------------------------------------------- *)
nipkow@5982
   207
(* Translate back a proof.                                                   *)
nipkow@5982
   208
(* ------------------------------------------------------------------------- *)
nipkow@5982
   209
nipkow@6056
   210
(* FIXME OPTIMIZE!!!!
nipkow@6056
   211
   Addition/Multiplication need i*t representation rather than t+t+...
nipkow@6056
   212
nipkow@6056
   213
Simplification may detect a contradiction 'prematurely' due to type
nipkow@6056
   214
information: n+1 <= 0 is simplified to False and does not need to be crossed
nipkow@6056
   215
with 0 <= n.
nipkow@6056
   216
*)
nipkow@6056
   217
local
nipkow@6056
   218
 exception FalseE of thm
nipkow@6056
   219
in
nipkow@6056
   220
fun mkproof sg asms just =
nipkow@6056
   221
  let val atoms = foldl (fn (ats,(lhs,_,_,rhs,_)) =>
nipkow@6056
   222
                            map fst lhs  union  (map fst rhs  union  ats))
nipkow@6056
   223
                        ([], mapfilter (LA_Data.decomp o concl_of) asms)
nipkow@6056
   224
nipkow@6056
   225
      fun addthms thm1 thm2 =
nipkow@5982
   226
        let val conj = thm1 RS (thm2 RS LA_Data.conjI)
nipkow@5982
   227
        in the(get_first (fn th => Some(conj RS th) handle _ => None)
nipkow@5982
   228
                         LA_Data.add_mono_thms)
nipkow@5982
   229
        end;
nipkow@5982
   230
nipkow@5982
   231
      fun multn(n,thm) =
nipkow@5982
   232
        let fun mul(i,th) = if i=1 then th else mul(i-1, addthms thm th)
nipkow@5982
   233
        in if n < 0 then mul(~n,thm) RS LA_Data.sym else mul(n,thm) end;
nipkow@5982
   234
nipkow@6056
   235
      fun simp thm =
nipkow@6056
   236
        let val thm' = LA_Data.simp thm
nipkow@6056
   237
        in if LA_Data.is_False thm' then raise FalseE thm' else thm' end
nipkow@6056
   238
nipkow@6056
   239
      fun mk(Asm i) = nth_elem(i,asms)
nipkow@6056
   240
        | mk(Nat(i)) = LA_Data.mk_nat_thm sg (nth_elem(i,atoms))
nipkow@5982
   241
        | mk(Fwd(j,thm)) = mk j RS thm
nipkow@6056
   242
        | mk(Added(j1,j2)) = simp(addthms (mk j1) (mk j2))
nipkow@5982
   243
        | mk(Multiplied(n,j)) = multn(n,mk j)
nipkow@5982
   244
nipkow@6056
   245
  in LA_Data.simp(mk just) handle FalseE thm => thm end
nipkow@6056
   246
end;
nipkow@5982
   247
nipkow@5982
   248
fun coeff poly atom = case assoc(poly,atom) of None => 0 | Some i => i;
nipkow@5982
   249
nipkow@5982
   250
fun mklineq atoms =
nipkow@5982
   251
  let val n = length atoms in
nipkow@5982
   252
    fn ((lhs,i,rel,rhs,j),k) =>
nipkow@5982
   253
    let val lhsa = map (coeff lhs) atoms
nipkow@5982
   254
        and rhsa = map (coeff rhs) atoms
nipkow@5982
   255
        val diff = map2 (op -) (rhsa,lhsa)
nipkow@5982
   256
        val c = i-j
nipkow@6056
   257
        val just = Asm k
nipkow@5982
   258
    in case rel of
nipkow@5982
   259
        "<="   => Some(Lineq(c,Le,diff,just))
nipkow@5982
   260
       | "~<=" => Some(Lineq(1-c,Le,map (op ~) diff,Fwd(just,LA_Data.not_leD)))
nipkow@5982
   261
       | "<"   => Some(Lineq(c+1,Le,diff,Fwd(just,LA_Data.lessD)))
nipkow@5982
   262
       | "~<"  => Some(Lineq(~c,Le,map (op~) diff,Fwd(just,LA_Data.not_lessD)))
nipkow@5982
   263
       | "="   => Some(Lineq(c,Eq,diff,just))
nipkow@5982
   264
       | "~="  => None
nipkow@5982
   265
       | _     => sys_error("mklineq" ^ rel)   
nipkow@5982
   266
    end
nipkow@5982
   267
  end;
nipkow@5982
   268
nipkow@6056
   269
fun mknat pTs ixs (atom,i) =
nipkow@6056
   270
  if LA_Data.is_nat(pTs,atom)
nipkow@6056
   271
  then let val l = map (fn j => if j=i then 1 else 0) ixs
nipkow@6056
   272
       in Some(Lineq(0,Le,l,Nat(i))) end
nipkow@6056
   273
  else None
nipkow@6056
   274
nipkow@6056
   275
fun abstract pTs items =
nipkow@5982
   276
  let val atoms = foldl (fn (ats,((lhs,_,_,rhs,_),_)) =>
nipkow@5982
   277
                            (map fst lhs) union ((map fst rhs) union ats))
nipkow@5982
   278
                        ([],items)
nipkow@6056
   279
      val ixs = 0 upto (length(atoms)-1)
nipkow@6056
   280
      val iatoms = atoms ~~ ixs
nipkow@6056
   281
  in mapfilter (mklineq atoms) items @ mapfilter (mknat pTs ixs) iatoms end;
nipkow@5982
   282
nipkow@5982
   283
(* Ordinary refutation *)
nipkow@6056
   284
fun refute1_tac pTs items =
nipkow@6056
   285
  let val lineqs = abstract pTs items
nipkow@5982
   286
  in case elim lineqs of
nipkow@5982
   287
       None => K no_tac
nipkow@5982
   288
     | Some(Lineq(_,_,_,j)) =>
nipkow@6056
   289
         fn i => fn state => 
nipkow@6056
   290
         let val sg = #sign(rep_thm state)
nipkow@6056
   291
         in resolve_tac [LA_Data.notI,LA_Data.ccontr] i THEN
nipkow@6056
   292
            METAHYPS (fn asms => rtac (mkproof sg asms j) 1) i
nipkow@6056
   293
         end state
nipkow@5982
   294
  end;
nipkow@5982
   295
nipkow@5982
   296
(* Double refutation caused by equality in conclusion *)
nipkow@6056
   297
fun refute2_tac pTs items (rhs,i,_,lhs,j) nHs =
nipkow@6056
   298
  (case elim (abstract pTs (items@[((rhs,i,"<",lhs,j),nHs)])) of
nipkow@5982
   299
    None => K no_tac
nipkow@5982
   300
  | Some(Lineq(_,_,_,j1)) =>
nipkow@6056
   301
      (case elim (abstract pTs (items@[((lhs,j,"<",rhs,i),nHs)])) of
nipkow@5982
   302
        None => K no_tac
nipkow@5982
   303
      | Some(Lineq(_,_,_,j2)) =>
nipkow@6056
   304
          fn i => fn state => 
nipkow@6056
   305
          let val sg = #sign(rep_thm state)
nipkow@6056
   306
          in rtac LA_Data.ccontr i THEN etac LA_Data.nat_neqE i THEN
nipkow@6056
   307
             METAHYPS (fn asms => rtac (mkproof sg asms j1) 1) i THEN
nipkow@6056
   308
             METAHYPS (fn asms => rtac (mkproof sg asms j2) 1) i
nipkow@6056
   309
          end state));
nipkow@5982
   310
nipkow@5982
   311
(*
nipkow@5982
   312
Fast but very incomplete decider. Only premises and conclusions
nipkow@5982
   313
that are already (negated) (in)equations are taken into account.
nipkow@5982
   314
*)
nipkow@5982
   315
val lin_arith_tac = SUBGOAL (fn (A,n) =>
nipkow@6056
   316
  let val pTs = rev(map snd (Logic.strip_params A))
nipkow@6056
   317
      val Hs = Logic.strip_assums_hyp A
nipkow@5982
   318
      val nHs = length Hs
nipkow@5982
   319
      val His = Hs ~~ (0 upto (nHs-1))
nipkow@5982
   320
      val Hitems = mapfilter (fn (h,i) => case LA_Data.decomp h of
nipkow@5982
   321
                                 None => None | Some(it) => Some(it,i)) His
nipkow@5982
   322
  in case LA_Data.decomp(Logic.strip_assums_concl A) of
nipkow@6056
   323
       None => if null Hitems then no_tac else refute1_tac pTs Hitems n
nipkow@5982
   324
     | Some(citem as (r,i,rel,l,j)) =>
nipkow@5982
   325
         if rel = "="
nipkow@6056
   326
         then refute2_tac pTs Hitems citem nHs n
nipkow@5982
   327
         else let val neg::rel0 = explode rel
nipkow@5982
   328
                  val nrel = if neg = "~" then implode rel0 else "~"^rel
nipkow@6056
   329
              in refute1_tac pTs (Hitems@[((r,i,nrel,l,j),nHs)]) n end
nipkow@5982
   330
  end);
nipkow@5982
   331
nipkow@5982
   332
fun cut_lin_arith_tac thms i = cut_facts_tac thms i THEN lin_arith_tac i;
nipkow@5982
   333
nipkow@5982
   334
end;