Complex/complex_grobner.ml

   1 (* ========================================================================= *)
   2 (* Grobner basis algorithm.                                                  *)
   3 (* ========================================================================= *)
   4
   5 needs "Complex/complexnumbers.ml";;
   6 needs "Complex/quelim.ml";;
   7
   8 prioritize_complex();;
   9
  10 (* ------------------------------------------------------------------------- *)
  11 (* Utility functions.                                                        *)
  12 (* ------------------------------------------------------------------------- *)
  13
  14 let allpairs f l1 l2 =
  15   itlist ((@) o C map l2 o f) l1 [];;
  16
  17 let rec merge ord l1 l2 =
  18   match l1 with
  19     [] -> l2
  20   | h1::t1 -> match l2 with
  21                 [] -> l1
  22               | h2::t2 -> if ord h1 h2 then h1::(merge ord t1 l2)
  23                           else h2::(merge ord l1 t2);;
  24
  25 let sort ord =
  26   let rec mergepairs l1 l2 =
  27     match (l1,l2) with
  28         ([s],[]) -> s
  29       | (l,[]) -> mergepairs [] l
  30       | (l,[s1]) -> mergepairs (s1::l) []
  31       | (l,(s1::s2::ss)) -> mergepairs ((merge ord s1 s2)::l) ss in
  32   fun l -> if l = [] then [] else mergepairs [] (map (fun x -> [x]) l);;
  33
  34 (* ------------------------------------------------------------------------- *)
  35 (* Type for recording history, i.e. how a polynomial was obtained.           *)
  36 (* ------------------------------------------------------------------------- *)
  37
  38 type history =
  39    Start of int
  40  | Mmul of (num * (int list)) * history
  41  | Add of history * history;;
  42
  43 (* ------------------------------------------------------------------------- *)
  44 (* Conversion of leaves, i.e. variables and constant rational constants.     *)
  45 (* ------------------------------------------------------------------------- *)
  46
  47 let grob_var vars tm =
  48   let res = map (fun i -> if i = tm then 1 else 0) vars in
  49   if exists (fun x -> x <> 0) res then [Int 1,res] else failwith "grob_var";;
  50
  51 let grob_const =
  52   let cx_tm = `Cx` in
  53   fun vars tm ->
  54     try let l,r = dest_comb tm in
  55         if l = cx_tm then
  56           let x = rat_of_term r in
  57           if x =/ Int 0 then [] else [x,map (fun v -> 0) vars]
  58         else failwith ""
  59     with Failure _ -> failwith "grob_const";;
  60
  61 (* ------------------------------------------------------------------------- *)
  62 (* Monomial ordering.                                                        *)
  63 (* ------------------------------------------------------------------------- *)
  64
  65 let morder_lt =
  66   let rec lexorder l1 l2 =
  67     match (l1,l2) with
  68         [],[] -> false
  69       | (x1::o1,x2::o2) -> x1 > x2 or x1 = x2 & lexorder o1 o2
  70       | _ -> failwith "morder: inconsistent monomial lengths" in
  71   fun m1 m2 -> let n1 = itlist (+) m1 0
  72                and n2 = itlist (+) m2 0 in
  73                n1 < n2 or n1 = n2 & lexorder m1 m2;;
  74
  75 let morder_le m1 m2 = morder_lt m1 m2 or (m1 = m2);;
  76
  77 let morder_gt m1 m2 = morder_lt m2 m1;;
  78
  79 (* ------------------------------------------------------------------------- *)
  80 (* Arithmetic on canonical polynomials.                                      *)
  81 (* ------------------------------------------------------------------------- *)
  82
  83 let grob_neg = map (fun (c,m) -> (minus_num c,m));;
  84
  85 let rec grob_add l1 l2 =
  86   match (l1,l2) with
  87     ([],l2) -> l2
  88   | (l1,[]) -> l1
  89   | ((c1,m1)::o1,(c2,m2)::o2) ->
  90         if m1 = m2 then
  91           let c = c1+/c2 and rest = grob_add o1 o2 in
  92           if c =/ Int 0 then rest else (c,m1)::rest
  93         else if morder_lt m2 m1 then (c1,m1)::(grob_add o1 l2)
  94         else (c2,m2)::(grob_add l1 o2);;
  95
  96 let grob_sub l1 l2 = grob_add l1 (grob_neg l2);;
  97
  98 let grob_mmul (c1,m1) (c2,m2) = (c1*/c2,map2 (+) m1 m2);;
  99
 100 let rec grob_cmul cm pol = map (grob_mmul cm) pol;;
 101
 102 let rec grob_mul l1 l2 =
 103   match l1 with
 104     [] -> []
 105   | (h1::t1) -> grob_add (grob_cmul h1 l2) (grob_mul t1 l2);;
 106
 107 let rec grob_pow vars l n =
 108   if n < 0 then failwith "grob_pow: negative power"
 109   else if n = 0 then [Int 1,map (fun v -> 0) vars]
 110   else grob_mul l (grob_pow vars l (n - 1));;
 111
 112 (* ------------------------------------------------------------------------- *)
 113 (* Monomial division operation.                                              *)
 114 (* ------------------------------------------------------------------------- *)
 115
 116 let mdiv (c1,m1) (c2,m2) =
 117   (c1//c2,
 118    map2 (fun n1 n2 -> if n1 < n2 then failwith "mdiv" else n1-n2) m1 m2);;
 119
 120 (* ------------------------------------------------------------------------- *)
 121 (* Lowest common multiple of two monomials.                                  *)
 122 (* ------------------------------------------------------------------------- *)
 123
 124 let mlcm (c1,m1) (c2,m2) = (Int 1,map2 max m1 m2);;
 125
 126 (* ------------------------------------------------------------------------- *)
 127 (* Reduce monomial cm by polynomial pol, returning replacement for cm.       *)
 128 (* ------------------------------------------------------------------------- *)
 129
 130 let reduce1 cm (pol,hpol) =
 131   match pol with
 132     [] -> failwith "reduce1"
 133   | cm1::cms -> try let (c,m) = mdiv cm cm1 in
 134                     (grob_cmul (minus_num c,m) cms,
 135                      Mmul((minus_num c,m),hpol))
 136                 with Failure _ -> failwith "reduce1";;
 137
 138 (* ------------------------------------------------------------------------- *)
 139 (* Try this for all polynomials in a basis.                                  *)
 140 (* ------------------------------------------------------------------------- *)
 141
 142 let reduceb cm basis = tryfind (fun p -> reduce1 cm p) basis;;
 143
 144 (* ------------------------------------------------------------------------- *)
 145 (* Reduction of a polynomial (always picking largest monomial possible).     *)
 146 (* ------------------------------------------------------------------------- *)
 147
 148 let rec reduce basis (pol,hist) =
 149   match pol with
 150     [] -> (pol,hist)
 151   | cm::ptl -> try let q,hnew = reduceb cm basis in
 152                    reduce basis (grob_add q ptl,Add(hnew,hist))
 153                with Failure _ ->
 154                    let q,hist' = reduce basis (ptl,hist) in
 155                    cm::q,hist';;
 156
 157 (* ------------------------------------------------------------------------- *)
 158 (* Check for orthogonality w.r.t. LCM.                                       *)
 159 (* ------------------------------------------------------------------------- *)
 160
 161 let orthogonal l p1 p2 =
 162   snd l = snd(grob_mmul (hd p1) (hd p2));;
 163
 164 (* ------------------------------------------------------------------------- *)
 165 (* Compute S-polynomial of two polynomials.                                  *)
 166 (* ------------------------------------------------------------------------- *)
 167
 168 let spoly cm ph1 ph2 =
 169   match (ph1,ph2) with
 170     ([],h),p -> ([],h)
 171   | p,([],h) -> ([],h)
 172   | (cm1::ptl1,his1),(cm2::ptl2,his2) ->
 173         (grob_sub (grob_cmul (mdiv cm cm1) ptl1)
 174                   (grob_cmul (mdiv cm cm2) ptl2),
 175          Add(Mmul(mdiv cm cm1,his1),
 176              Mmul(mdiv (minus_num(fst cm),snd cm) cm2,his2)));;
 177
 178 (* ------------------------------------------------------------------------- *)
 179 (* Make a polynomial monic.                                                  *)
 180 (* ------------------------------------------------------------------------- *)
 181
 182 let monic (pol,hist) =
 183   if pol = [] then (pol,hist) else
 184   let c',m' = hd pol in
 185   (map (fun (c,m) -> (c//c',m)) pol,
 186    Mmul((Int 1 // c',map (K 0) m'),hist));;
 187
 188 (* ------------------------------------------------------------------------- *)
 189 (* The most popular heuristic is to order critical pairs by LCM monomial.    *)
 190 (* ------------------------------------------------------------------------- *)
 191
 192 let forder ((c1,m1),_) ((c2,m2),_) = morder_lt m1 m2;;
 193
 194 (* ------------------------------------------------------------------------- *)
 195 (* Stupid stuff forced on us by lack of equality test on num type.           *)
 196 (* ------------------------------------------------------------------------- *)
 197
 198 let rec poly_lt p q =
 199   match (p,q) with
 200     p,[] -> false
 201   | [],q -> true
 202   | (c1,m1)::o1,(c2,m2)::o2 ->
 203         c1 </ c2 or
 204         c1 =/ c2 & (m1 < m2 or m1 = m2 & poly_lt o1 o2);;
 205
 206 let align ((p,hp),(q,hq)) =
 207   if poly_lt p q then ((p,hp),(q,hq)) else ((q,hq),(p,hp));;
 208
 209 let poly_eq p1 p2 =
 210   forall2 (fun (c1,m1) (c2,m2) -> c1 =/ c2 & m1 = m2) p1 p2;;
 211
 212 let memx ((p1,h1),(p2,h2)) ppairs =
 213   not (exists (fun ((q1,_),(q2,_)) -> poly_eq p1 q1 & poly_eq p2 q2) ppairs);;
 214
 215 (* ------------------------------------------------------------------------- *)
 216 (* Buchberger's second criterion.                                            *)
 217 (* ------------------------------------------------------------------------- *)
 218
 219 let criterion2 basis (lcm,((p1,h1),(p2,h2))) opairs =
 220   exists (fun g -> not(poly_eq (fst g) p1) & not(poly_eq (fst g) p2) &
 221                    can (mdiv lcm) (hd(fst g)) &
 222                    not(memx (align(g,(p1,h1))) (map snd opairs)) &
 223                    not(memx (align(g,(p2,h2))) (map snd opairs))) basis;;
 224
 225 (* ------------------------------------------------------------------------- *)
 226 (* Test for hitting constant polynomial.                                     *)
 227 (* ------------------------------------------------------------------------- *)
 228
 229 let constant_poly p =
 230   length p = 1 & forall ((=) 0) (snd(hd p));;
 231
 232 (* ------------------------------------------------------------------------- *)
 233 (* Grobner basis algorithm.                                                  *)
 234 (* ------------------------------------------------------------------------- *)
 235
 236 let rec grobner histories basis pairs =
 237   print_string(string_of_int(length basis)^" basis elements and "^
 238                string_of_int(length pairs)^" critical pairs");
 239   print_newline();
 240   match pairs with
 241     [] -> rev histories,basis
 242   | (l,(p1,p2))::opairs ->
 243         let (sp,hist) = monic (reduce basis (spoly l p1 p2)) in
 244         if sp = [] or criterion2 basis (l,(p1,p2)) opairs
 245         then grobner histories basis opairs else
 246         let sph = sp,Start(length histories) in
 247         if constant_poly sp
 248         then grobner ((sp,hist)::histories) (sph::basis) [] else
 249         let rawcps =
 250           map (fun p -> mlcm (hd(fst p)) (hd sp),align(p,sph)) basis in
 251         let newcps = filter
 252           (fun (l,(p,q)) -> not(orthogonal l (fst p) (fst q))) rawcps in
 253         grobner ((sp,hist)::histories) (sph::basis)
 254                 (merge forder opairs (sort forder newcps));;
 255
 256 (* ------------------------------------------------------------------------- *)
 257 (* Overall function.                                                         *)
 258 (* ------------------------------------------------------------------------- *)
 259
 260 let groebner pols =
 261   let npols = map2 (fun p n -> p,Start n) pols (0--(length pols - 1)) in
 262   let phists = filter (fun (p,_) -> p <> []) npols in
 263   let bas0 = map monic phists in
 264   let bas = map2 (fun (p,h) n -> p,Start n) bas0
 265                  ((length bas0)--(2 * length bas0 - 1)) in
 266   let prs0 = allpairs (fun x y -> x,y) bas bas in
 267   let prs1 = filter (fun ((x,_),(y,_)) -> poly_lt x y) prs0 in
 268   let prs2 = map (fun (p,q) -> mlcm (hd(fst p)) (hd(fst q)),(p,q)) prs1 in
 269   let prs3 = filter (fun (l,(p,q)) -> not(orthogonal l (fst p) (fst q))) prs2 in
 270   grobner (rev bas0 @ rev phists) bas (mergesort forder prs3);;
 271
 272 (* ------------------------------------------------------------------------- *)
 273 (* Alternative orthography.                                                  *)
 274 (* ------------------------------------------------------------------------- *)
 275
 276 let gr'o'bner = groebner;;
 277
 278 (* ------------------------------------------------------------------------- *)
 279 (* Conversion from HOL term.                                                 *)
 280 (* ------------------------------------------------------------------------- *)
 281
 282 let grobify_term =
 283   let neg_tm = `(--):complex->complex`
 284   and add_tm = `(+):complex->complex->complex`
 285   and sub_tm = `(-):complex->complex->complex`
 286   and mul_tm = `(*):complex->complex->complex`
 287   and pow_tm = `(pow):complex->num->complex` in
 288   let rec grobify_term vars tm =
 289     try grob_var vars tm with Failure _ ->
 290     try grob_const vars tm with Failure _ ->
 291     let lop,r = dest_comb tm in
 292     if lop = neg_tm then grob_neg(grobify_term vars r) else
 293     let op,l = dest_comb lop in
 294     if op = pow_tm then
 295       grob_pow vars (grobify_term vars l) (dest_small_numeral r)
 296     else
 297      (if op = add_tm then grob_add else if op = sub_tm then grob_sub
 298       else if op = mul_tm then grob_mul else failwith "unknown term")
 299      (grobify_term vars l) (grobify_term vars r) in
 300   fun vars tm ->
 301     try grobify_term vars tm with Failure _ -> failwith "grobify_term";;
 302
 303 let grobvars =
 304   let neg_tm = `(--):complex->complex`
 305   and add_tm = `(+):complex->complex->complex`
 306   and sub_tm = `(-):complex->complex->complex`
 307   and mul_tm = `(*):complex->complex->complex`
 308   and pow_tm = `(pow):complex->num->complex` in
 309   let rec grobvars tm acc =
 310     if is_complex_const tm then acc
 311     else if not (is_comb tm) then tm::acc else
 312     let lop,r = dest_comb tm in
 313     if lop = neg_tm then grobvars r acc
 314     else if not (is_comb lop) then tm::acc else
 315     let op,l = dest_comb lop in
 316     if op = pow_tm then grobvars l acc
 317     else if op = mul_tm or op = sub_tm or op = add_tm
 318     then grobvars l (grobvars r acc)
 319     else tm::acc in
 320   fun tm -> setify(grobvars tm []);;
 321
 322 let grobify_equations =
 323   let zero_tm = `Cx(&0)`
 324   and sub_tm = `(-):complex->complex->complex`
 325   and complex_ty = `:complex` in
 326   let grobify_equation vars tm =
 327   let l,r = dest_eq tm in
 328     if r <> zero_tm then grobify_term vars (mk_comb(mk_comb(sub_tm,l),r))
 329     else grobify_term vars l in
 330   fun tm ->
 331     let cjs = conjuncts tm in
 332     let rawvars = itlist
 333        (fun eq acc -> let l,r = dest_eq eq in
 334                      union (union (grobvars l) (grobvars r)) acc) cjs [] in
 335     let vars = sort (fun x y -> x < y) rawvars in
 336     vars,map(grobify_equation vars) cjs;;
 337
 338 (* ------------------------------------------------------------------------- *)
 339 (* Printer.                                                                  *)
 340 (* ------------------------------------------------------------------------- *)
 341
 342 let string_of_monomial vars (c,m) =
 343   let xns = filter (fun (x,y) -> y <> 0) (map2 (fun x y -> x,y) vars m) in
 344   let xnstrs = map
 345     (fun (x,n) -> x^(if n = 1 then "" else "^"^(string_of_int n))) xns in
 346   if xns = [] then Num.string_of_num c else
 347   let basstr = if c =/ Int 1 then "" else (Num.string_of_num c)^" * " in
 348   basstr ^ end_itlist (fun s t -> s^" * "^t) xnstrs;;
 349
 350 let string_of_polynomial vars l =
 351   if l = [] then "0" else
 352   end_itlist (fun s t -> s^" + "^t) (map (string_of_monomial vars) l);;
 353
 354 (* ------------------------------------------------------------------------- *)
 355 (* Resolve a proof.                                                          *)
 356 (* ------------------------------------------------------------------------- *)
 357
 358 let rec resolve_proof vars prf =
 359   match prf with
 360     Start n ->
 361         [n,[Int 1,map (K 0) vars]]
 362   | Mmul(pol,lin) ->
 363         let lis = resolve_proof vars lin in
 364         map (fun (n,p) -> n,grob_cmul pol p) lis
 365   | Add(lin1,lin2) ->
 366         let lis1 = resolve_proof vars lin1
 367         and lis2 = resolve_proof vars lin2 in
 368         let dom = setify(union (map fst lis1) (map fst lis2)) in
 369         map (fun n -> let a = try assoc n lis1 with Failure _ -> []
 370                       and b = try assoc n lis2 with Failure _ -> [] in
 371                       n,grob_add a b) dom;;
 372
 373 (* ------------------------------------------------------------------------- *)
 374 (* Convert a polynomial back to HOL.                                         *)
 375 (* ------------------------------------------------------------------------- *)
 376
 377 let holify_polynomial =
 378   let complex_ty = `:complex`
 379   and pow_tm = `(pow):complex->num->complex`
 380   and mk_mul = mk_binop `(*):complex->complex->complex`
 381   and mk_add = mk_binop `(+):complex->complex->complex`
 382   and zero_tm = `Cx(&0)`
 383   and add_tm = `(+):complex->complex->complex`
 384   and mul_tm = `(*):complex->complex->complex`
 385   and complex_term_of_rat = curry mk_comb `Cx` o term_of_rat in
 386   let holify_varpow (v,n) =
 387     if n = 1 then v
 388     else list_mk_comb(pow_tm,[v;mk_small_numeral n]) in
 389   let holify_monomial vars (c,m) =
 390     let xps = map holify_varpow (filter (fun (_,n) -> n <> 0) (zip vars m)) in
 391     end_itlist mk_mul (complex_term_of_rat c :: xps) in
 392   let holify_polynomial vars p =
 393     if p = [] then zero_tm
 394     else end_itlist mk_add (map (holify_monomial vars) p) in
 395   holify_polynomial;;
 396
 397 (* ------------------------------------------------------------------------- *)
 398 (* Recursively find the set of basis elements involved.                      *)
 399 (* ------------------------------------------------------------------------- *)
 400
 401 let dependencies =
 402   let rec dependencies prf acc =
 403     match prf with
 404       Start n -> n::acc
 405     | Mmul(pol,lin) -> dependencies lin acc
 406     | Add(lin1,lin2) -> dependencies lin1 (dependencies lin2 acc) in
 407   fun prf -> setify(dependencies prf []);;
 408
 409 let rec involved deps sofar todo =
 410   match todo with
 411     [] -> sort (<) sofar
 412   | (h::hs) ->
 413         if mem h sofar then involved deps sofar hs
 414         else involved deps (h::sofar) (el h deps @ hs);;
 415
 416 (* ------------------------------------------------------------------------- *)
 417 (* Refute a conjunction of equations in HOL.                                 *)
 418 (* ------------------------------------------------------------------------- *)
 419
 420 let GROBNER_REFUTE =
 421   let add_tm = `(+):complex->complex->complex`
 422   and mul_tm = `(*):complex->complex->complex` in
 423   let APPLY_pth = MATCH_MP(SIMPLE_COMPLEX_ARITH
 424     `(x = y) ==> (x + Cx(--(&1)) * (y + Cx(&1)) = Cx(&0)) ==> F`)
 425   and MK_ADD th1 th2 = MK_COMB(AP_TERM add_tm th1,th2) in
 426   let rec holify_lincombs vars cjs prfs =
 427     match prfs with
 428       [] -> cjs
 429     | (p::ps) ->
 430         if p = [] then holify_lincombs vars (cjs @ [TRUTH]) ps else
 431         let holis = map (fun (n,q) -> n,holify_polynomial vars q) p in
 432         let ths =
 433           map (fun (n,m) -> AP_TERM (mk_comb(mul_tm,m)) (el n cjs)) holis in
 434         let th = CONV_RULE(BINOP_CONV COMPLEX_POLY_CONV)
 435                           (end_itlist MK_ADD ths) in
 436         holify_lincombs vars (cjs @ [th]) ps in
 437   fun tm ->
 438     let vars,pols = grobify_equations tm in
 439     let (prfs,gb) = groebner pols in
 440     let (_,prf) =
 441       find (fun (p,h) -> length p = 1 & forall ((=)0) (snd(hd p))) gb in
 442     let deps = map (dependencies o snd) prfs
 443     and depl = dependencies prf in
 444     let need = involved deps [] depl in
 445     let pprfs =
 446       map2 (fun p n -> if mem n need then resolve_proof vars (snd p) else [])
 447            prfs (0--(length prfs - 1))
 448     and ppr = resolve_proof vars prf in
 449     let ths = CONJUNCTS(ASSUME tm) in
 450     let th = last
 451      (holify_lincombs vars ths (snd(chop_list(length ths) (pprfs @ [ppr])))) in
 452     CONV_RULE COMPLEX_RAT_EQ_CONV th;;
 453
 454 (* ------------------------------------------------------------------------- *)
 455 (* Overall conversion.                                                       *)
 456 (* ------------------------------------------------------------------------- *)
 457
 458 let COMPLEX_ARITH =
 459   let pth0 = SIMPLE_COMPLEX_ARITH `(x = y) <=> (x - y = Cx(&0))`
 460   and pth1 = prove
 461    (`!x. ~(x = Cx(&0)) <=> ?z. z * x + Cx(&1) = Cx(&0)`,
 462     CONV_TAC(time FULL_COMPLEX_QUELIM_CONV))
 463   and pth2a = prove
 464    (`!x y u v. ~(x = y) \/ ~(u = v) <=>
 465                 ?w z. w * (x - y) + z * (u - v) + Cx(&1) = Cx(&0)`,
 466     CONV_TAC(time FULL_COMPLEX_QUELIM_CONV))
 467   and pth2b = prove
 468    (`!x y. ~(x = y) <=> ?z. z * (x - y) + Cx(&1) = Cx(&0)`,
 469     CONV_TAC(time FULL_COMPLEX_QUELIM_CONV))
 470   and pth3 = TAUT `(p ==> F) ==> (~q <=> p) ==> q` in
 471   let GEN_PRENEX_CONV =
 472     GEN_REWRITE_CONV REDEPTH_CONV
 473      [AND_FORALL_THM;
 474       LEFT_AND_FORALL_THM;
 475       RIGHT_AND_FORALL_THM;
 476       LEFT_OR_FORALL_THM;
 477       RIGHT_OR_FORALL_THM;
 478       OR_EXISTS_THM;
 479       LEFT_OR_EXISTS_THM;
 480       RIGHT_OR_EXISTS_THM;
 481       LEFT_AND_EXISTS_THM;
 482       RIGHT_AND_EXISTS_THM] in
 483   let INITIAL_CONV =
 484     NNF_CONV THENC
 485     GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth1] THENC
 486     GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth2a] THENC
 487     GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth2b] THENC
 488     ONCE_DEPTH_CONV(GEN_REWRITE_CONV I [pth0] o
 489                     check ((<>) `Cx(&0)` o rand)) THENC
 490     GEN_PRENEX_CONV THENC
 491     DNF_CONV in
 492   fun tm ->
 493     let avs = frees tm in
 494     let tm' = list_mk_forall(avs,tm) in
 495     let th1 = INITIAL_CONV(mk_neg tm') in
 496     let evs,bod = strip_exists(rand(concl th1)) in
 497     if is_forall bod then failwith "COMPLEX_ARITH: non-universal formula" else
 498     let djs = disjuncts bod in
 499     let th2 = end_itlist SIMPLE_DISJ_CASES(map GROBNER_REFUTE djs) in
 500     let th3 = itlist SIMPLE_CHOOSE evs th2 in
 501     SPECL avs (MATCH_MP (MATCH_MP pth3 (DISCH_ALL th3)) th1);;