grobner.ml

   1 (* ========================================================================= *)
   2 (* Generic Grobner basis algorithm.                                          *)
   3 (*                                                                           *)
   4 (* Whatever the instantiation, it basically solves the universal theory of   *)
   5 (* the complex numbers, or equivalently something like the theory of all     *)
   6 (* commutative cancellation semirings with no nilpotent elements and having  *)
   7 (* characteristic zero. We could do "all rings" by a more elaborate integer  *)
   8 (* version of Grobner bases, but I don't have any useful applications.       *)
   9 (*                                                                           *)
  10 (*              (c) Copyright, John Harrison 1998-2007                       *)
  11 (* ========================================================================= *)
  12
  13 needs "normalizer.ml";;
  14
  15 (* ------------------------------------------------------------------------- *)
  16 (* Type for recording history, i.e. how a polynomial was obtained.           *)
  17 (* ------------------------------------------------------------------------- *)
  18
  19 type history =
  20    Start of int
  21  | Mmul of (num * (int list)) * history
  22  | Add of history * history;;
  23
  24 (* ------------------------------------------------------------------------- *)
  25 (* Overall function; everything else is local.                               *)
  26 (* ------------------------------------------------------------------------- *)
  27
  28 let RING_AND_IDEAL_CONV =
  29
  30   (* ----------------------------------------------------------------------- *)
  31   (* Monomial ordering.                                                      *)
  32   (* ----------------------------------------------------------------------- *)
  33
  34   let morder_lt =
  35     let rec lexorder l1 l2 =
  36       match (l1,l2) with
  37           [],[] -> false
  38         | (x1::o1,x2::o2) -> x1 > x2 or x1 = x2 & lexorder o1 o2
  39         | _ -> failwith "morder: inconsistent monomial lengths" in
  40     fun m1 m2 -> let n1 = itlist (+) m1 0
  41                  and n2 = itlist (+) m2 0 in
  42                  n1 < n2 or n1 = n2 & lexorder m1 m2 in
  43
  44   (* ----------------------------------------------------------------------- *)
  45   (* Arithmetic on canonical polynomials.                                    *)
  46   (* ----------------------------------------------------------------------- *)
  47
  48   let grob_neg = map (fun (c,m) -> (minus_num c,m)) in
  49
  50   let rec grob_add l1 l2 =
  51     match (l1,l2) with
  52       ([],l2) -> l2
  53     | (l1,[]) -> l1
  54     | ((c1,m1)::o1,(c2,m2)::o2) ->
  55           if m1 = m2 then
  56             let c = c1+/c2 and rest = grob_add o1 o2 in
  57             if c =/ num_0 then rest else (c,m1)::rest
  58           else if morder_lt m2 m1 then (c1,m1)::(grob_add o1 l2)
  59           else (c2,m2)::(grob_add l1 o2) in
  60
  61   let grob_sub l1 l2 = grob_add l1 (grob_neg l2) in
  62
  63   let grob_mmul (c1,m1) (c2,m2) = (c1*/c2,map2 (+) m1 m2) in
  64
  65   let rec grob_cmul cm pol = map (grob_mmul cm) pol in
  66
  67   let rec grob_mul l1 l2 =
  68     match l1 with
  69       [] -> []
  70     | (h1::t1) -> grob_add (grob_cmul h1 l2) (grob_mul t1 l2) in
  71
  72   let grob_inv l =
  73     match l with
  74       [c,vs] when forall (fun x -> x = 0) vs ->
  75           if c =/ num_0 then failwith "grob_inv: division by zero"
  76           else [num_1 // c,vs]
  77     | _ -> failwith "grob_inv: non-constant divisor polynomial" in
  78
  79   let grob_div l1 l2 =
  80     match l2 with
  81       [c,l] when forall (fun x -> x = 0) l ->
  82           if c =/ num_0 then failwith "grob_div: division by zero"
  83           else grob_cmul (num_1 // c,l) l1
  84     | _ -> failwith "grob_div: non-constant divisor polynomial" in
  85
  86   let rec grob_pow vars l n =
  87     if n < 0 then failwith "grob_pow: negative power"
  88     else if n = 0 then [num_1,map (fun v -> 0) vars]
  89     else grob_mul l (grob_pow vars l (n - 1)) in
  90
  91   (* ----------------------------------------------------------------------- *)
  92   (* Monomial division operation.                                            *)
  93   (* ----------------------------------------------------------------------- *)
  94
  95   let mdiv (c1,m1) (c2,m2) =
  96     (c1//c2,
  97      map2 (fun n1 n2 -> if n1 < n2 then failwith "mdiv" else n1-n2) m1 m2) in
  98
  99   (* ----------------------------------------------------------------------- *)
 100   (* Lowest common multiple of two monomials.                                *)
 101   (* ----------------------------------------------------------------------- *)
 102
 103   let mlcm (c1,m1) (c2,m2) = (num_1,map2 max m1 m2) in
 104
 105   (* ----------------------------------------------------------------------- *)
 106   (* Reduce monomial cm by polynomial pol, returning replacement for cm.     *)
 107   (* ----------------------------------------------------------------------- *)
 108
 109   let reduce1 cm (pol,hpol) =
 110     match pol with
 111       [] -> failwith "reduce1"
 112     | cm1::cms -> try let (c,m) = mdiv cm cm1 in
 113                       (grob_cmul (minus_num c,m) cms,
 114                        Mmul((minus_num c,m),hpol))
 115                   with Failure _ -> failwith "reduce1" in
 116
 117   (* ----------------------------------------------------------------------- *)
 118   (* Try this for all polynomials in a basis.                                *)
 119   (* ----------------------------------------------------------------------- *)
 120
 121   let reduceb cm basis = tryfind (fun p -> reduce1 cm p) basis in
 122
 123   (* ----------------------------------------------------------------------- *)
 124   (* Reduction of a polynomial (always picking largest monomial possible).   *)
 125   (* ----------------------------------------------------------------------- *)
 126
 127   let rec reduce basis (pol,hist) =
 128     match pol with
 129       [] -> (pol,hist)
 130     | cm::ptl -> try let q,hnew = reduceb cm basis in
 131                      reduce basis (grob_add q ptl,Add(hnew,hist))
 132                  with Failure _ ->
 133                      let q,hist' = reduce basis (ptl,hist) in
 134                      cm::q,hist' in
 135
 136   (* ----------------------------------------------------------------------- *)
 137   (* Check for orthogonality w.r.t. LCM.                                     *)
 138   (* ----------------------------------------------------------------------- *)
 139
 140   let orthogonal l p1 p2 =
 141     snd l = snd(grob_mmul (hd p1) (hd p2)) in
 142
 143   (* ----------------------------------------------------------------------- *)
 144   (* Compute S-polynomial of two polynomials.                                *)
 145   (* ----------------------------------------------------------------------- *)
 146
 147   let spoly cm ph1 ph2 =
 148     match (ph1,ph2) with
 149       ([],h),p -> ([],h)
 150     | p,([],h) -> ([],h)
 151     | (cm1::ptl1,his1),(cm2::ptl2,his2) ->
 152           (grob_sub (grob_cmul (mdiv cm cm1) ptl1)
 153                     (grob_cmul (mdiv cm cm2) ptl2),
 154            Add(Mmul(mdiv cm cm1,his1),
 155                Mmul(mdiv (minus_num(fst cm),snd cm) cm2,his2))) in
 156
 157   (* ----------------------------------------------------------------------- *)
 158   (* Make a polynomial monic.                                                *)
 159   (* ----------------------------------------------------------------------- *)
 160
 161   let monic (pol,hist) =
 162     if pol = [] then (pol,hist) else
 163     let c',m' = hd pol in
 164     (map (fun (c,m) -> (c//c',m)) pol,
 165      Mmul((num_1 // c',map (K 0) m'),hist)) in
 166
 167   (* ----------------------------------------------------------------------- *)
 168   (* The most popular heuristic is to order critical pairs by LCM monomial.  *)
 169   (* ----------------------------------------------------------------------- *)
 170
 171   let forder ((c1,m1),_) ((c2,m2),_) = morder_lt m1 m2 in
 172
 173   (* ----------------------------------------------------------------------- *)
 174   (* Stupid stuff forced on us by lack of equality test on num type.         *)
 175   (* ----------------------------------------------------------------------- *)
 176
 177   let rec poly_lt p q =
 178     match (p,q) with
 179       p,[] -> false
 180     | [],q -> true
 181     | (c1,m1)::o1,(c2,m2)::o2 ->
 182           c1 </ c2 or
 183           c1 =/ c2 & (m1 < m2 or m1 = m2 & poly_lt o1 o2) in
 184
 185   let align ((p,hp),(q,hq)) =
 186     if poly_lt p q then ((p,hp),(q,hq)) else ((q,hq),(p,hp)) in
 187
 188   let poly_eq p1 p2 =
 189     forall2 (fun (c1,m1) (c2,m2) -> c1 =/ c2 & m1 = m2) p1 p2 in
 190
 191   let memx ((p1,h1),(p2,h2)) ppairs =
 192     not (exists (fun ((q1,_),(q2,_)) -> poly_eq p1 q1 & poly_eq p2 q2)
 193                 ppairs) in
 194
 195   (* ----------------------------------------------------------------------- *)
 196   (* Buchberger's second criterion.                                          *)
 197   (* ----------------------------------------------------------------------- *)
 198
 199   let criterion2 basis (lcm,((p1,h1),(p2,h2))) opairs =
 200     exists (fun g -> not(poly_eq (fst g) p1) & not(poly_eq (fst g) p2) &
 201                      can (mdiv lcm) (hd(fst g)) &
 202                      not(memx (align(g,(p1,h1))) (map snd opairs)) &
 203                      not(memx (align(g,(p2,h2))) (map snd opairs))) basis in
 204
 205   (* ----------------------------------------------------------------------- *)
 206   (* Test for hitting constant polynomial.                                   *)
 207   (* ----------------------------------------------------------------------- *)
 208
 209   let constant_poly p =
 210     length p = 1 & forall ((=) 0) (snd(hd p)) in
 211
 212   (* ----------------------------------------------------------------------- *)
 213   (* Grobner basis algorithm.                                                *)
 214   (* ----------------------------------------------------------------------- *)
 215
 216   let rec grobner_basis basis pairs =
 217     Format.print_string(string_of_int(length basis)^" basis elements and "^
 218                         string_of_int(length pairs)^" critical pairs");
 219     Format.print_newline();
 220     match pairs with
 221       [] -> basis
 222     | (l,(p1,p2))::opairs ->
 223           let (sp,hist as sph) = monic (reduce basis (spoly l p1 p2)) in
 224           if sp = [] or criterion2 basis (l,(p1,p2)) opairs
 225           then grobner_basis basis opairs else
 226           if constant_poly sp then grobner_basis (sph::basis) [] else
 227           let rawcps =
 228             map (fun p -> mlcm (hd(fst p)) (hd sp),align(p,sph)) basis in
 229           let newcps = filter
 230             (fun (l,(p,q)) -> not(orthogonal l (fst p) (fst q))) rawcps in
 231           grobner_basis (sph::basis)
 232                   (merge forder opairs (mergesort forder newcps)) in
 233
 234   (* ----------------------------------------------------------------------- *)
 235   (* Interreduce initial polynomials.                                        *)
 236   (* ----------------------------------------------------------------------- *)
 237
 238   let rec grobner_interreduce rpols ipols =
 239     match ipols with
 240       [] -> map monic (rev rpols)
 241     | p::ps -> let p' = reduce (rpols @ ps) p in
 242                if fst p' = [] then grobner_interreduce rpols ps
 243                else grobner_interreduce (p'::rpols) ps in
 244
 245   (* ----------------------------------------------------------------------- *)
 246   (* Overall function.                                                       *)
 247   (* ----------------------------------------------------------------------- *)
 248
 249   let grobner pols =
 250     let npols = map2 (fun p n -> p,Start n) pols (0--(length pols - 1)) in
 251     let phists = filter (fun (p,_) -> p <> []) npols in
 252     let bas = grobner_interreduce [] (map monic phists) in
 253     let prs0 = allpairs (fun x y -> x,y) bas bas in
 254     let prs1 = filter (fun ((x,_),(y,_)) -> poly_lt x y) prs0 in
 255     let prs2 = map (fun (p,q) -> mlcm (hd(fst p)) (hd(fst q)),(p,q)) prs1 in
 256     let prs3 =
 257       filter (fun (l,(p,q)) -> not(orthogonal l (fst p) (fst q))) prs2 in
 258     grobner_basis bas (mergesort forder prs3) in
 259
 260   (* ----------------------------------------------------------------------- *)
 261   (* Get proof of contradiction from Grobner basis.                          *)
 262   (* ----------------------------------------------------------------------- *)
 263
 264   let grobner_refute pols =
 265     let gb = grobner pols in
 266     snd(find (fun (p,h) -> length p = 1 & forall ((=)0) (snd(hd p))) gb) in
 267
 268   (* ----------------------------------------------------------------------- *)
 269   (* Turn proof into a certificate as sum of multipliers.                    *)
 270   (*                                                                         *)
 271   (* In principle this is very inefficient: in a heavily shared proof it may *)
 272   (* make the same calculation many times. Could add a cache or something.   *)
 273   (* ----------------------------------------------------------------------- *)
 274
 275   let rec resolve_proof vars prf =
 276     match prf with
 277       Start(-1) -> []
 278     | Start m -> [m,[num_1,map (K 0) vars]]
 279     | Mmul(pol,lin) ->
 280           let lis = resolve_proof vars lin in
 281           map (fun (n,p) -> n,grob_cmul pol p) lis
 282     | Add(lin1,lin2) ->
 283           let lis1 = resolve_proof vars lin1
 284           and lis2 = resolve_proof vars lin2 in
 285           let dom = setify(union (map fst lis1) (map fst lis2)) in
 286           map (fun n -> let a = try assoc n lis1 with Failure _ -> []
 287                         and b = try assoc n lis2 with Failure _ -> [] in
 288                         n,grob_add a b) dom in
 289
 290   (* ----------------------------------------------------------------------- *)
 291   (* Run the procedure and produce Weak Nullstellensatz certificate.         *)
 292   (* ----------------------------------------------------------------------- *)
 293
 294   let grobner_weak vars pols =
 295     let cert = resolve_proof vars (grobner_refute pols) in
 296     let l =
 297       itlist (itlist (lcm_num o denominator o fst) o snd) cert (num_1) in
 298     l,map (fun (i,p) -> i,map (fun (d,m) -> (l*/d,m)) p) cert in
 299
 300   (* ----------------------------------------------------------------------- *)
 301   (* Prove polynomial is in ideal generated by others, using Grobner basis.  *)
 302   (* ----------------------------------------------------------------------- *)
 303
 304   let grobner_ideal vars pols pol =
 305     let pol',h = reduce (grobner pols) (grob_neg pol,Start(-1)) in
 306     if pol' <> [] then failwith "grobner_ideal: not in the ideal" else
 307     resolve_proof vars h in
 308
 309   (* ----------------------------------------------------------------------- *)
 310   (* Produce Strong Nullstellensatz certificate for a power of pol.          *)
 311   (* ----------------------------------------------------------------------- *)
 312
 313   let grobner_strong vars pols pol =
 314     if pol = [] then 1,num_1,[] else
 315     let vars' = (concl TRUTH)::vars in
 316     let grob_z = [num_1,1::(map (fun x -> 0) vars)]
 317     and grob_1 = [num_1,(map (fun x -> 0) vars')]
 318     and augment = map (fun (c,m) -> (c,0::m)) in
 319     let pols' = map augment pols
 320     and pol' = augment pol in
 321     let allpols = (grob_sub (grob_mul grob_z pol') grob_1)::pols' in
 322     let l,cert = grobner_weak vars' allpols in
 323     let d = itlist (itlist (max o hd o snd) o snd) cert 0 in
 324     let transform_monomial (c,m) =
 325       grob_cmul (c,tl m) (grob_pow vars pol (d - hd m)) in
 326     let transform_polynomial q = itlist (grob_add o transform_monomial) q [] in
 327     let cert' = map (fun (c,q) -> c-1,transform_polynomial q)
 328                     (filter (fun (k,_) -> k <> 0) cert) in
 329     d,l,cert' in
 330
 331   (* ----------------------------------------------------------------------- *)
 332   (* Overall parametrized universal procedure for (semi)rings.               *)
 333   (* We return an IDEAL_CONV and the actual ring prover.                     *)
 334   (* ----------------------------------------------------------------------- *)
 335
 336   let pth_step = prove
 337    (`!(add:A->A->A) (mul:A->A->A) (n0:A).
 338           (!x. mul n0 x = n0) /\
 339           (!x y z. (add x y = add x z) <=> (y = z)) /\
 340           (!w x y z. (add (mul w y) (mul x z) = add (mul w z) (mul x y)) <=>
 341                      (w = x) \/ (y = z))
 342           ==> (!a b c d. ~(a = b) /\ ~(c = d) <=>
 343                          ~(add (mul a c) (mul b d) =
 344                            add (mul a d) (mul b c))) /\
 345               (!n a b c d. ~(n = n0)
 346                            ==> (a = b) /\ ~(c = d)
 347                                ==> ~(add a (mul n c) = add b (mul n d)))`,
 348     REPEAT GEN_TAC THEN STRIP_TAC THEN
 349     ASM_REWRITE_TAC[GSYM DE_MORGAN_THM] THEN
 350     REPEAT GEN_TAC THEN DISCH_TAC THEN STRIP_TAC THEN
 351     FIRST_X_ASSUM(MP_TAC o SPECL [`n0:A`; `n:A`; `d:A`; `c:A`]) THEN
 352     ONCE_REWRITE_TAC[GSYM CONTRAPOS_THM] THEN ASM_SIMP_TAC[])
 353   and FINAL_RULE = MATCH_MP(TAUT `(p ==> F) ==> (~q = p) ==> q`)
 354   and false_tm = `F` in
 355   let rec refute_disj rfn tm =
 356     match tm with
 357       Comb(Comb(Const("\\/",_),l),r) ->
 358         DISJ_CASES (ASSUME tm) (refute_disj rfn l) (refute_disj rfn r)
 359     | _ -> rfn tm in
 360   fun (ring_dest_const,ring_mk_const,RING_EQ_CONV,
 361        ring_neg_tm,ring_add_tm,ring_sub_tm,
 362        ring_inv_tm,ring_mul_tm,ring_div_tm,ring_pow_tm,
 363        RING_INTEGRAL,RABINOWITSCH_THM,RING_NORMALIZE_CONV) ->
 364     let INITIAL_CONV =
 365       TOP_DEPTH_CONV BETA_CONV THENC
 366       PRESIMP_CONV THENC
 367       CONDS_ELIM_CONV THENC
 368       NNF_CONV THENC
 369       (if is_iff(snd(strip_forall(concl RABINOWITSCH_THM)))
 370        then GEN_REWRITE_CONV ONCE_DEPTH_CONV [RABINOWITSCH_THM]
 371        else ALL_CONV) THENC
 372       GEN_REWRITE_CONV REDEPTH_CONV
 373        [AND_FORALL_THM;
 374         LEFT_AND_FORALL_THM;
 375         RIGHT_AND_FORALL_THM;
 376         LEFT_OR_FORALL_THM;
 377         RIGHT_OR_FORALL_THM;
 378         OR_EXISTS_THM;
 379         LEFT_OR_EXISTS_THM;
 380         RIGHT_OR_EXISTS_THM;
 381         LEFT_AND_EXISTS_THM;
 382         RIGHT_AND_EXISTS_THM] in
 383     let ring_dest_neg t =
 384       let l,r = dest_comb t in
 385       if l = ring_neg_tm then r else failwith "ring_dest_neg"
 386     and ring_dest_inv t =
 387       let l,r = dest_comb t in
 388       if l = ring_inv_tm then r else failwith "ring_dest_inv"
 389     and ring_dest_add = dest_binop ring_add_tm
 390     and ring_mk_add = mk_binop ring_add_tm
 391     and ring_dest_sub = dest_binop ring_sub_tm
 392     and ring_dest_mul = dest_binop ring_mul_tm
 393     and ring_mk_mul = mk_binop ring_mul_tm
 394     and ring_dest_div = dest_binop ring_div_tm
 395     and ring_dest_pow = dest_binop ring_pow_tm
 396     and ring_mk_pow = mk_binop ring_pow_tm in
 397     let rec grobvars tm acc =
 398       if can ring_dest_const tm then acc
 399       else if can ring_dest_neg tm then grobvars (rand tm) acc
 400       else if can ring_dest_pow tm & is_numeral (rand tm)
 401            then grobvars (lhand tm) acc
 402       else if can ring_dest_add tm or can ring_dest_sub tm
 403            or can ring_dest_mul tm
 404            then grobvars (lhand tm) (grobvars (rand tm) acc)
 405       else if can ring_dest_inv tm then
 406            let gvs = grobvars (rand tm) [] in
 407            if gvs = [] then acc else tm::acc
 408       else if can ring_dest_div tm then
 409            let lvs = grobvars (lhand tm) acc
 410            and gvs = grobvars (rand tm) [] in
 411            if gvs = [] then lvs else tm::acc
 412       else tm::acc in
 413     let rec grobify_term vars tm =
 414       try if not(mem tm vars) then failwith "" else
 415           [num_1,map (fun i -> if i = tm then 1 else 0) vars]
 416       with Failure _ -> try
 417           let x = ring_dest_const tm in
 418           if x =/ num_0 then [] else [x,map (fun v -> 0) vars]
 419       with Failure _ -> try
 420           grob_neg(grobify_term vars (ring_dest_neg tm))
 421       with Failure _ -> try
 422           grob_inv(grobify_term vars (ring_dest_inv tm))
 423       with Failure _ -> try
 424           let l,r = ring_dest_add tm in
 425           grob_add (grobify_term vars l) (grobify_term vars r)
 426       with Failure _ -> try
 427           let l,r = ring_dest_sub tm in
 428           grob_sub (grobify_term vars l) (grobify_term vars r)
 429       with Failure _ -> try
 430           let l,r = ring_dest_mul tm in
 431           grob_mul (grobify_term vars l) (grobify_term vars r)
 432       with Failure _ -> try
 433           let l,r = ring_dest_div tm in
 434           grob_div (grobify_term vars l) (grobify_term vars r)
 435       with Failure _ -> try
 436           let l,r = ring_dest_pow tm in
 437           grob_pow vars (grobify_term vars l) (dest_small_numeral r)
 438       with Failure _ ->
 439             failwith "grobify_term: unknown or invalid term" in
 440     let grobify_equation vars tm =
 441       let l,r = dest_eq tm in
 442       grob_sub (grobify_term vars l) (grobify_term vars r) in
 443     let grobify_equations tm =
 444       let cjs = conjuncts tm in
 445       let rawvars =
 446         itlist (fun eq a -> grobvars (lhand eq) (grobvars (rand eq) a))
 447                cjs [] in
 448       let vars = sort (fun x y -> x < y) (setify rawvars) in
 449       vars,map (grobify_equation vars) cjs in
 450     let holify_polynomial =
 451       let holify_varpow (v,n) =
 452         if n = 1 then v else ring_mk_pow v (mk_small_numeral n) in
 453       let holify_monomial vars (c,m) =
 454         let xps = map holify_varpow
 455           (filter (fun (_,n) -> n <> 0) (zip vars m)) in
 456         end_itlist ring_mk_mul (ring_mk_const c :: xps) in
 457       let holify_polynomial vars p =
 458         if p = [] then ring_mk_const (num_0)
 459         else end_itlist ring_mk_add (map (holify_monomial vars) p) in
 460       holify_polynomial in
 461     let (pth_idom,pth_ine) = CONJ_PAIR(MATCH_MP pth_step RING_INTEGRAL) in
 462     let IDOM_RULE = CONV_RULE(REWR_CONV pth_idom) in
 463     let PROVE_NZ n = EQF_ELIM(RING_EQ_CONV
 464                 (mk_eq(ring_mk_const n,ring_mk_const(num_0)))) in
 465     let NOT_EQ_01 = PROVE_NZ (num_1)
 466     and INE_RULE n = MATCH_MP(MATCH_MP pth_ine (PROVE_NZ n))
 467     and MK_ADD th1 th2 = MK_COMB(AP_TERM ring_add_tm th1,th2) in
 468     let execute_proof vars eths prf =
 469       let x,th1 = SPEC_VAR(CONJUNCT1(CONJUNCT2 RING_INTEGRAL)) in
 470       let y,th2 = SPEC_VAR th1 in
 471       let z,th3 = SPEC_VAR th2 in
 472       let SUB_EQ_RULE = GEN_REWRITE_RULE I
 473         [SYM(INST [mk_comb(ring_neg_tm,z),x] th3)] in
 474       let initpols = map (CONV_RULE(BINOP_CONV RING_NORMALIZE_CONV) o
 475                           SUB_EQ_RULE) eths in
 476       let ADD_RULE th1 th2 =
 477          CONV_RULE (BINOP_CONV RING_NORMALIZE_CONV)
 478                    (MK_COMB(AP_TERM ring_add_tm th1,th2))
 479       and MUL_RULE vars m th =
 480          CONV_RULE (BINOP_CONV RING_NORMALIZE_CONV)
 481                    (AP_TERM (mk_comb(ring_mul_tm,holify_polynomial vars [m]))
 482                             th) in
 483       let execache = ref [] in
 484       let memoize prf x = (execache := (prf,x)::(!execache)); x in
 485       let rec assoceq a l =
 486         match l with
 487          [] -> failwith "assoceq"
 488         | (x,y)::t -> if x==a then y else assoceq a t in
 489       let rec run_proof vars prf =
 490         try assoceq prf (!execache) with Failure _ ->
 491         (match prf with
 492            Start m -> el m initpols
 493          | Add(p1,p2) ->
 494             memoize prf (ADD_RULE (run_proof vars p1) (run_proof vars p2))
 495          | Mmul(m,p2) ->
 496             memoize prf (MUL_RULE vars m (run_proof vars p2))) in
 497       let th = run_proof vars prf in
 498       execache := []; CONV_RULE RING_EQ_CONV th in
 499     let REFUTE tm =
 500       if tm = false_tm then ASSUME tm else
 501       let nths0,eths0 = partition (is_neg o concl) (CONJUNCTS(ASSUME tm)) in
 502       let nths = filter (is_eq o rand o concl) nths0
 503       and eths = filter (is_eq o concl) eths0 in
 504       if eths = [] then
 505         let th1 = end_itlist (fun th1 th2 -> IDOM_RULE(CONJ th1 th2)) nths in
 506         let th2 = CONV_RULE(RAND_CONV(BINOP_CONV RING_NORMALIZE_CONV)) th1 in
 507         let l,r = dest_eq(rand(concl th2)) in
 508         EQ_MP (EQF_INTRO th2) (REFL l)
 509       else if nths = [] & not(is_var ring_neg_tm) then
 510         let vars,pols = grobify_equations(list_mk_conj(map concl eths)) in
 511         execute_proof vars eths (grobner_refute pols)
 512       else
 513       let vars,l,cert,noteqth =
 514         if nths = [] then
 515           let vars,pols = grobify_equations(list_mk_conj(map concl eths)) in
 516           let l,cert = grobner_weak vars pols in
 517           vars,l,cert,NOT_EQ_01
 518         else
 519           let nth = end_itlist
 520            (fun th1 th2 -> IDOM_RULE(CONJ th1 th2)) nths in
 521           let vars,pol::pols =
 522            grobify_equations(list_mk_conj(rand(concl nth)::map concl eths)) in
 523           let deg,l,cert = grobner_strong vars pols pol in
 524           let th1 =
 525             CONV_RULE(RAND_CONV(BINOP_CONV RING_NORMALIZE_CONV)) nth in
 526           let th2 = funpow deg (IDOM_RULE o CONJ th1) NOT_EQ_01 in
 527           vars,l,cert,th2 in
 528       Format.print_string("Translating certificate to HOL inferences");
 529       Format.print_newline();
 530       let cert_pos = map
 531         (fun (i,p) -> i,filter (fun (c,m) -> c >/ num_0) p) cert
 532       and cert_neg = map
 533         (fun (i,p) -> i,map (fun (c,m) -> minus_num c,m)
 534                             (filter (fun (c,m) -> c </ num_0) p)) cert in
 535       let herts_pos =
 536         map (fun (i,p) -> i,holify_polynomial vars p) cert_pos
 537       and herts_neg =
 538         map (fun (i,p) -> i,holify_polynomial vars p) cert_neg in
 539       let thm_fn pols =
 540         if pols = [] then REFL(ring_mk_const num_0) else
 541         end_itlist MK_ADD
 542         (map (fun (i,p) -> AP_TERM(mk_comb(ring_mul_tm,p)) (el i eths))
 543              pols) in
 544       let th1 = thm_fn herts_pos and th2 = thm_fn herts_neg in
 545       let th3 = CONJ(MK_ADD (SYM th1) th2) noteqth in
 546       let th4 = CONV_RULE (RAND_CONV(BINOP_CONV RING_NORMALIZE_CONV))
 547                           (INE_RULE l th3) in
 548       let l,r = dest_eq(rand(concl th4)) in
 549       EQ_MP (EQF_INTRO th4) (REFL l) in
 550   let RING tm =
 551     let avs = frees tm in
 552     let tm' = list_mk_forall(avs,tm) in
 553     let th1 = INITIAL_CONV(mk_neg tm') in
 554     let evs,bod = strip_exists(rand(concl th1)) in
 555     if is_forall bod then failwith "RING: non-universal formula" else
 556     let th1a = WEAK_DNF_CONV bod in
 557     let boda = rand(concl th1a) in
 558     let th2a = refute_disj REFUTE boda in
 559     let th2b = TRANS th1a (EQF_INTRO(NOT_INTRO(DISCH boda th2a))) in
 560     let th2 = UNDISCH(NOT_ELIM(EQF_ELIM th2b)) in
 561     let th3 = itlist SIMPLE_CHOOSE evs th2 in
 562     SPECL avs (MATCH_MP (FINAL_RULE (DISCH_ALL th3)) th1)
 563   and ideal tms tm =
 564     let rawvars = itlist grobvars (tm::tms) [] in
 565     let vars = sort (fun x y -> x < y) (setify rawvars) in
 566     let pols = map (grobify_term vars) tms and pol = grobify_term vars tm in
 567     let cert = grobner_ideal vars pols pol in
 568     map (fun n -> let p = assocd n cert [] in holify_polynomial vars p)
 569         (0--(length pols-1)) in
 570   RING,ideal;;
 571
 572 (* ----------------------------------------------------------------------- *)
 573 (* Separate out the cases.                                                 *)
 574 (* ----------------------------------------------------------------------- *)
 575
 576 let RING parms = fst(RING_AND_IDEAL_CONV parms);;
 577
 578 let ideal_cofactors parms = snd(RING_AND_IDEAL_CONV parms);;
 579
 580 (* ------------------------------------------------------------------------- *)
 581 (* Simplify a natural number assertion to eliminate conditionals, DIV, MOD,  *)
 582 (* PRE, cutoff subtraction, EVEN and ODD. Try to do it in a way that makes   *)
 583 (* new quantifiers universal. At the moment we don't split "<=>" which would *)
 584 (* make this quantifier selection work there too; better to do NNF first if  *)
 585 (* you care. This also applies to EVEN and ODD.                              *)
 586 (* ------------------------------------------------------------------------- *)
 587
 588 let NUM_SIMPLIFY_CONV =
 589   let pre_tm = `PRE`
 590   and div_tm = `(DIV):num->num->num`
 591   and mod_tm = `(MOD):num->num->num`
 592   and p_tm = `P:num->bool` and n_tm = `n:num` and m_tm = `m:num`
 593   and q_tm = `P:num->num->bool` and a_tm = `a:num` and b_tm = `b:num` in
 594   let is_pre tm = is_comb tm & rator tm = pre_tm
 595   and is_sub = is_binop `(-):num->num->num`
 596   and is_divmod =
 597     let is_div = is_binop div_tm and is_mod = is_binop mod_tm in
 598     fun tm -> is_div tm or is_mod tm
 599   and contains_quantifier =
 600     can (find_term (fun t -> is_forall t or is_exists t or is_uexists t))
 601   and BETA2_CONV = RATOR_CONV BETA_CONV THENC BETA_CONV
 602   and PRE_ELIM_THM'' = CONV_RULE (RAND_CONV NNF_CONV) PRE_ELIM_THM
 603   and SUB_ELIM_THM'' = CONV_RULE (RAND_CONV NNF_CONV) SUB_ELIM_THM
 604   and DIVMOD_ELIM_THM'' = CONV_RULE (RAND_CONV NNF_CONV) DIVMOD_ELIM_THM
 605   and pth_evenodd = prove
 606    (`(EVEN(x) <=> (!y. ~(x = SUC(2 * y)))) /\
 607      (ODD(x) <=> (!y. ~(x = 2 * y))) /\
 608      (~EVEN(x) <=> (!y. ~(x = 2 * y))) /\
 609      (~ODD(x) <=> (!y. ~(x = SUC(2 * y))))`,
 610     REWRITE_TAC[GSYM NOT_EXISTS_THM; GSYM EVEN_EXISTS; GSYM ODD_EXISTS] THEN
 611     REWRITE_TAC[NOT_EVEN; NOT_ODD]) in
 612   let rec NUM_MULTIPLY_CONV pos tm =
 613     if is_forall tm or is_exists tm or is_uexists tm then
 614        BINDER_CONV (NUM_MULTIPLY_CONV pos) tm
 615     else if is_imp tm & contains_quantifier tm then
 616         COMB2_CONV (RAND_CONV(NUM_MULTIPLY_CONV(not pos)))
 617                    (NUM_MULTIPLY_CONV pos) tm
 618     else if (is_conj tm or is_disj tm or is_iff tm) &
 619             contains_quantifier tm
 620          then BINOP_CONV (NUM_MULTIPLY_CONV pos) tm
 621     else if is_neg tm & not pos & contains_quantifier tm then
 622        RAND_CONV (NUM_MULTIPLY_CONV (not pos)) tm
 623     else
 624        try let t = find_term (fun t -> is_pre t & free_in t tm) tm in
 625            let ty = type_of t in
 626            let v = genvar ty in
 627            let p = mk_abs(v,subst [v,t] tm) in
 628            let th0 = if pos then PRE_ELIM_THM'' else PRE_ELIM_THM' in
 629            let th1 = INST [p,p_tm; rand t,n_tm] th0 in
 630            let th2 = CONV_RULE(COMB2_CONV (RAND_CONV BETA_CONV)
 631                       (BINDER_CONV(RAND_CONV BETA_CONV))) th1 in
 632            CONV_RULE(RAND_CONV (NUM_MULTIPLY_CONV pos)) th2
 633        with Failure _ -> try
 634            let t = find_term (fun t -> is_sub t & free_in t tm) tm in
 635            let ty = type_of t in
 636            let v = genvar ty in
 637            let p = mk_abs(v,subst [v,t] tm) in
 638            let th0 = if pos then SUB_ELIM_THM'' else SUB_ELIM_THM' in
 639            let th1 = INST [p,p_tm; lhand t,a_tm; rand t,b_tm] th0 in
 640            let th2 = CONV_RULE(COMB2_CONV (RAND_CONV BETA_CONV)
 641                       (BINDER_CONV(RAND_CONV BETA_CONV))) th1 in
 642            CONV_RULE(RAND_CONV (NUM_MULTIPLY_CONV pos)) th2
 643        with Failure _ -> try
 644            let t = find_term (fun t -> is_divmod t & free_in t tm) tm in
 645            let x = lhand t and y = rand t in
 646            let dtm = mk_comb(mk_comb(div_tm,x),y)
 647            and mtm = mk_comb(mk_comb(mod_tm,x),y) in
 648            let vd = genvar(type_of dtm)
 649            and vm = genvar(type_of mtm) in
 650            let p = list_mk_abs([vd;vm],subst[vd,dtm; vm,mtm] tm) in
 651            let th0 = if pos then DIVMOD_ELIM_THM'' else DIVMOD_ELIM_THM' in
 652            let th1 = INST [p,q_tm; x,m_tm; y,n_tm] th0 in
 653            let th2 = CONV_RULE(COMB2_CONV(RAND_CONV BETA2_CONV)
 654                 (funpow 2 BINDER_CONV(RAND_CONV BETA2_CONV))) th1 in
 655            CONV_RULE(RAND_CONV (NUM_MULTIPLY_CONV pos)) th2
 656        with Failure _ -> REFL tm in
 657   NUM_REDUCE_CONV THENC
 658   CONDS_CELIM_CONV THENC
 659   NNF_CONV THENC
 660   NUM_MULTIPLY_CONV true THENC
 661   NUM_REDUCE_CONV THENC
 662   GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth_evenodd];;
 663
 664 (* ----------------------------------------------------------------------- *)
 665 (* Natural number version of ring procedure with this normalization.       *)
 666 (* ----------------------------------------------------------------------- *)
 667
 668 let NUM_RING =
 669   let NUM_INTEGRAL_LEMMA = prove
 670    (`(w = x + d) /\ (y = z + e)
 671      ==> ((w * y + x * z = w * z + x * y) <=> (w = x) \/ (y = z))`,
 672     DISCH_THEN(fun th -> REWRITE_TAC[th]) THEN
 673     REWRITE_TAC[LEFT_ADD_DISTRIB; RIGHT_ADD_DISTRIB; GSYM ADD_ASSOC] THEN
 674     ONCE_REWRITE_TAC[AC ADD_AC
 675      `a + b + c + d + e = a + c + e + b + d`] THEN
 676     REWRITE_TAC[EQ_ADD_LCANCEL; EQ_ADD_LCANCEL_0; MULT_EQ_0]) in
 677   let NUM_INTEGRAL = prove
 678    (`(!x. 0 * x = 0) /\
 679      (!x y z. (x + y = x + z) <=> (y = z)) /\
 680      (!w x y z. (w * y + x * z = w * z + x * y) <=> (w = x) \/ (y = z))`,
 681     REWRITE_TAC[MULT_CLAUSES; EQ_ADD_LCANCEL] THEN
 682     REPEAT GEN_TAC THEN
 683     DISJ_CASES_TAC (SPECL [`w:num`; `x:num`] LE_CASES) THEN
 684     DISJ_CASES_TAC (SPECL [`y:num`; `z:num`] LE_CASES) THEN
 685     REPEAT(FIRST_X_ASSUM
 686      (CHOOSE_THEN SUBST1_TAC o REWRITE_RULE[LE_EXISTS])) THEN
 687     ASM_MESON_TAC[NUM_INTEGRAL_LEMMA; ADD_SYM; MULT_SYM]) in
 688   let rawring =
 689     RING(dest_numeral,mk_numeral,NUM_EQ_CONV,
 690          genvar bool_ty,`(+):num->num->num`,genvar bool_ty,
 691          genvar bool_ty,`(*):num->num->num`,genvar bool_ty,
 692          `(EXP):num->num->num`,
 693          NUM_INTEGRAL,TRUTH,NUM_NORMALIZE_CONV) in
 694   let initconv = NUM_SIMPLIFY_CONV THENC GEN_REWRITE_CONV DEPTH_CONV [ADD1]
 695   and t_tm = `T` in
 696   fun tm -> let th = initconv tm in
 697             if rand(concl th) = t_tm then th
 698             else EQ_MP (SYM th) (rawring(rand(concl th)));;