author  chaieb 
Mon, 08 Oct 2007 20:20:13 +0200  
changeset 24913  eb6fd8f78d56 
parent 24630  351a308ab58d 
child 25251  759bffe1d416 
permissions  rwrr 
23252  1 
(* Title: HOL/Tools/Groebner_Basis/groebner.ML 
2 
ID: $Id$ 

3 
Author: Amine Chaieb, TU Muenchen 

4 
*) 

5 

6 
signature GROEBNER = 

7 
sig 

8 
val ring_and_ideal_conv : 

9 
{idom: thm list, ring: cterm list * thm list, vars: cterm list, 

10 
semiring: Thm.cterm list * thm list} > 

23487  11 
(cterm > Rat.rat) > (Rat.rat > cterm) > 
12 
conv > conv > 

13 
conv * (cterm list > cterm > (cterm * cterm > order) > cterm list) 

23579  14 
val ring_tac: thm list > thm list > Proof.context > int > tactic 
23252  15 
end 
16 

17 
structure Groebner: GROEBNER = 

18 
struct 

19 

23557  20 
open Conv Misc Normalizer; 
21 

23252  22 
val rat_0 = Rat.zero; 
23 
val rat_1 = Rat.one; 

24 
val minus_rat = Rat.neg; 

25 
val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int; 

26 
fun int_of_rat a = 

27 
case Rat.quotient_of_rat a of (i,1) => i  _ => error "int_of_rat: not an int"; 

23514  28 
val lcm_rat = fn x => fn y => Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y)); 
23252  29 

30 
val (eqF_intr, eqF_elim) = 

31 
let val [th1,th2] = thms "PFalse" 

32 
in (fn th => th COMP th2, fn th => th COMP th1) end; 

33 

34 
val (PFalse, PFalse') = 

35 
let val PFalse_eq = nth (thms "simp_thms") 13 

36 
in (PFalse_eq RS iffD1, PFalse_eq RS iffD2) end; 

37 

38 

39 
(*  *) 

40 
(* Type for recording history, i.e. how a polynomial was obtained. *) 

41 
(*  *) 

42 

43 
datatype history = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

44 
Start of int 
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

45 
 Mmul of (Rat.rat * int list) * history 
23252  46 
 Add of history * history; 
47 

48 

49 
(* Monomial ordering. *) 

50 

51 
fun morder_lt m1 m2= 

52 
let fun lexorder l1 l2 = 

53 
case (l1,l2) of 

54 
([],[]) => false 

55 
 (x1::o1,x2::o2) => x1 > x2 orelse x1 = x2 andalso lexorder o1 o2 

56 
 _ => error "morder: inconsistent monomial lengths" 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

57 
val n1 = Integer.sum m1 
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

58 
val n2 = Integer.sum m2 in 
23252  59 
n1 < n2 orelse n1 = n2 andalso lexorder m1 m2 
60 
end; 

61 

62 
fun morder_le m1 m2 = morder_lt m1 m2 orelse (m1 = m2); 

63 

64 
fun morder_gt m1 m2 = morder_lt m2 m1; 

65 

66 
(* Arithmetic on canonical polynomials. *) 

67 

68 
fun grob_neg l = map (fn (c,m) => (minus_rat c,m)) l; 

69 

70 
fun grob_add l1 l2 = 

71 
case (l1,l2) of 

72 
([],l2) => l2 

73 
 (l1,[]) => l1 

74 
 ((c1,m1)::o1,(c2,m2)::o2) => 

75 
if m1 = m2 then 

76 
let val c = c1+/c2 val rest = grob_add o1 o2 in 

77 
if c =/ rat_0 then rest else (c,m1)::rest end 

78 
else if morder_lt m2 m1 then (c1,m1)::(grob_add o1 l2) 

79 
else (c2,m2)::(grob_add l1 o2); 

80 

81 
fun grob_sub l1 l2 = grob_add l1 (grob_neg l2); 

82 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

83 
fun grob_mmul (c1,m1) (c2,m2) = (c1*/c2, ListPair.map (op +) (m1, m2)); 
23252  84 

85 
fun grob_cmul cm pol = map (grob_mmul cm) pol; 

86 

87 
fun grob_mul l1 l2 = 

88 
case l1 of 

89 
[] => [] 

90 
 (h1::t1) => grob_add (grob_cmul h1 l2) (grob_mul t1 l2); 

91 

92 
fun grob_inv l = 

93 
case l of 

94 
[(c,vs)] => if (forall (fn x => x = 0) vs) then 

95 
if (c =/ rat_0) then error "grob_inv: division by zero" 

96 
else [(rat_1 // c,vs)] 

97 
else error "grob_inv: nonconstant divisor polynomial" 

98 
 _ => error "grob_inv: nonconstant divisor polynomial"; 

99 

100 
fun grob_div l1 l2 = 

101 
case l2 of 

102 
[(c,l)] => if (forall (fn x => x = 0) l) then 

103 
if c =/ rat_0 then error "grob_div: division by zero" 

104 
else grob_cmul (rat_1 // c,l) l1 

105 
else error "grob_div: nonconstant divisor polynomial" 

106 
 _ => error "grob_div: nonconstant divisor polynomial"; 

107 

108 
fun grob_pow vars l n = 

109 
if n < 0 then error "grob_pow: negative power" 

110 
else if n = 0 then [(rat_1,map (fn v => 0) vars)] 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

111 
else grob_mul l (grob_pow vars l (n  1)); 
23252  112 

113 
fun degree vn p = 

114 
case p of 

115 
[] => error "Zero polynomial" 

116 
 [(c,ns)] => nth ns vn 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

117 
 (c,ns)::p' => Int.max (nth ns vn, degree vn p'); 
23252  118 

119 
fun head_deg vn p = let val d = degree vn p in 

120 
(d,fold (fn (c,r) => fn q => grob_add q [(c, map_index (fn (i,n) => if i = vn then 0 else n) r)]) (filter (fn (c,ns) => c <>/ rat_0 andalso nth ns vn = d) p) []) end; 

121 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

122 
val is_zerop = forall (fn (c,ns) => c =/ rat_0 andalso forall (curry (op =) 0) ns); 
23252  123 
val grob_pdiv = 
124 
let fun pdiv_aux vn (n,a) p k s = 

125 
if is_zerop s then (k,s) else 

126 
let val (m,b) = head_deg vn s 

127 
in if m < n then (k,s) else 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

128 
let val p' = grob_mul p [(rat_1, map_index (fn (i,v) => if i = vn then m  n else 0) 
23252  129 
(snd (hd s)))] 
130 
in if a = b then pdiv_aux vn (n,a) p k (grob_sub s p') 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

131 
else pdiv_aux vn (n,a) p (k + 1) (grob_sub (grob_mul a s) (grob_mul b p')) 
23252  132 
end 
133 
end 

134 
in fn vn => fn s => fn p => pdiv_aux vn (head_deg vn p) p 0 s 

135 
end; 

136 

137 
(* Monomial division operation. *) 

138 

139 
fun mdiv (c1,m1) (c2,m2) = 

140 
(c1//c2, 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

141 
map2 (fn n1 => fn n2 => if n1 < n2 then error "mdiv" else n1  n2) m1 m2); 
23252  142 

143 
(* Lowest common multiple of two monomials. *) 

144 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

145 
fun mlcm (c1,m1) (c2,m2) = (rat_1, ListPair.map Int.max (m1, m2)); 
23252  146 

147 
(* Reduce monomial cm by polynomial pol, returning replacement for cm. *) 

148 

149 
fun reduce1 cm (pol,hpol) = 

150 
case pol of 

151 
[] => error "reduce1" 

152 
 cm1::cms => ((let val (c,m) = mdiv cm cm1 in 

153 
(grob_cmul (minus_rat c,m) cms, 

154 
Mmul((minus_rat c,m),hpol)) end) 

155 
handle ERROR _ => error "reduce1"); 

156 

157 
(* Try this for all polynomials in a basis. *) 

158 
fun tryfind f l = 

159 
case l of 

160 
[] => error "tryfind" 

161 
 (h::t) => ((f h) handle ERROR _ => tryfind f t); 

162 

163 
fun reduceb cm basis = tryfind (fn p => reduce1 cm p) basis; 

164 

165 
(* Reduction of a polynomial (always picking largest monomial possible). *) 

166 

167 
fun reduce basis (pol,hist) = 

168 
case pol of 

169 
[] => (pol,hist) 

170 
 cm::ptl => ((let val (q,hnew) = reduceb cm basis in 

171 
reduce basis (grob_add q ptl,Add(hnew,hist)) end) 

172 
handle (ERROR _) => 

173 
(let val (q,hist') = reduce basis (ptl,hist) in 

174 
(cm::q,hist') end)); 

175 

176 
(* Check for orthogonality w.r.t. LCM. *) 

177 

178 
fun orthogonal l p1 p2 = 

179 
snd l = snd(grob_mmul (hd p1) (hd p2)); 

180 

181 
(* Compute Spolynomial of two polynomials. *) 

182 

183 
fun spoly cm ph1 ph2 = 

184 
case (ph1,ph2) of 

185 
(([],h),p) => ([],h) 

186 
 (p,([],h)) => ([],h) 

187 
 ((cm1::ptl1,his1),(cm2::ptl2,his2)) => 

188 
(grob_sub (grob_cmul (mdiv cm cm1) ptl1) 

189 
(grob_cmul (mdiv cm cm2) ptl2), 

190 
Add(Mmul(mdiv cm cm1,his1), 

191 
Mmul(mdiv (minus_rat(fst cm),snd cm) cm2,his2))); 

192 

193 
(* Make a polynomial monic. *) 

194 

195 
fun monic (pol,hist) = 

23579  196 
if null pol then (pol,hist) else 
23252  197 
let val (c',m') = hd pol in 
198 
(map (fn (c,m) => (c//c',m)) pol, 

199 
Mmul((rat_1 // c',map (K 0) m'),hist)) end; 

200 

201 
(* The most popular heuristic is to order critical pairs by LCM monomial. *) 

202 

203 
fun forder ((c1,m1),_) ((c2,m2),_) = morder_lt m1 m2; 

204 

205 
fun poly_lt p q = 

206 
case (p,q) of 

207 
(p,[]) => false 

208 
 ([],q) => true 

209 
 ((c1,m1)::o1,(c2,m2)::o2) => 

210 
c1 </ c2 orelse 

211 
c1 =/ c2 andalso ((morder_lt m1 m2) orelse m1 = m2 andalso poly_lt o1 o2); 

212 

213 
fun align ((p,hp),(q,hq)) = 

214 
if poly_lt p q then ((p,hp),(q,hq)) else ((q,hq),(p,hp)); 

215 
fun forall2 p l1 l2 = 

216 
case (l1,l2) of 

217 
([],[]) => true 

218 
 (h1::t1,h2::t2) => p h1 h2 andalso forall2 p t1 t2 

219 
 _ => false; 

220 

221 
fun poly_eq p1 p2 = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

222 
forall2 (fn (c1,m1) => fn (c2,m2) => c1 =/ c2 andalso (m1: int list) = m2) p1 p2; 
23252  223 

224 
fun memx ((p1,h1),(p2,h2)) ppairs = 

225 
not (exists (fn ((q1,_),(q2,_)) => poly_eq p1 q1 andalso poly_eq p2 q2) ppairs); 

226 

227 
(* Buchberger's second criterion. *) 

228 

229 
fun criterion2 basis (lcm,((p1,h1),(p2,h2))) opairs = 

230 
exists (fn g => not(poly_eq (fst g) p1) andalso not(poly_eq (fst g) p2) andalso 

231 
can (mdiv lcm) (hd(fst g)) andalso 

232 
not(memx (align (g,(p1,h1))) (map snd opairs)) andalso 

233 
not(memx (align (g,(p2,h2))) (map snd opairs))) basis; 

234 

235 
(* Test for hitting constant polynomial. *) 

236 

237 
fun constant_poly p = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

238 
length p = 1 andalso forall (fn x => x = 0) (snd(hd p)); 
23252  239 

240 
(*  *) 

241 
(* Grobner basis algorithm. *) 

242 
(*  *) 

243 
(* FIXME: try to get rid of mergesort? *) 

244 
fun merge ord l1 l2 = 

245 
case l1 of 

246 
[] => l2 

247 
 h1::t1 => 

248 
case l2 of 

249 
[] => l1 

250 
 h2::t2 => if ord h1 h2 then h1::(merge ord t1 l2) 

251 
else h2::(merge ord l1 t2); 

252 
fun mergesort ord l = 

253 
let 

254 
fun mergepairs l1 l2 = 

255 
case (l1,l2) of 

256 
([s],[]) => s 

257 
 (l,[]) => mergepairs [] l 

258 
 (l,[s1]) => mergepairs (s1::l) [] 

259 
 (l,(s1::s2::ss)) => mergepairs ((merge ord s1 s2)::l) ss 

23579  260 
in if null l then [] else mergepairs [] (map (fn x => [x]) l) 
23252  261 
end; 
262 

263 

264 
fun grobner_basis basis pairs = 

23334  265 
((* writeln (Int.toString(length basis)^" basis elements and "^ 
266 
Int.toString(length pairs)^" critical pairs"); *) 

23252  267 
case pairs of 
268 
[] => basis 

269 
 (l,(p1,p2))::opairs => 

270 
let val (sph as (sp,hist)) = monic (reduce basis (spoly l p1 p2)) 

23579  271 
in if null sp orelse criterion2 basis (l,(p1,p2)) opairs 
23252  272 
then grobner_basis basis opairs 
273 
else if constant_poly sp then grobner_basis (sph::basis) [] 

274 
else let val rawcps = map (fn p => (mlcm (hd(fst p)) (hd sp),align(p,sph))) 

275 
basis 

276 
val newcps = filter 

277 
(fn (l,(p,q)) => not(orthogonal l (fst p) (fst q))) 

278 
rawcps 

279 
in grobner_basis (sph::basis) 

280 
(merge forder opairs (mergesort forder newcps)) 

281 
end 

282 
end); 

283 

284 
(*  *) 

285 
(* Interreduce initial polynomials. *) 

286 
(*  *) 

287 

288 
fun grobner_interreduce rpols ipols = 

289 
case ipols of 

290 
[] => map monic (rev rpols) 

291 
 p::ps => let val p' = reduce (rpols @ ps) p in 

23579  292 
if null (fst p') then grobner_interreduce rpols ps 
23252  293 
else grobner_interreduce (p'::rpols) ps end; 
294 

295 
(*  *) 

296 
(* Overall function. *) 

297 
(*  *) 

298 

299 
fun grobner pols = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

300 
let val npols = map2 (fn p => fn n => (p,Start n)) pols (0 upto (length pols  1)) 
23579  301 
val phists = filter (fn (p,_) => not (null p)) npols 
23252  302 
val bas = grobner_interreduce [] (map monic phists) 
303 
val prs0 = product bas bas 

304 
val prs1 = filter (fn ((x,_),(y,_)) => poly_lt x y) prs0 

305 
val prs2 = map (fn (p,q) => (mlcm (hd(fst p)) (hd(fst q)),(p,q))) prs1 

306 
val prs3 = 

307 
filter (fn (l,(p,q)) => not(orthogonal l (fst p) (fst q))) prs2 in 

308 
grobner_basis bas (mergesort forder prs3) end; 

309 

310 
(*  *) 

311 
(* Get proof of contradiction from Grobner basis. *) 

312 
(*  *) 

313 
fun find p l = 

314 
case l of 

315 
[] => error "find" 

316 
 (h::t) => if p(h) then h else find p t; 

317 

318 
fun grobner_refute pols = 

319 
let val gb = grobner pols in 

320 
snd(find (fn (p,h) => length p = 1 andalso forall (fn x=> x=0) (snd(hd p))) gb) 

321 
end; 

322 

323 
(*  *) 

324 
(* Turn proof into a certificate as sum of multipliers. *) 

325 
(* *) 

326 
(* In principle this is very inefficient: in a heavily shared proof it may *) 

327 
(* make the same calculation many times. Could put in a cache or something. *) 

328 
(*  *) 

329 
fun resolve_proof vars prf = 

330 
case prf of 

331 
Start(~1) => [] 

332 
 Start m => [(m,[(rat_1,map (K 0) vars)])] 

333 
 Mmul(pol,lin) => 

334 
let val lis = resolve_proof vars lin in 

335 
map (fn (n,p) => (n,grob_cmul pol p)) lis end 

336 
 Add(lin1,lin2) => 

337 
let val lis1 = resolve_proof vars lin1 

338 
val lis2 = resolve_proof vars lin2 

339 
val dom = distinct (op =) ((map fst lis1) union (map fst lis2)) 

340 
in 

23557  341 
map (fn n => let val a = these (AList.lookup (op =) lis1 n) 
342 
val b = these (AList.lookup (op =) lis2 n) 

343 
in (n,grob_add a b) end) dom end; 

23252  344 

345 
(*  *) 

346 
(* Run the procedure and produce Weak Nullstellensatz certificate. *) 

347 
(*  *) 

348 
fun grobner_weak vars pols = 

349 
let val cert = resolve_proof vars (grobner_refute pols) 

350 
val l = 

351 
fold_rev (fold_rev (lcm_rat o denominator_rat o fst) o snd) cert (rat_1) in 

352 
(l,map (fn (i,p) => (i,map (fn (d,m) => (l*/d,m)) p)) cert) end; 

353 

354 
(*  *) 

355 
(* Prove a polynomial is in ideal generated by others, using Grobner basis. *) 

356 
(*  *) 

357 

358 
fun grobner_ideal vars pols pol = 

359 
let val (pol',h) = reduce (grobner pols) (grob_neg pol,Start(~1)) in 

24913  360 
if not (null pol') then error "grobner_ideal: not in the ideal" else 
23252  361 
resolve_proof vars h end; 
362 

363 
(*  *) 

364 
(* Produce Strong Nullstellensatz certificate for a power of pol. *) 

365 
(*  *) 

366 

367 
fun grobner_strong vars pols pol = 

368 
let val vars' = @{cterm "True"}::vars 

369 
val grob_z = [(rat_1,1::(map (fn x => 0) vars))] 

370 
val grob_1 = [(rat_1,(map (fn x => 0) vars'))] 

371 
fun augment p= map (fn (c,m) => (c,0::m)) p 

372 
val pols' = map augment pols 

373 
val pol' = augment pol 

374 
val allpols = (grob_sub (grob_mul grob_z pol') grob_1)::pols' 

375 
val (l,cert) = grobner_weak vars' allpols 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

376 
val d = fold_rev (fold_rev (curry Int.max o hd o snd) o snd) cert 0 
23252  377 
fun transform_monomial (c,m) = 
378 
grob_cmul (c,tl m) (grob_pow vars pol (d  hd m)) 

379 
fun transform_polynomial q = fold_rev (grob_add o transform_monomial) q [] 

380 
val cert' = map (fn (c,q) => (c1,transform_polynomial q)) 

381 
(filter (fn (k,_) => k <> 0) cert) in 

382 
(d,l,cert') end; 

383 

384 

385 
(*  *) 

386 
(* Overall parametrized universal procedure for (semi)rings. *) 

387 
(* We return an ideal_conv and the actual ring prover. *) 

388 
(*  *) 

389 
fun refute_disj rfn tm = 

390 
case term_of tm of 

391 
Const("op ",_)$l$r => 

392 
Drule.compose_single(refute_disj rfn (Thm.dest_arg tm),2,Drule.compose_single(refute_disj rfn (Thm.dest_arg1 tm),2,disjE)) 

393 
 _ => rfn tm ; 

394 

395 
val notnotD = @{thm "notnotD"}; 

396 
fun mk_binop ct x y = 

397 
Thm.capply (Thm.capply ct x) y 

398 

399 
val mk_comb = Thm.capply; 

400 
fun is_neg t = 

401 
case term_of t of 

402 
(Const("Not",_)$p) => true 

403 
 _ => false; 

404 
fun is_eq t = 

405 
case term_of t of 

406 
(Const("op =",_)$_$_) => true 

407 
 _ => false; 

408 

409 
fun end_itlist f l = 

410 
case l of 

411 
[] => error "end_itlist" 

412 
 [x] => x 

413 
 (h::t) => f h (end_itlist f t); 

414 

415 
val list_mk_binop = fn b => end_itlist (mk_binop b); 

416 

417 
val list_dest_binop = fn b => 

418 
let fun h acc t = 

419 
((let val (l,r) = dest_binop b t in h (h acc r) l end) 

420 
handle CTERM _ => (t::acc)) (* Why had I handle _ => ? *) 

421 
in h [] 

422 
end; 

423 

424 
val strip_exists = 

425 
let fun h (acc, t) = 

426 
case (term_of t) of 

427 
Const("Ex",_)$Abs(x,T,p) => h (Thm.dest_abs NONE (Thm.dest_arg t) >> (fn v => v::acc)) 

428 
 _ => (acc,t) 

429 
in fn t => h ([],t) 

430 
end; 

431 

432 
fun is_forall t = 

433 
case term_of t of 

434 
(Const("All",_)$Abs(_,_,_)) => true 

435 
 _ => false; 

436 

437 
val mk_object_eq = fn th => th COMP meta_eq_to_obj_eq; 

438 
val bool_simps = @{thms "bool_simps"}; 

439 
val nnf_simps = @{thms "nnf_simps"}; 

440 
val nnf_conv = Simplifier.rewrite (HOL_basic_ss addsimps bool_simps addsimps nnf_simps) 

23557  441 
val weak_dnf_conv = Simplifier.rewrite (HOL_basic_ss addsimps @{thms "weak_dnf_simps"}); 
23252  442 
val initial_conv = 
443 
Simplifier.rewrite 

444 
(HOL_basic_ss addsimps nnf_simps 

445 
addsimps [not_all, not_ex] addsimps map (fn th => th RS sym) (ex_simps @ all_simps)); 

446 

447 
val specl = fold_rev (fn x => fn th => instantiate' [] [SOME x] (th RS spec)); 

448 

449 
val cTrp = @{cterm "Trueprop"}; 

450 
val cConj = @{cterm "op &"}; 

451 
val (cNot,false_tm) = (@{cterm "Not"}, @{cterm "False"}); 

23557  452 
val assume_Trueprop = mk_comb cTrp #> assume; 
23252  453 
val list_mk_conj = list_mk_binop cConj; 
454 
val conjs = list_dest_binop cConj; 

455 
val mk_neg = mk_comb cNot; 

456 

457 

458 

459 
(** main **) 

460 

461 
fun ring_and_ideal_conv 

462 
{vars, semiring = (sr_ops, sr_rules), ring = (r_ops, r_rules), idom} 

463 
dest_const mk_const ring_eq_conv ring_normalize_conv = 

464 
let 

465 
val [add_pat, mul_pat, pow_pat, zero_tm, one_tm] = sr_ops; 

466 
val [ring_add_tm, ring_mul_tm, ring_pow_tm] = 

467 
map (Thm.dest_fun o Thm.dest_fun) [add_pat, mul_pat, pow_pat]; 

468 

469 
val (ring_sub_tm, ring_neg_tm) = 

470 
(case r_ops of 

471 
[] => (@{cterm "True"}, @{cterm "True"}) 

472 
 [sub_pat, neg_pat] => (Thm.dest_fun (Thm.dest_fun sub_pat), Thm.dest_fun neg_pat)); 

473 

474 
val [idom_thm, neq_thm] = idom; 

475 

476 
val ring_dest_neg = 

477 
fn t => let val (l,r) = Thm.dest_comb t in 

478 
if could_unify(term_of l,term_of ring_neg_tm) then r else raise CTERM ("ring_dest_neg", [t]) 

479 
end 

480 

481 
val ring_mk_neg = fn tm => mk_comb (ring_neg_tm) (tm); 

482 
(* 

483 
fun ring_dest_inv t = 

484 
let val (l,r) = Thm.dest_comb t in 

485 
if could_unify(term_of l, term_of ring_inv_tm) then r else raise CTERM "ring_dest_inv" 

486 
end 

487 
*) 

488 
val ring_dest_add = dest_binop ring_add_tm; 

489 
val ring_mk_add = mk_binop ring_add_tm; 

490 
val ring_dest_sub = dest_binop ring_sub_tm; 

491 
val ring_mk_sub = mk_binop ring_sub_tm; 

492 
val ring_dest_mul = dest_binop ring_mul_tm; 

493 
val ring_mk_mul = mk_binop ring_mul_tm; 

494 
(* val ring_dest_div = dest_binop ring_div_tm; 

495 
val ring_mk_div = mk_binop ring_div_tm;*) 

496 
val ring_dest_pow = dest_binop ring_pow_tm; 

497 
val ring_mk_pow = mk_binop ring_pow_tm ; 

498 
fun grobvars tm acc = 

499 
if can dest_const tm then acc 

500 
else if can ring_dest_neg tm then grobvars (Thm.dest_arg tm) acc 

501 
else if can ring_dest_pow tm then grobvars (Thm.dest_arg1 tm) acc 

502 
else if can ring_dest_add tm orelse can ring_dest_sub tm 

503 
orelse can ring_dest_mul tm 

504 
then grobvars (Thm.dest_arg1 tm) (grobvars (Thm.dest_arg tm) acc) 

505 
(* else if can ring_dest_inv tm 

506 
then 

507 
let val gvs = grobvars (Thm.dest_arg tm) [] in 

508 
if gvs = [] then acc else tm::acc 

509 
end 

510 
else if can ring_dest_div tm then 

511 
let val lvs = grobvars (Thm.dest_arg1 tm) acc 

512 
val gvs = grobvars (Thm.dest_arg tm) [] 

513 
in if gvs = [] then lvs else tm::acc 

514 
end *) 

515 
else tm::acc ; 

516 

517 
fun grobify_term vars tm = 

518 
((if not (member (op aconvc) vars tm) then raise CTERM ("Not a variable", [tm]) else 

519 
[(rat_1,map (fn i => if i aconvc tm then 1 else 0) vars)]) 

520 
handle CTERM _ => 

521 
((let val x = dest_const tm 

522 
in if x =/ rat_0 then [] else [(x,map (fn v => 0) vars)] 

523 
end) 

524 
handle ERROR _ => 

525 
((grob_neg(grobify_term vars (ring_dest_neg tm))) 

526 
handle CTERM _ => 

527 
( 

528 
(* (grob_inv(grobify_term vars (ring_dest_inv tm))) 

529 
handle CTERM _ => *) 

530 
((let val (l,r) = ring_dest_add tm 

531 
in grob_add (grobify_term vars l) (grobify_term vars r) 

532 
end) 

533 
handle CTERM _ => 

534 
((let val (l,r) = ring_dest_sub tm 

535 
in grob_sub (grobify_term vars l) (grobify_term vars r) 

536 
end) 

537 
handle CTERM _ => 

538 
((let val (l,r) = ring_dest_mul tm 

539 
in grob_mul (grobify_term vars l) (grobify_term vars r) 

540 
end) 

541 
handle CTERM _ => 

542 
( 

543 
(* (let val (l,r) = ring_dest_div tm 

544 
in grob_div (grobify_term vars l) (grobify_term vars r) 

545 
end) 

546 
handle CTERM _ => *) 

547 
((let val (l,r) = ring_dest_pow tm 

548 
in grob_pow vars (grobify_term vars l) ((term_of #> HOLogic.dest_number #> snd) r) 

549 
end) 

550 
handle CTERM _ => error "grobify_term: unknown or invalid term"))))))))); 

551 
val eq_tm = idom_thm > concl > Thm.dest_arg > Thm.dest_arg > Thm.dest_fun > Thm.dest_fun ; 

552 
(*ring_integral > hd > concl > Thm.dest_arg 

553 
> Thm.dest_abs NONE > snd > Thm.dest_fun > Thm.dest_fun; *) 

554 
val dest_eq = dest_binop eq_tm; 

555 

556 
fun grobify_equation vars tm = 

557 
let val (l,r) = dest_binop eq_tm tm 

558 
in grob_sub (grobify_term vars l) (grobify_term vars r) 

559 
end; 

560 

561 
fun grobify_equations tm = 

562 
let 

563 
val cjs = conjs tm 

564 
val rawvars = fold_rev (fn eq => fn a => 

565 
grobvars (Thm.dest_arg1 eq) (grobvars (Thm.dest_arg eq) a)) cjs [] 

566 
val vars = sort (fn (x, y) => Term.term_ord(term_of x,term_of y)) 

567 
(distinct (op aconvc) rawvars) 

568 
in (vars,map (grobify_equation vars) cjs) 

569 
end; 

570 

571 
val holify_polynomial = 

572 
let fun holify_varpow (v,n) = 

23579  573 
if n = 1 then v else ring_mk_pow v (Numeral.mk_cnumber @{ctyp "nat"} n) (* FIXME *) 
23252  574 
fun holify_monomial vars (c,m) = 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

575 
let val xps = map holify_varpow (filter (fn (_,n) => n <> 0) (vars ~~ m)) 
23252  576 
in end_itlist ring_mk_mul (mk_const c :: xps) 
577 
end 

578 
fun holify_polynomial vars p = 

23579  579 
if null p then mk_const (rat_0) 
23252  580 
else end_itlist ring_mk_add (map (holify_monomial vars) p) 
581 
in holify_polynomial 

582 
end ; 

583 
val idom_rule = simplify (HOL_basic_ss addsimps [idom_thm]); 

584 
fun prove_nz n = eqF_elim 

585 
(ring_eq_conv(mk_binop eq_tm (mk_const n) (mk_const(rat_0)))); 

586 
val neq_01 = prove_nz (rat_1); 

587 
fun neq_rule n th = [prove_nz n, th] MRS neq_thm; 

588 
fun mk_add th1 = combination(Drule.arg_cong_rule ring_add_tm th1); 

589 

590 
fun refute tm = 

23557  591 
if tm aconvc false_tm then assume_Trueprop tm else 
23252  592 
let 
23557  593 
val (nths0,eths0) = List.partition (is_neg o concl) (HOLogic.conj_elims (assume_Trueprop tm)) 
23252  594 
val nths = filter (is_eq o Thm.dest_arg o concl) nths0 
595 
val eths = filter (is_eq o concl) eths0 

596 
in 

597 
if null eths then 

598 
let 

23557  599 
val th1 = end_itlist (fn th1 => fn th2 => idom_rule(HOLogic.conj_intr th1 th2)) nths 
23252  600 
val th2 = Conv.fconv_rule 
601 
((arg_conv #> arg_conv) 

602 
(binop_conv ring_normalize_conv)) th1 

603 
val conc = th2 > concl > Thm.dest_arg 

604 
val (l,r) = conc > dest_eq 

605 
in implies_intr (mk_comb cTrp tm) 

606 
(equal_elim (Drule.arg_cong_rule cTrp (eqF_intr th2)) 

607 
(reflexive l > mk_object_eq)) 

608 
end 

609 
else 

610 
let 

611 
val (vars,l,cert,noteqth) =( 

612 
if null nths then 

613 
let val (vars,pols) = grobify_equations(list_mk_conj(map concl eths)) 

614 
val (l,cert) = grobner_weak vars pols 

615 
in (vars,l,cert,neq_01) 

616 
end 

617 
else 

618 
let 

23557  619 
val nth = end_itlist (fn th1 => fn th2 => idom_rule(HOLogic.conj_intr th1 th2)) nths 
23252  620 
val (vars,pol::pols) = 
621 
grobify_equations(list_mk_conj(Thm.dest_arg(concl nth)::map concl eths)) 

622 
val (deg,l,cert) = grobner_strong vars pols pol 

623 
val th1 = Conv.fconv_rule((arg_conv o arg_conv)(binop_conv ring_normalize_conv)) nth 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

624 
val th2 = funpow deg (idom_rule o HOLogic.conj_intr th1) neq_01 
23252  625 
in (vars,l,cert,th2) 
626 
end) 

23334  627 
(* val _ = writeln ("Translating certificate to HOL inferences") *) 
23252  628 
val cert_pos = map (fn (i,p) => (i,filter (fn (c,m) => c >/ rat_0) p)) cert 
629 
val cert_neg = map (fn (i,p) => (i,map (fn (c,m) => (minus_rat c,m)) 

630 
(filter (fn (c,m) => c </ rat_0) p))) cert 

631 
val herts_pos = map (fn (i,p) => (i,holify_polynomial vars p)) cert_pos 

632 
val herts_neg = map (fn (i,p) => (i,holify_polynomial vars p)) cert_neg 

633 
fun thm_fn pols = 

634 
if null pols then reflexive(mk_const rat_0) else 

635 
end_itlist mk_add 

23259  636 
(map (fn (i,p) => Drule.arg_cong_rule (mk_comb ring_mul_tm p) 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

637 
(nth eths i > mk_meta_eq)) pols) 
23252  638 
val th1 = thm_fn herts_pos 
639 
val th2 = thm_fn herts_neg 

23557  640 
val th3 = HOLogic.conj_intr(mk_add (symmetric th1) th2 > mk_object_eq) noteqth 
23252  641 
val th4 = Conv.fconv_rule ((arg_conv o arg_conv o binop_conv) ring_normalize_conv) 
642 
(neq_rule l th3) 

643 
val (l,r) = dest_eq(Thm.dest_arg(concl th4)) 

644 
in implies_intr (mk_comb cTrp tm) 

645 
(equal_elim (Drule.arg_cong_rule cTrp (eqF_intr th4)) 

646 
(reflexive l > mk_object_eq)) 

647 
end 

648 
end 

649 

650 
fun ring tm = 

651 
let 

652 
fun mk_forall x p = 

653 
mk_comb (Drule.cterm_rule (instantiate' [SOME (ctyp_of_term x)] []) @{cpat "All:: (?'a => bool) => _"}) (Thm.cabs x p) 

23487  654 
val avs = Thm.add_cterm_frees tm [] 
23252  655 
val P' = fold mk_forall avs tm 
656 
val th1 = initial_conv(mk_neg P') 

657 
val (evs,bod) = strip_exists(concl th1) in 

658 
if is_forall bod then error "ring: nonuniversal formula" 

659 
else 

660 
let 

661 
val th1a = weak_dnf_conv bod 

662 
val boda = concl th1a 

663 
val th2a = refute_disj refute boda 

664 
val th2b = [mk_object_eq th1a, (th2a COMP notI) COMP PFalse'] MRS trans 

665 
val th2 = fold (fn v => fn th => (forall_intr v th) COMP allI) evs (th2b RS PFalse) 

666 
val th3 = equal_elim 

667 
(Simplifier.rewrite (HOL_basic_ss addsimps [not_ex RS sym]) 

668 
(th2 > cprop_of)) th2 

669 
in specl avs 

670 
([[[mk_object_eq th1, th3 RS PFalse'] MRS trans] MRS PFalse] MRS notnotD) 

671 
end 

672 
end 

673 
fun ideal tms tm ord = 

674 
let 

675 
val rawvars = fold_rev grobvars (tm::tms) [] 

676 
val vars = sort ord (distinct (fn (x,y) => (term_of x) aconv (term_of y)) rawvars) 

677 
val pols = map (grobify_term vars) tms 

678 
val pol = grobify_term vars tm 

679 
val cert = grobner_ideal vars pols pol 

23557  680 
in map (fn n => these (AList.lookup (op =) cert n) > holify_polynomial vars) 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
23579
diff
changeset

681 
(0 upto (length pols  1)) 
23252  682 
end 
683 
in (ring,ideal) 

684 
end; 

685 

686 

687 
fun find_term bounds tm = 

688 
(case term_of tm of 

689 
Const ("op =", T) $ _ $ _ => 

690 
if domain_type T = HOLogic.boolT then find_args bounds tm 

691 
else Thm.dest_arg tm 

692 
 Const ("Not", _) $ _ => find_term bounds (Thm.dest_arg tm) 

693 
 Const ("All", _) $ _ => find_body bounds (Thm.dest_arg tm) 

694 
 Const ("Ex", _) $ _ => find_body bounds (Thm.dest_arg tm) 

695 
 Const ("op &", _) $ _ $ _ => find_args bounds tm 

696 
 Const ("op ", _) $ _ $ _ => find_args bounds tm 

697 
 Const ("op >", _) $ _ $ _ => find_args bounds tm 

698 
 _ => raise TERM ("find_term", [])) 

699 
and find_args bounds tm = 

700 
let val (t, u) = Thm.dest_binop tm 

701 
in (find_term bounds t handle TERM _ => find_term bounds u) end 

702 
and find_body bounds b = 

703 
let val (_, b') = Thm.dest_abs (SOME (Name.bound bounds)) b 

704 
in find_term (bounds + 1) b' end; 

705 

23579  706 
fun ring_solve ctxt form = 
23252  707 
(case try (find_term 0 (* FIXME !? *)) form of 
708 
NONE => reflexive form 

709 
 SOME tm => 

710 
(case NormalizerData.match ctxt tm of 

711 
NONE => reflexive form 

712 
 SOME (res as (theory, {is_const, dest_const, mk_const, conv = ring_eq_conv})) => 

713 
fst (ring_and_ideal_conv theory 

23330
01c09922ce59
Conversion for computation on constants now depends on the context
chaieb
parents:
23259
diff
changeset

714 
dest_const (mk_const (Thm.ctyp_of_term tm)) (ring_eq_conv ctxt) 
01c09922ce59
Conversion for computation on constants now depends on the context
chaieb
parents:
23259
diff
changeset

715 
(semiring_normalize_wrapper ctxt res)) form)); 
23252  716 

23579  717 
fun ring_tac add_ths del_ths ctxt = 
718 
ObjectLogic.full_atomize_tac 

719 
THEN' simp_tac (fst (NormalizerData.get ctxt) delsimps del_ths addsimps add_ths) 

720 
THEN' CSUBGOAL (fn (p, i) => 

721 
rtac (ring_solve ctxt (ObjectLogic.dest_judgment p)) i 

722 
handle TERM _ => no_tac 

723 
 CTERM _ => no_tac 

724 
 THM _ => no_tac); 

23334  725 

23252  726 
end; 