src/Provers/Arith/fast_lin_arith.ML
author nipkow
Tue Sep 21 14:16:08 1999 +0200 (1999-09-21)
changeset 7551 8e934d1a9ac6
parent 6128 2acc5d36610c
child 7552 0d6d1f50b86d
permissions -rw-r--r--
Now distinguishes discrete from non-distrete types.
     1 (*  Title:      Provers/Arith/fast_lin_arith.ML
     2     ID:         $Id$
     3     Author:     Tobias Nipkow
     4     Copyright   1998  TU Munich
     5 
     6 A generic linear arithmetic package.
     7 It provides two tactics
     8 
     9     lin_arith_tac:         int -> tactic
    10 cut_lin_arith_tac: thms -> int -> tactic
    11 
    12 and a simplification procedure
    13 
    14     lin_arith_prover: Sign.sg -> thm list -> term -> thm option
    15 
    16 Only take premises and conclusions into account that are already (negated)
    17 (in)equations. lin_arith_prover tries to prove or disprove the term.
    18 *)
    19 
    20 (*** Data needed for setting up the linear arithmetic package ***)
    21 
    22 signature LIN_ARITH_LOGIC =
    23 sig
    24   val conjI:		thm
    25   val ccontr:           thm (* (~ P ==> False) ==> P *)
    26   val neqE:             thm (* [| m ~= n; m < n ==> P; n < m ==> P |] ==> P *)
    27   val notI:             thm (* (P ==> False) ==> ~ P *)
    28   val not_lessD:        thm (* ~(m < n) ==> n <= m *)
    29   val not_leD:          thm (* ~(m <= n) ==> n < m *)
    30   val sym:		thm (* x = y ==> y = x *)
    31   val mk_Eq: thm -> thm
    32   val mk_Trueprop: term -> term
    33   val neg_prop: term -> term
    34   val is_False: thm -> bool
    35   val is_nat: typ list * term -> bool
    36   val mk_nat_thm: Sign.sg -> term -> thm
    37 end;
    38 (*
    39 mk_Eq(~in) = `in == False'
    40 mk_Eq(in) = `in == True'
    41 where `in' is an (in)equality.
    42 
    43 neg_prop(t) = neg if t is wrapped up in Trueprop and
    44   nt is the (logically) negated version of t, where the negation
    45   of a negative term is the term itself (no double negation!);
    46 
    47 is_nat(parameter-types,t) =  t:nat
    48 mk_nat_thm(t) = "0 <= t"
    49 *)
    50 
    51 signature LIN_ARITH_DATA =
    52 sig
    53   val add_mono_thms:    thm list ref
    54                             (* [| i rel1 j; m rel2 n |] ==> i + m rel3 j + n *)
    55   val lessD:            thm list ref (* m < n ==> m+1 <= n *)
    56   val decomp:
    57     term ->
    58       ((term * int)list * int * string * (term * int)list * int * bool)option
    59   val simp: (thm -> thm) ref
    60 end;
    61 (*
    62 decomp(`x Rel y') should yield (p,i,Rel,q,j,d)
    63    where Rel is one of "<", "~<", "<=", "~<=" and "=" and
    64          p/q is the decomposition of the sum terms x/y into a list
    65          of summand * multiplicity pairs and a constant summand and
    66          d indicates if the domain is discrete.
    67 
    68 simp must reduce contradictory <= to False.
    69    It should also cancel common summands to keep <= reduced;
    70    otherwise <= can grow to massive proportions.
    71 *)
    72 
    73 signature FAST_LIN_ARITH =
    74 sig
    75   val lin_arith_prover: Sign.sg -> thm list -> term -> thm option
    76   val     lin_arith_tac:             int -> tactic
    77   val cut_lin_arith_tac: thm list -> int -> tactic
    78 end;
    79 
    80 functor Fast_Lin_Arith(structure LA_Logic:LIN_ARITH_LOGIC 
    81                        and       LA_Data:LIN_ARITH_DATA) : FAST_LIN_ARITH =
    82 struct
    83 
    84 (*** A fast decision procedure ***)
    85 (*** Code ported from HOL Light ***)
    86 (* possible optimizations:
    87    use (var,coeff) rep or vector rep  tp save space;
    88    treat non-negative atoms separately rather than adding 0 <= atom
    89 *)
    90 
    91 datatype lineq_type = Eq | Le | Lt;
    92 
    93 datatype injust = Asm of int
    94                 | Nat of int (* index of atom *)
    95                 | LessD of injust
    96                 | NotLessD of injust
    97                 | NotLeD of injust
    98                 | NotLeDD of injust
    99                 | Multiplied of int * injust
   100                 | Added of injust * injust;
   101 
   102 datatype lineq = Lineq of int * lineq_type * int list * injust;
   103 
   104 (* ------------------------------------------------------------------------- *)
   105 (* Calculate new (in)equality type after addition.                           *)
   106 (* ------------------------------------------------------------------------- *)
   107 
   108 fun find_add_type(Eq,x) = x
   109   | find_add_type(x,Eq) = x
   110   | find_add_type(_,Lt) = Lt
   111   | find_add_type(Lt,_) = Lt
   112   | find_add_type(Le,Le) = Le;
   113 
   114 (* ------------------------------------------------------------------------- *)
   115 (* Multiply out an (in)equation.                                             *)
   116 (* ------------------------------------------------------------------------- *)
   117 
   118 fun multiply_ineq n (i as Lineq(k,ty,l,just)) =
   119   if n = 1 then i
   120   else if n = 0 andalso ty = Lt then sys_error "multiply_ineq"
   121   else if n < 0 andalso (ty=Le orelse ty=Lt) then sys_error "multiply_ineq"
   122   else Lineq(n * k,ty,map (apl(n,op * )) l,Multiplied(n,just));
   123 
   124 (* ------------------------------------------------------------------------- *)
   125 (* Add together (in)equations.                                               *)
   126 (* ------------------------------------------------------------------------- *)
   127 
   128 fun add_ineq (i1 as Lineq(k1,ty1,l1,just1)) (i2 as Lineq(k2,ty2,l2,just2)) =
   129   let val l = map2 (op +) (l1,l2)
   130   in Lineq(k1+k2,find_add_type(ty1,ty2),l,Added(just1,just2)) end;
   131 
   132 (* ------------------------------------------------------------------------- *)
   133 (* Elimination of variable between a single pair of (in)equations.           *)
   134 (* If they're both inequalities, 1st coefficient must be +ve, 2nd -ve.       *)
   135 (* ------------------------------------------------------------------------- *)
   136 
   137 fun gcd x y =
   138   let fun gxd x y =
   139     if y = 0 then x else gxd y (x mod y)
   140   in if x < y then gxd y x else gxd x y end;
   141 
   142 fun lcm x y = (x * y) div gcd x y;
   143 
   144 fun el 0 (h::_) = h
   145   | el n (_::t) = el (n - 1) t
   146   | el _ _  = sys_error "el";
   147 
   148 fun elim_var v (i1 as Lineq(k1,ty1,l1,just1)) (i2 as Lineq(k2,ty2,l2,just2)) =
   149   let val c1 = el v l1 and c2 = el v l2
   150       val m = lcm (abs c1) (abs c2)
   151       val m1 = m div (abs c1) and m2 = m div (abs c2)
   152       val (n1,n2) =
   153         if (c1 >= 0) = (c2 >= 0)
   154         then if ty1 = Eq then (~m1,m2)
   155              else if ty2 = Eq then (m1,~m2)
   156                   else sys_error "elim_var"
   157         else (m1,m2)
   158       val (p1,p2) = if ty1=Eq andalso ty2=Eq andalso (n1 = ~1 orelse n2 = ~1)
   159                     then (~n1,~n2) else (n1,n2)
   160   in add_ineq (multiply_ineq n1 i1) (multiply_ineq n2 i2) end;
   161 
   162 (* ------------------------------------------------------------------------- *)
   163 (* The main refutation-finding code.                                         *)
   164 (* ------------------------------------------------------------------------- *)
   165 
   166 fun is_trivial (Lineq(_,_,l,_)) = forall (fn i => i=0) l;
   167 
   168 fun is_answer (ans as Lineq(k,ty,l,_)) =
   169   case ty  of Eq => k <> 0 | Le => k > 0 | Lt => k >= 0;
   170 
   171 fun calc_blowup l =
   172   let val (p,n) = partition (apl(0,op<)) (filter (apl(0,op<>)) l)
   173   in (length p) * (length n) end;
   174 
   175 (* ------------------------------------------------------------------------- *)
   176 (* Main elimination code:                                                    *)
   177 (*                                                                           *)
   178 (* (1) Looks for immediate solutions (false assertions with no variables).   *)
   179 (*                                                                           *)
   180 (* (2) If there are any equations, picks a variable with the lowest absolute *)
   181 (* coefficient in any of them, and uses it to eliminate.                     *)
   182 (*                                                                           *)
   183 (* (3) Otherwise, chooses a variable in the inequality to minimize the       *)
   184 (* blowup (number of consequences generated) and eliminates it.              *)
   185 (* ------------------------------------------------------------------------- *)
   186 
   187 fun allpairs f xs ys =
   188   flat(map (fn x => map (fn y => f x y) ys) xs);
   189 
   190 fun extract_first p =
   191   let fun extract xs (y::ys) = if p y then (Some y,xs@ys)
   192                                else extract (y::xs) ys
   193         | extract xs []      = (None,xs)
   194   in extract [] end;
   195 
   196 
   197 fun print_ineqs ineqs =
   198  writeln(cat_lines(""::map (fn Lineq(c,t,l,_) =>
   199    string_of_int c ^
   200    (case t of Eq => " =  " | Lt=> " <  " | Le => " <= ") ^
   201    commas(map string_of_int l)) ineqs));
   202 
   203 
   204 fun elim ineqs =
   205   let (*val dummy = print_ineqs ineqs;*)
   206       val (triv,nontriv) = partition is_trivial ineqs in
   207   if not(null triv)
   208   then case find_first is_answer triv of
   209          None => elim nontriv | some => some
   210   else
   211   if null nontriv then None else
   212   let val (eqs,noneqs) = partition (fn (Lineq(_,ty,_,_)) => ty=Eq) nontriv in
   213   if not(null eqs) then
   214      let val clist = foldl (fn (cs,Lineq(_,_,l,_)) => l union cs) ([],eqs)
   215          val sclist = sort (fn (x,y) => int_ord(abs(x),abs(y)))
   216                            (filter (fn i => i<>0) clist)
   217          val c = hd sclist
   218          val (Some(eq as Lineq(_,_,ceq,_)),othereqs) =
   219                extract_first (fn Lineq(_,_,l,_) => c mem l) eqs
   220          val v = find_index (fn k => k=c) ceq
   221          val (ioth,roth) = partition (fn (Lineq(_,_,l,_)) => el v l = 0)
   222                                      (othereqs @ noneqs)
   223          val others = map (elim_var v eq) roth @ ioth
   224      in elim others end
   225   else
   226   let val lists = map (fn (Lineq(_,_,l,_)) => l) noneqs
   227       val numlist = 0 upto (length(hd lists) - 1)
   228       val coeffs = map (fn i => map (el i) lists) numlist
   229       val blows = map calc_blowup coeffs
   230       val iblows = blows ~~ numlist
   231       val nziblows = filter (fn (i,_) => i<>0) iblows
   232   in if null nziblows then None else
   233      let val (c,v) = hd(sort (fn (x,y) => int_ord(fst(x),fst(y))) nziblows)
   234          val (no,yes) = partition (fn (Lineq(_,_,l,_)) => el v l = 0) ineqs
   235          val (pos,neg) = partition(fn (Lineq(_,_,l,_)) => el v l > 0) yes
   236      in elim (no @ allpairs (elim_var v) pos neg) end
   237   end
   238   end
   239   end;
   240 
   241 (* ------------------------------------------------------------------------- *)
   242 (* Translate back a proof.                                                   *)
   243 (* ------------------------------------------------------------------------- *)
   244 
   245 (* FIXME OPTIMIZE!!!!
   246    Addition/Multiplication need i*t representation rather than t+t+...
   247 
   248 Simplification may detect a contradiction 'prematurely' due to type
   249 information: n+1 <= 0 is simplified to False and does not need to be crossed
   250 with 0 <= n.
   251 *)
   252 local
   253  exception FalseE of thm
   254 in
   255 fun mkthm sg asms just =
   256   let val atoms = foldl (fn (ats,(lhs,_,_,rhs,_,_)) =>
   257                             map fst lhs  union  (map fst rhs  union  ats))
   258                         ([], mapfilter (LA_Data.decomp o concl_of) asms)
   259 
   260       fun addthms thm1 thm2 =
   261         let val conj = thm1 RS (thm2 RS LA_Logic.conjI)
   262         in the(get_first (fn th => Some(conj RS th) handle _ => None)
   263                          (!LA_Data.add_mono_thms))
   264         end;
   265 
   266       fun multn(n,thm) =
   267         let fun mul(i,th) = if i=1 then th else mul(i-1, addthms thm th)
   268         in if n < 0 then mul(~n,thm) RS LA_Logic.sym else mul(n,thm) end;
   269 
   270       fun simp thm =
   271         let val thm' = !LA_Data.simp thm
   272         in if LA_Logic.is_False thm' then raise FalseE thm' else thm' end
   273 
   274       fun mk(Asm i) = ((*writeln"Asm";prth*)(nth_elem(i,asms)))
   275         | mk(Nat(i)) = ((*writeln"N";LA_Logic.mk_nat_thm sg (nth_elem(i,atoms)))
   276         | mk(LessD(j)) = ((*writeln"L";prth*)(hd([mk j] RL !LA_Data.lessD)))
   277         | mk(NotLeD(j)) = ((*writeln"NLe";prth*)(mk j RS LA_Logic.not_leD))
   278         | mk(NotLeDD(j)) = ((*writeln"NLeD";prth*)(hd([mk j RS LA_Logic.not_leD] RL !LA_Data.lessD)))
   279         | mk(NotLessD(j)) = ((*writeln"NL";prth*)(mk j RS LA_Logic.not_lessD))
   280         | mk(Added(j1,j2)) = ((*writeln"+";prth*)(simp(prth(addthms (mk j1) (mk j2)))))
   281         | mk(Multiplied(n,j)) = ((*writeln"*";multn(n,mk j))
   282 
   283   in (*writeln"mkthm";*)!LA_Data.simp(mk just) handle FalseE thm => thm end
   284 end;
   285 
   286 fun coeff poly atom = case assoc(poly,atom) of None => 0 | Some i => i;
   287 
   288 fun mklineq atoms =
   289   let val n = length atoms in
   290     fn ((lhs,i,rel,rhs,j,discrete),k) =>
   291     let val lhsa = map (coeff lhs) atoms
   292         and rhsa = map (coeff rhs) atoms
   293         val diff = map2 (op -) (rhsa,lhsa)
   294         val c = i-j
   295         val just = Asm k
   296     in case rel of
   297         "<="   => Some(Lineq(c,Le,diff,just))
   298        | "~<=" => if discrete
   299                   then Some(Lineq(1-c,Le,map (op ~) diff,NotLeDD(just)))
   300                   else Some(Lineq(~c,Lt,map (op ~) diff,NotLeD(just)))
   301        | "<"   => if discrete
   302                   then Some(Lineq(c+1,Le,diff,LessD(just)))
   303                   else Some(Lineq(c,Lt,diff,just))
   304        | "~<"  => Some(Lineq(~c,Le,map (op~) diff,NotLessD(just)))
   305        | "="   => Some(Lineq(c,Eq,diff,just))
   306        | "~="  => None
   307        | _     => sys_error("mklineq" ^ rel)   
   308     end
   309   end;
   310 
   311 fun mknat pTs ixs (atom,i) =
   312   if LA_Logic.is_nat(pTs,atom)
   313   then let val l = map (fn j => if j=i then 1 else 0) ixs
   314        in Some(Lineq(0,Le,l,Nat(i))) end
   315   else None
   316 
   317 fun abstract pTs items =
   318   let val atoms = foldl (fn (ats,((lhs,_,_,rhs,_,_),_)) =>
   319                             (map fst lhs) union ((map fst rhs) union ats))
   320                         ([],items)
   321       val ixs = 0 upto (length(atoms)-1)
   322       val iatoms = atoms ~~ ixs
   323   in mapfilter (mklineq atoms) items @ mapfilter (mknat pTs ixs) iatoms end;
   324 
   325 (* Ordinary refutation *)
   326 fun refute1(pTs,items) =
   327   (case elim (abstract pTs items) of
   328        None => []
   329      | Some(Lineq(_,_,_,j)) => [j]);
   330 
   331 fun refute1_tac(i,just) =
   332   fn state =>
   333     let val sg = #sign(rep_thm state)
   334     in resolve_tac [LA_Logic.notI,LA_Logic.ccontr] i THEN
   335        METAHYPS (fn asms => rtac (mkthm sg asms just) 1) i
   336     end
   337     state;
   338 
   339 (* Double refutation caused by equality in conclusion *)
   340 fun refute2(pTs,items, (rhs,i,_,lhs,j,d), nHs) =
   341   (case elim (abstract pTs (items@[((rhs,i,"<",lhs,j,d),nHs)])) of
   342     None => []
   343   | Some(Lineq(_,_,_,j1)) =>
   344       (case elim (abstract pTs (items@[((lhs,j,"<",rhs,i,d),nHs)])) of
   345         None => []
   346       | Some(Lineq(_,_,_,j2)) => [j1,j2]));
   347 
   348 fun refute2_tac(i,just1,just2) =
   349   fn state => 
   350     let val sg = #sign(rep_thm state)
   351     in rtac LA_Logic.ccontr i THEN rotate_tac ~1 i THEN etac LA_Logic.neqE i THEN
   352        METAHYPS (fn asms => rtac (mkthm sg asms just1) 1) i THEN
   353        METAHYPS (fn asms => rtac (mkthm sg asms just2) 1) i
   354     end
   355     state;
   356 
   357 fun prove(pTs,Hs,concl) =
   358 let val nHs = length Hs
   359     val ixHs = Hs ~~ (0 upto (nHs-1))
   360     val Hitems = mapfilter (fn (h,i) => case LA_Data.decomp h of
   361                                  None => None | Some(it) => Some(it,i)) ixHs
   362 in case LA_Data.decomp concl of
   363      None => if null Hitems then [] else refute1(pTs,Hitems)
   364    | Some(citem as (r,i,rel,l,j,d)) =>
   365        if rel = "="
   366        then refute2(pTs,Hitems,citem,nHs)
   367        else let val neg::rel0 = explode rel
   368                 val nrel = if neg = "~" then implode rel0 else "~"^rel
   369             in refute1(pTs, Hitems@[((r,i,nrel,l,j,d),nHs)]) end
   370 end;
   371 
   372 (*
   373 Fast but very incomplete decider. Only premises and conclusions
   374 that are already (negated) (in)equations are taken into account.
   375 *)
   376 val lin_arith_tac = SUBGOAL (fn (A,n) =>
   377   let val pTs = rev(map snd (Logic.strip_params A))
   378       val Hs = Logic.strip_assums_hyp A
   379       val concl = Logic.strip_assums_concl A
   380   in case prove(pTs,Hs,concl) of
   381        [j] => refute1_tac(n,j)
   382      | [j1,j2] => refute2_tac(n,j1,j2)
   383      | _ => no_tac
   384   end);
   385 
   386 fun cut_lin_arith_tac thms i = cut_facts_tac thms i THEN lin_arith_tac i;
   387 
   388 fun prover1(just,sg,thms,concl,pos) =
   389 let val nconcl = LA_Logic.neg_prop concl
   390     val cnconcl = cterm_of sg nconcl
   391     val Fthm = mkthm sg (thms @ [assume cnconcl]) just
   392     val contr = if pos then LA_Logic.ccontr else LA_Logic.notI
   393 in Some(LA_Logic.mk_Eq ((implies_intr cnconcl Fthm) COMP contr)) end
   394 handle _ => None;
   395 
   396 (* handle thm with equality conclusion *)
   397 fun prover2(just1,just2,sg,thms,concl) =
   398 let val nconcl = LA_Logic.neg_prop concl (* m ~= n *)
   399     val cnconcl = cterm_of sg nconcl
   400     val neqthm = assume cnconcl
   401     val casethm = neqthm COMP LA_Logic.neqE (* [|m<n ==> R; n<m ==> R|] ==> R *)
   402     val [lessimp1,lessimp2] = prems_of casethm
   403     val less1 = fst(Logic.dest_implies lessimp1) (* m<n *)
   404     and less2 = fst(Logic.dest_implies lessimp2) (* n<m *)
   405     val cless1 = cterm_of sg less1 and cless2 = cterm_of sg less2
   406     val thm1 = mkthm sg (thms @ [assume cless1]) just1
   407     and thm2 = mkthm sg (thms @ [assume cless2]) just2
   408     val dthm1 = implies_intr cless1 thm1 and dthm2 = implies_intr cless2 thm2
   409     val thm = dthm2 COMP (dthm1 COMP casethm)
   410 in Some(LA_Logic.mk_Eq ((implies_intr cnconcl thm) COMP LA_Logic.ccontr)) end
   411 handle _ => None;
   412 
   413 (* PRE: concl is not negated! *)
   414 fun lin_arith_prover sg thms concl =
   415 let val Hs = map (#prop o rep_thm) thms
   416     val Tconcl = LA_Logic.mk_Trueprop concl
   417 in case prove([],Hs,Tconcl) of
   418      [j] => prover1(j,sg,thms,Tconcl,true)
   419    | [j1,j2] => prover2(j1,j2,sg,thms,Tconcl)
   420    | _ => let val nTconcl = LA_Logic.neg_prop Tconcl
   421           in case prove([],Hs,nTconcl) of
   422                [j] => prover1(j,sg,thms,nTconcl,false)
   423                (* [_,_] impossible because of negation *)
   424              | _ => None
   425           end
   426 end;
   427 
   428 end;