1 (* ========================================================================= *)
2 (* Grobner basis algorithm. *)
3 (* ========================================================================= *)
5 needs "Complex/complexnumbers.ml";;
6 needs "Complex/quelim.ml";;
10 (* ------------------------------------------------------------------------- *)
11 (* Utility functions. *)
12 (* ------------------------------------------------------------------------- *)
14 let allpairs f l1 l2 =
15 itlist ((@) o C map l2 o f) l1 [];;
17 let rec merge ord l1 l2 =
20 | h1::t1 -> match l2 with
22 | h2::t2 -> if ord h1 h2 then h1::(merge ord t1 l2)
23 else h2::(merge ord l1 t2);;
26 let rec mergepairs l1 l2 =
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);;
34 (* ------------------------------------------------------------------------- *)
35 (* Type for recording history, i.e. how a polynomial was obtained. *)
36 (* ------------------------------------------------------------------------- *)
40 | Mmul of (num * (int list)) * history
41 | Add of history * history;;
43 (* ------------------------------------------------------------------------- *)
44 (* Conversion of leaves, i.e. variables and constant rational constants. *)
45 (* ------------------------------------------------------------------------- *)
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";;
54 try let l,r = dest_comb tm in
56 let x = rat_of_term r in
57 if x =/ Int 0 then [] else [x,map (fun v -> 0) vars]
59 with Failure _ -> failwith "grob_const";;
61 (* ------------------------------------------------------------------------- *)
62 (* Monomial ordering. *)
63 (* ------------------------------------------------------------------------- *)
66 let rec lexorder l1 l2 =
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;;
75 let morder_le m1 m2 = morder_lt m1 m2 or (m1 = m2);;
77 let morder_gt m1 m2 = morder_lt m2 m1;;
79 (* ------------------------------------------------------------------------- *)
80 (* Arithmetic on canonical polynomials. *)
81 (* ------------------------------------------------------------------------- *)
83 let grob_neg = map (fun (c,m) -> (minus_num c,m));;
85 let rec grob_add l1 l2 =
89 | ((c1,m1)::o1,(c2,m2)::o2) ->
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);;
96 let grob_sub l1 l2 = grob_add l1 (grob_neg l2);;
98 let grob_mmul (c1,m1) (c2,m2) = (c1*/c2,map2 (+) m1 m2);;
100 let rec grob_cmul cm pol = map (grob_mmul cm) pol;;
102 let rec grob_mul l1 l2 =
105 | (h1::t1) -> grob_add (grob_cmul h1 l2) (grob_mul t1 l2);;
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));;
112 (* ------------------------------------------------------------------------- *)
113 (* Monomial division operation. *)
114 (* ------------------------------------------------------------------------- *)
116 let mdiv (c1,m1) (c2,m2) =
118 map2 (fun n1 n2 -> if n1 < n2 then failwith "mdiv" else n1-n2) m1 m2);;
120 (* ------------------------------------------------------------------------- *)
121 (* Lowest common multiple of two monomials. *)
122 (* ------------------------------------------------------------------------- *)
124 let mlcm (c1,m1) (c2,m2) = (Int 1,map2 max m1 m2);;
126 (* ------------------------------------------------------------------------- *)
127 (* Reduce monomial cm by polynomial pol, returning replacement for cm. *)
128 (* ------------------------------------------------------------------------- *)
130 let reduce1 cm (pol,hpol) =
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";;
138 (* ------------------------------------------------------------------------- *)
139 (* Try this for all polynomials in a basis. *)
140 (* ------------------------------------------------------------------------- *)
142 let reduceb cm basis = tryfind (fun p -> reduce1 cm p) basis;;
144 (* ------------------------------------------------------------------------- *)
145 (* Reduction of a polynomial (always picking largest monomial possible). *)
146 (* ------------------------------------------------------------------------- *)
148 let rec reduce basis (pol,hist) =
151 | cm::ptl -> try let q,hnew = reduceb cm basis in
152 reduce basis (grob_add q ptl,Add(hnew,hist))
154 let q,hist' = reduce basis (ptl,hist) in
157 (* ------------------------------------------------------------------------- *)
158 (* Check for orthogonality w.r.t. LCM. *)
159 (* ------------------------------------------------------------------------- *)
161 let orthogonal l p1 p2 =
162 snd l = snd(grob_mmul (hd p1) (hd p2));;
164 (* ------------------------------------------------------------------------- *)
165 (* Compute S-polynomial of two polynomials. *)
166 (* ------------------------------------------------------------------------- *)
168 let spoly cm ph1 ph2 =
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)));;
178 (* ------------------------------------------------------------------------- *)
179 (* Make a polynomial monic. *)
180 (* ------------------------------------------------------------------------- *)
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));;
188 (* ------------------------------------------------------------------------- *)
189 (* The most popular heuristic is to order critical pairs by LCM monomial. *)
190 (* ------------------------------------------------------------------------- *)
192 let forder ((c1,m1),_) ((c2,m2),_) = morder_lt m1 m2;;
194 (* ------------------------------------------------------------------------- *)
195 (* Stupid stuff forced on us by lack of equality test on num type. *)
196 (* ------------------------------------------------------------------------- *)
198 let rec poly_lt p q =
202 | (c1,m1)::o1,(c2,m2)::o2 ->
204 c1 =/ c2 & (m1 < m2 or m1 = m2 & poly_lt o1 o2);;
206 let align ((p,hp),(q,hq)) =
207 if poly_lt p q then ((p,hp),(q,hq)) else ((q,hq),(p,hp));;
210 forall2 (fun (c1,m1) (c2,m2) -> c1 =/ c2 & m1 = m2) p1 p2;;
212 let memx ((p1,h1),(p2,h2)) ppairs =
213 not (exists (fun ((q1,_),(q2,_)) -> poly_eq p1 q1 & poly_eq p2 q2) ppairs);;
215 (* ------------------------------------------------------------------------- *)
216 (* Buchberger's second criterion. *)
217 (* ------------------------------------------------------------------------- *)
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;;
225 (* ------------------------------------------------------------------------- *)
226 (* Test for hitting constant polynomial. *)
227 (* ------------------------------------------------------------------------- *)
229 let constant_poly p =
230 length p = 1 & forall ((=) 0) (snd(hd p));;
232 (* ------------------------------------------------------------------------- *)
233 (* Grobner basis algorithm. *)
234 (* ------------------------------------------------------------------------- *)
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");
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
248 then grobner ((sp,hist)::histories) (sph::basis) [] else
250 map (fun p -> mlcm (hd(fst p)) (hd sp),align(p,sph)) basis in
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));;
256 (* ------------------------------------------------------------------------- *)
257 (* Overall function. *)
258 (* ------------------------------------------------------------------------- *)
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);;
272 (* ------------------------------------------------------------------------- *)
273 (* Alternative orthography. *)
274 (* ------------------------------------------------------------------------- *)
276 let gr'o'bner = groebner;;
278 (* ------------------------------------------------------------------------- *)
279 (* Conversion from HOL term. *)
280 (* ------------------------------------------------------------------------- *)
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
295 grob_pow vars (grobify_term vars l) (dest_small_numeral r)
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
301 try grobify_term vars tm with Failure _ -> failwith "grobify_term";;
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)
320 fun tm -> setify(grobvars tm []);;
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
331 let cjs = conjuncts tm in
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;;
338 (* ------------------------------------------------------------------------- *)
340 (* ------------------------------------------------------------------------- *)
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
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;;
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);;
354 (* ------------------------------------------------------------------------- *)
355 (* Resolve a proof. *)
356 (* ------------------------------------------------------------------------- *)
358 let rec resolve_proof vars prf =
361 [n,[Int 1,map (K 0) vars]]
363 let lis = resolve_proof vars lin in
364 map (fun (n,p) -> n,grob_cmul pol p) lis
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;;
373 (* ------------------------------------------------------------------------- *)
374 (* Convert a polynomial back to HOL. *)
375 (* ------------------------------------------------------------------------- *)
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) =
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
397 (* ------------------------------------------------------------------------- *)
398 (* Recursively find the set of basis elements involved. *)
399 (* ------------------------------------------------------------------------- *)
402 let rec dependencies prf acc =
405 | Mmul(pol,lin) -> dependencies lin acc
406 | Add(lin1,lin2) -> dependencies lin1 (dependencies lin2 acc) in
407 fun prf -> setify(dependencies prf []);;
409 let rec involved deps sofar todo =
413 if mem h sofar then involved deps sofar hs
414 else involved deps (h::sofar) (el h deps @ hs);;
416 (* ------------------------------------------------------------------------- *)
417 (* Refute a conjunction of equations in HOL. *)
418 (* ------------------------------------------------------------------------- *)
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 =
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
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
438 let vars,pols = grobify_equations tm in
439 let (prfs,gb) = groebner pols in
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
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
451 (holify_lincombs vars ths (snd(chop_list(length ths) (pprfs @ [ppr])))) in
452 CONV_RULE COMPLEX_RAT_EQ_CONV th;;
454 (* ------------------------------------------------------------------------- *)
455 (* Overall conversion. *)
456 (* ------------------------------------------------------------------------- *)
459 let pth0 = SIMPLE_COMPLEX_ARITH `(x = y) <=> (x - y = Cx(&0))`
461 (`!x. ~(x = Cx(&0)) <=> ?z. z * x + Cx(&1) = Cx(&0)`,
462 CONV_TAC(time FULL_COMPLEX_QUELIM_CONV))
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))
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
475 RIGHT_AND_FORALL_THM;
482 RIGHT_AND_EXISTS_THM] in
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
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);;