1(* Copyright (c) Michael Norrish *) 2 3structure mlibOmega :> mlibOmega = 4struct 5 6(* ---------------------------------------------------------------------- 7 Implementation of "shadow chasing" in ML. This procedure can find 8 a real shadow refutation without forcing all of the shadow to be 9 computed in the logic; a path to false can be found that only 10 involves a small proportion of the possible logical inferences. 11 12 Alternatively, if the shadow is exact, then this procedure can 13 also generate witnesses for the logic to exploit, again saving it 14 from having to do a lot of calculation internally. 15 ---------------------------------------------------------------------- *) 16 17type arbint = mlibOmegaint.int 18type 'a table = 'a mlibPatricia.t 19infix == 20val op == = mlibOmegaint.eq 21 22(* ---------------------------------------------------------------------- 23 The factoid datatype 24 25 A factoid c[0..n] represents the term 26 0 <= c[0] * v_0 + .. + c[n-1] * v_{n-1} + c[n] * 1 27 The factoid "key" is the vector of variable coefficients, leaving 28 out the constant part. 29 ---------------------------------------------------------------------- *) 30 31 32datatype factoid = 33 ALT of arbint Vector.vector 34 35structure V = Vector 36structure VS = VectorSlice 37structure AS = ArraySlice 38 39abstype factoid_key = KEY of arbint V.vector * int 40with 41 local open mlibOmegaint 42 val alpha = fromInt 33 43 val beta = fromInt 720 44 in 45 fun khash v len = 46 hash(if zero <= V.sub(v, 0) 47 then V.foldli (fn(_,i,h) => alpha * h + beta + i) zero v 48 else ~(V.foldli (fn(_,i,h)=> alpha * h + beta + ~i) zero v)) 49 end 50 51 fun factoid_key (ALT av) = KEY(av, V.length av - 1) 52 53 fun keyhash (KEY(av, len)) = khash av (SOME len) 54 55 fun map_key f (KEY(v, len)) = 56 VS.mapi (fn(_, x) => f x) (VS.slice(v, 0, SOME len)) 57 58 fun foldl_key f c (KEY(v, len)) = 59 VS.foldli (fn(_, x, c) => f(x, c)) c (VS.slice(v, 0, SOME len)) 60 61 fun foldli_key f c (KEY(v, len)) = VS.foldli f c (VS.slice(v, 0, SOME len)) 62 63 fun appi_key f (KEY(v, len)) = VS.appi f (VS.slice(v, 0, SOME len)) 64 65 fun sub_key (KEY(v, len), i) = if i < len then V.sub(v, i) 66 else raise Subscript 67 68 fun negate_key (KEY(v, len)) = KEY(V.map mlibOmegaint.~ v, len) 69 70 fun key_eq (KEY(v1, len1)) (KEY(v2, len2)) = 71 len1 = len2 72 andalso 73 VS.foldli (fn(i,e1,res) => res andalso e1 == V.sub(v2,i)) true 74 (VS.slice(v1, 0, SOME len1)) 75end 76 77 78fun factoid v = ALT v 79 80val fromArbList = factoid o V.fromList 81val fromList = fromArbList o map mlibOmegaint.fromInt 82fun toList (ALT av) = V.foldr op:: [] av 83 84 85 86fun factoid_constant f = 87 case f of 88 ALT av => V.sub(av, V.length av - 1) 89 90fun factoid_size f = 91 case f of ALT av => V.length av 92 93fun split_factoid f = (factoid_key f, factoid_constant f) 94 95fun zero_var_factoid f = 96 foldl_key (fn (ai, b) => b andalso ai == mlibOmegaint.zero) true 97 (factoid_key f) 98 99val the_false_factoid = factoid (V.fromList [mlibOmegaint.fromInt ~1]) 100fun false_factoid f = 101 zero_var_factoid f andalso mlibOmegaint.<(factoid_constant f, mlibOmegaint.zero) 102fun true_factoid f = 103 zero_var_factoid f andalso mlibOmegaint.<=(mlibOmegaint.zero, factoid_constant f) 104 105fun eval_factoid_rhs vmap f = let 106 val maxdex = factoid_size f - 1 107 open mlibOmegaint Vector 108 fun foldthis (i, ai, acc) = 109 if ai == zero then acc 110 else if i = maxdex then ai + acc 111 else ai * mlibPatricia.find i vmap + acc 112in 113 case f of 114 ALT fv => foldli foldthis zero fv 115end 116 117fun eval_factoid_except vmap i f = let 118 val maxdex = factoid_size f - 1 119 open mlibOmegaint Vector 120 fun foldthis (j, ai, acc) = 121 if ai == zero orelse i = j then acc 122 else if j = maxdex then ai + acc 123 else ai * mlibPatricia.find j vmap + acc 124in 125 case f of 126 ALT fv => foldli foldthis zero fv 127end 128 129 130(* "prints" factoids *) 131fun factoid_pairings vars f = 132 case f of 133 ALT av => V.foldri 134 (fn (i, ai, list) => 135 (mlibOmegaint.toString ai, List.nth(vars, i)) :: list) 136 [] av 137 138 139fun ERR f s = Fail ("Function: "^f^", message: "^s) 140 141fun gcdl l = 142 case l of 143 [] => raise ERR "gcdl" "empty list" 144 | (h::t) => foldl mlibOmegaint.gcd h t 145 146(*---------------------------------------------------------------------------*) 147(* Function to compute the Lowest Common Multiple of two integers. *) 148(*---------------------------------------------------------------------------*) 149 150fun lcm (i,j) = let open mlibOmegaint in (i * j) div gcd (i,j) end 151(*handle _ => raise ERR "lcm" "negative arguments to lcm"; 152*) 153(* ---------------------------------------------------------------------- 154 combine_real_factoids i f1 f2 155 156 takes two factoids and produces a fresh one by "variable 157 elimination" over the i'th variable in the vector. It requires 158 that the coefficient of v_i is strictly positive in f1, and 159 strictly negative in f2. The combination may produce a factoid where 160 the coefficients can all be divided through by some number, but 161 this factor won't be the gcd of the coefficients of v_i; this 162 factor is eliminated from the outset. 163 ---------------------------------------------------------------------- *) 164 165fun combine_real_factoids i f1 f2 = 166 case (f1, f2) of 167 (ALT av1, ALT av2) => let 168 open mlibOmegaint Vector 169 val c0 = sub(av1, i) 170 val d0 = ~ (sub(av2, i)) 171 val l = lcm (c0, d0) 172 val c = l div d0 173 val d = l div c0 174 fun gen j = c * sub(av2, j) + sub(av1, j) * d 175 in 176 factoid (tabulate(length av1, gen)) 177 end 178 179(* ---------------------------------------------------------------------- 180 combine_dark_factoids i lower upper 181 182 takes two factoids and combines them to produce a new, "dark shadow" 183 factoid. As above, the first one must have a positive coefficient of 184 i, and the second a negative coefficient. 185 ---------------------------------------------------------------------- *) 186 187fun combine_dark_factoids i low up = 188 case (low, up) of 189 (ALT L, ALT U) => let 190 open mlibOmegaint Vector 191 (* have ~L <= ax /\ bx <= U *) 192 val a = sub(L, i) 193 val b = ~(sub(U, i)) 194 val maxdex = Int.-(length L, 1) 195 fun gen j = let 196 val base = a * sub(U, j) + b * sub(L, j) 197 in 198 if j = maxdex then base - ((a - one) * (b - one)) 199 else base 200 end 201 in 202 factoid (tabulate(Int.+(maxdex, 1), gen)) 203 end 204 205 206(* ---------------------------------------------------------------------- 207 factoid_gcd f 208 209 eliminates common factors in the variable coefficients of a factoid, 210 or raises the exception no_gcd, if there is no common factor. 211 ---------------------------------------------------------------------- *) 212 213exception no_gcd 214fun factoid_gcd f = 215 case f of 216 ALT av => let 217 val g = VS.foldli 218 (fn (_, c, g0) => mlibOmegaint.gcd (mlibOmegaint.abs c, g0)) 219 mlibOmegaint.zero 220 (VS.slice(av, 0, SOME (V.length av - 1))) 221 in 222 if mlibOmegaint.<=(g, mlibOmegaint.one) then raise no_gcd 223 else factoid (V.map (fn c => mlibOmegaint.div(c, g)) av) 224 end 225 226(* ---------------------------------------------------------------------- 227 term_to_factoid vars tm 228 229 returns the factoid corresponding to tm. tm is thus of the form 230 0 <= c1 * v1 + ... + cn * vn + num 231 Assumes that the variables are in the "correct" order (as given in the 232 list vars), but that all are not necessarily present. Omission 233 indicates a coefficient of zero, of course. 234 ---------------------------------------------------------------------- *) 235 236(* 237fun term_to_factoid vars tm = let 238 val summands = strip_plus (rand tm) 239 fun mk_coeffs vlist slist = 240 case (vlist,slist) of 241 ([],[]) => [mlibOmegaint.zero] 242 | ([],[s]) => [int_of_term s] 243 | ([], _) => raise HOL_ERR { origin_function = "term_to_factoid", 244 origin_structure = "mlibOmegaMLShadow", 245 message = "Too many summands in term"} 246 | (v::vs, []) => mlibOmegaint.zero :: mk_coeffs vs [] 247 | (v::vs, s::ss) => 248 if is_mult s then let 249 val (c, mv) = dest_mult s 250 in 251 if mv = v then int_of_term c :: mk_coeffs vs ss 252 else mlibOmegaint.zero :: mk_coeffs vs slist 253 end 254 else mlibOmegaint.zero :: mk_coeffs vs slist 255in 256 ALT (V.fromList (mk_coeffs vars summands)) 257end 258*) 259 260(* ---------------------------------------------------------------------- 261 factoid_to_term vars f 262 263 returns the term corresponding to f, interpreting f over the list of 264 variables vars. 265 ---------------------------------------------------------------------- *) 266 267(* 268fun factoid_to_term vars f = 269 case f of 270 ALT av => let 271 open Vector 272 fun varcoeff (i, c, t0) = 273 if c = mlibOmegaint.zero then t0 274 else let 275 val cv = mk_mult(term_of_int c, List.nth(vars, i)) 276 in 277 case t0 of NONE => SOME cv | SOME t => SOME (mk_plus(t, cv)) 278 end 279 val varpart = foldli varcoeff NONE (av, 0, SOME (length av - 1)) 280 val numpart = term_of_int (sub(av,length av - 1)) 281 in 282 case varpart of 283 NONE => numpart 284 | SOME t => mk_plus(t, numpart) 285 end 286*) 287 288(* ---------------------------------------------------------------------- 289 The derivation datatype 290 ---------------------------------------------------------------------- *) 291 292(* a derivation represents a proof of a factoid *) 293datatype 'a derivation = 294 ASM of 'a 295 | REAL_COMBIN of int * 'a derivation * 'a derivation 296 | GCD_CHECK of 'a derivation 297 | DIRECT_CONTR of 'a derivation * 'a derivation 298type 'a dfactoid = (factoid * 'a derivation) 299 300datatype 'a result = CONTR of 'a derivation 301 | SATISFIABLE of mlibOmegaint.int table 302 | NO_CONCL 303 304fun direct_contradiction(d1, d2) = CONTR (DIRECT_CONTR(d1, d2)) 305 306fun gcd_check_dfactoid (df as (f, d)) = 307 (factoid_gcd f, GCD_CHECK d) handle no_gcd => df 308 309fun split_dfactoid (f, d) = (factoid_key f, (factoid_constant f, d)) 310fun dfactoid_key (f, d) = factoid_key f 311fun dfactoid_derivation (f, d) = d 312 313(* fun term_to_dfactoid vars t = (term_to_factoid vars t, ASM t) *) 314 315(* ---------------------------------------------------------------------- 316 The "elimination mode" datatype. 317 318 This records what sort of shadow we're currently working on. 319 ---------------------------------------------------------------------- *) 320 321datatype elimmode = REAL | DARK | EXACT | EDARK 322 (* REAL when we're looking for a contradiction (only) 323 DARK when we're looking for satisfiability and the problem is not 324 exact 325 EXACT when we're looking for either a contradiction or satisfiability. 326 EDARK when we're looking for satisfiability (i.e., have switched from 327 a REAL search, but where the problem is still EXACT) *) 328 329 330 331 332fun inexactify EXACT = REAL 333 | inexactify EDARK = DARK 334 | inexactify x = x 335 336fun mode_result em result = 337 case em of 338 EXACT => result 339 | REAL => (case result of SATISFIABLE _ => NO_CONCL | x => x) 340 | DARK => (case result of CONTR _ => NO_CONCL | x => x) 341 | EDARK => (case result of CONTR _ => NO_CONCL | x => x) 342 343fun combine_dfactoids em i ((f1, d1), (f2, d2)) = 344 case em of 345 DARK => (combine_dark_factoids i f1 f2, d1) 346 | _ => (combine_real_factoids i f1 f2, REAL_COMBIN(i, d1, d2)) 347 348 349(* ---------------------------------------------------------------------- 350 The "db" datatype 351 352 So far, I'm using a mlibPatricia tree to store my sets of constraints, 353 so the function parameters of type "database" are often called 354 ptree. The keys are the hashes of the constraint keys. The items 355 are lists (buckets) of dfactoids. 356 ---------------------------------------------------------------------- *) 357local structure V = Vector in 358fun veceq eq ( v1) ( v2) = 359 if V.length v1 <> V.length v2 then false 360 else V.foldli (fn(i,e1,res) => res andalso eq(e1, V.sub(v2,i))) true v1 361end 362 363fun aveceq x = veceq op== x; 364 365 366fun fkassoc k alist = 367 case List.find (fn (f, data) => key_eq k (factoid_key f)) alist of 368 NONE => raise mlibPatricia.NotFound 369 | SOME x => x 370 371fun lookup_fkey (ptree,w) fk = 372 fkassoc fk (mlibPatricia.find (keyhash fk) ptree) 373 374type 'a cstdb = ((factoid * 'a derivation) list) mlibPatricia.t * int 375 376fun dbempty i = (mlibPatricia.empty, i) 377 378fun dbfold (f : ((factoid * 'b derivation) * 'a) -> 'a) (acc:'a) (ptree,w) = 379 mlibPatricia.fold (fn (i,b,acc) => List.foldl f acc b) acc ptree 380 381fun dbapp (f : factoid * 'a derivation -> unit) (ptree,w) = 382 mlibPatricia.app (fn (i,b) => List.app f b) ptree 383 384(* does a raw insert, with no checking; the dfactoid really should have 385 width w *) 386fun dbinsert (df as (f,d)) (ptree,w) = 387 (mlibPatricia.addf (fn b => df :: b) (keyhash (factoid_key f)) [df] ptree, 388 w) 389 390fun dbchoose (ptree,w) = 391 case mlibPatricia.getItem ptree of 392 SOME((_, a_bucket),_) => List.hd a_bucket 393 | NONE => raise mlibPatricia.NotFound 394 395fun dbwidth (ptree,w) = w 396 397fun dbsize (ptree,w) = mlibPatricia.size ptree 398 399(* ---------------------------------------------------------------------- 400 dbadd df ptree 401 402 Precondition: df is neither a trivially true nor false factoid. 403 404 adds a dfactoid (i.e., a factoid along with its derivation) to a 405 database of accumulating factoids. If the factoid is 406 actually redundant, because there is already a factoid with the 407 same constraint key in the tree with a less or equal constant, 408 then the exception RedundantAddition is raised. Otherwise the 409 return value is the new tree. 410 ---------------------------------------------------------------------- *) 411 412exception Pluck 413fun pluck P = 414 let fun pl _ [] = raise Pluck 415 | pl A (h::t) = if P h then (h, List.foldl (op::) t A) else pl (h::A) t 416 in pl [] 417 end; 418 419exception RedundantAddition 420fun dbadd df (ptree, w) = let 421 fun merge alist = let 422 val (fk, (fc, _)) = split_dfactoid df 423 val ((f',_), samehashes) = 424 pluck (fn (f', _) => key_eq fk (factoid_key f')) alist 425 in 426 if mlibOmegaint.<(fc, factoid_constant f') then 427 df :: samehashes 428 else raise RedundantAddition 429 end handle Pluck => df :: alist 430in 431 (mlibPatricia.addf merge (keyhash (dfactoid_key df)) [df] ptree, w) 432end 433 434(* ---------------------------------------------------------------------- 435 add_check_factoid ptree df next kont 436 437 Precondition: df is neither a trivially true nor false factoid. 438 Types: 439 next: ptree -> (result -> 'a) -> 'a 440 kont: result -> 'a 441 442 ---------------------------------------------------------------------- *) 443 444fun add_check_factoid ptree0 (df as (f,d)) next kont = let 445 val ptree = dbadd df ptree0 446in 447 let 448 val (fk,fc) = split_factoid f 449 val (negdf as (negf, negd)) = lookup_fkey ptree (negate_key fk) 450 val negc = factoid_constant negf 451 open mlibOmegaint 452 in 453 if negc < ~fc then kont (direct_contradiction(d, negd)) 454 else next ptree kont 455 end handle mlibPatricia.NotFound => next ptree kont 456end handle RedundantAddition => next ptree0 kont 457 458 459(* ---------------------------------------------------------------------- 460 has_one_var ptree 461 462 returns true if all of ptree's factoids are over just one variable 463 (and they're all the same variable) 464 ---------------------------------------------------------------------- *) 465 466fun has_one_var ptree = let 467 exception return_false 468 fun checkkey (i,c,acc) = 469 if c == mlibOmegaint.zero then acc 470 else 471 case acc of 472 NONE => SOME i 473 | _ => raise return_false 474 fun check ((f,d),a) = let 475 val fk = factoid_key f 476 val fk_single = valOf (foldli_key checkkey NONE fk) 477 in 478 case a of 479 NONE => SOME fk_single 480 | SOME i => if i = fk_single then a else raise return_false 481 end 482in 483 (dbfold check NONE ptree; true) handle return_false => false 484end 485 486 487(* ---------------------------------------------------------------------- 488 one_var_analysis ptree em 489 490 Precondition: the dfactoids in ptree have just one variable with a 491 non-zero coefficient, and its everywhere the same variable. Our aim 492 is to either derive a contradiction, or to return a satisfying 493 assignment. Our gcd_checks will have ensured that our variable 494 (call it x) will only have coefficients of one at this point, so 495 we just need to check that the maximum of the lower bounds is less 496 than or equal to the minimum of the lower bounds. If it is then 497 return either as a satisfying valuation for x. Otherwise return 498 false, combining the maximum lower and the minimum upper constraint. 499 ---------------------------------------------------------------------- *) 500 501fun one_var_analysis ptree em = let 502 val (a_constraint, _) = dbchoose ptree 503 fun find_var (i, ai, NONE) = if not (ai == mlibOmegaint.zero) then SOME i else NONE 504 | find_var (_, _, v as SOME _) = v 505 val x_var = 506 valOf (foldli_key find_var NONE (factoid_key a_constraint)) 507 open mlibOmegaint 508 fun assign_factoid (df, acc as (upper, lower)) = let 509 val (fk as k, (fc, d)) = split_dfactoid df 510 in 511 if sub_key(k, x_var) < zero then 512 case upper of 513 NONE => (SOME (fc, d), lower) 514 | SOME (c', d') => if fc < c' then (SOME (fc, d), lower) else acc 515 else 516 case lower of 517 NONE => (upper, SOME (~fc, d)) 518 | SOME (c', d') => if ~fc > c' then (SOME (~fc,d), lower) else acc 519 end 520 val (upper, lower) = dbfold assign_factoid (NONE,NONE) ptree 521in 522 case (upper, lower) of 523 (NONE, NONE) => raise Fail "one_var_analysis: this can't happen" 524 | (SOME (c, _), NONE) => if em = REAL then NO_CONCL 525 else SATISFIABLE (mlibPatricia.add x_var c mlibPatricia.empty) 526 | (NONE, SOME (c, _)) => if em = REAL then NO_CONCL 527 else SATISFIABLE (mlibPatricia.add x_var c mlibPatricia.empty) 528 | (SOME (u,du), SOME (l, dl)) => 529 if u < l then 530 if em = DARK orelse em = EDARK then NO_CONCL 531 else direct_contradiction(du,dl) 532 else 533 if em = REAL then NO_CONCL 534 else SATISFIABLE (mlibPatricia.add x_var u mlibPatricia.empty) 535end 536 537(* ---------------------------------------------------------------------- 538 throwaway_redundant_factoids ptree nextstage kont 539 540 checks ptree for variables that are constrained only in one sense 541 (i.e., with upper or lower bounds only). 542 543 The function nextstage takes a ptree and a continuation; it is 544 called when ptree has run out of constraints that can be thrown 545 away. 546 547 The continuation function kont is of type result -> 'a, and will 548 be called when the whole process eventually gets an answer. This 549 code will not call it directly, but if it does throw anything 550 away, it will modify it so that a satisfying value can be 551 calculated for the variables that are chucked. 552 ---------------------------------------------------------------------- *) 553 554fun throwaway_redundant_factoids ptree nextstage kont = let 555 val dwidth = dbwidth ptree 556 val numvars = dwidth - 1 557 val has_low = Array.array(numvars, false) 558 val has_up = Array.array(numvars, false) 559 560 fun assess fk = let 561 fun appthis (i, ai) = let 562 open mlibOmegaint 563 in 564 case compare(ai,zero) of 565 LESS => Array.update(has_up, i, true) 566 | EQUAL => () 567 | GREATER => Array.update(has_low, i, true) 568 end 569 in 570 appi_key appthis fk 571 end 572 val () = dbapp (fn (f,d) => assess (factoid_key f)) ptree 573 574 fun find_redundant_var i = 575 if i = numvars then NONE 576 else if Array.sub(has_up, i) = Array.sub(has_low, i) then 577 find_redundant_var (i + 1) 578 else SOME (i, Array.sub(has_up, i)) 579in 580 case find_redundant_var 0 of 581 NONE => nextstage ptree kont 582 | SOME(i, hasupper) => let 583 fun partition (df as (f,d), (pt, elim)) = let 584 open mlibOmegaint 585 val fk = factoid_key f 586 in 587 if sub_key(fk, i) == zero then 588 (dbinsert df pt, elim) 589 else 590 (pt, df :: elim) 591 end 592 val (newptree, elim) = dbfold partition (dbempty dwidth, []) ptree 593 fun handle_result r = 594 case r of 595 SATISFIABLE vmap => let 596 open mlibOmegaint 597 fun mapthis (f,_) = 598 let val fk = factoid_key f 599 in (eval_factoid_except vmap i f) div 600 abs (sub_key(fk, i)) 601 end 602 val evaluated = if hasupper then map mapthis elim 603 else map (~ o mapthis) elim 604 val foldfn = if hasupper then mlibOmegaint.min else mlibOmegaint.max 605 in 606 SATISFIABLE (mlibPatricia.add i (foldl foldfn (hd evaluated) 607 (tl evaluated)) 608 vmap) 609 end 610 | x => x 611 in 612 throwaway_redundant_factoids newptree nextstage (kont o handle_result) 613 end 614end 615 616(* ---------------------------------------------------------------------- 617 exact_var ptree 618 619 An exact var is one that has coefficients of one in either all of its 620 upper bounds or all of its lower bounds. This function returns 621 SOME v if v is an exact var in ptree, or NONE if there is no exact 622 var. 623 ---------------------------------------------------------------------- *) 624 625fun exact_var ptree = let 626 val up_coeffs_unit = Array.array(dbwidth ptree - 1, true) 627 val low_coeffs_unit = Array.array(dbwidth ptree - 1, true) 628 val coeffs_all_zero = Array.array(dbwidth ptree - 1, true) 629 fun appthis (f,d) = let 630 open mlibOmegaint 631 fun examine_key (i, ai) = 632 case compare(ai,zero) of 633 LESS => (if not(abs ai == one) then 634 Array.update(up_coeffs_unit, i, false) 635 else (); 636 Array.update(coeffs_all_zero, i, false)) 637 | EQUAL => () 638 | GREATER => (if not(ai == one) then 639 Array.update(low_coeffs_unit, i, false) 640 else (); 641 Array.update(coeffs_all_zero, i, false)) 642 val fk = factoid_key f 643 in 644 appi_key examine_key fk 645 end 646 val () = dbapp appthis ptree 647 fun check_index (i, b, acc) = 648 if b andalso not (Array.sub(coeffs_all_zero,i)) then SOME i else acc 649in 650 case Array.foldli check_index NONE low_coeffs_unit of 651 NONE => Array.foldli check_index NONE up_coeffs_unit 652 | x => x 653end 654 655(* ---------------------------------------------------------------------- 656 least_coeff_var ptree 657 658 Returns the variable whose coefficients' absolute values sum to the 659 least amount (that isn't zero). 660 ---------------------------------------------------------------------- *) 661 662fun least_coeff_var ptree = let 663 val sums = Array.array(dbwidth ptree - 1, mlibOmegaint.zero) 664 fun appthis df = let 665 open mlibOmegaint 666 fun add_in_c (i, ai) = 667 Array.update(sums, i, Array.sub(sums, i) + abs ai) 668 in 669 appi_key add_in_c (dfactoid_key df) 670 end 671 val () = dbapp appthis ptree 672 fun find_min (i,ai,NONE) = if not(ai == mlibOmegaint.zero) then SOME(i,ai) else NONE 673 | find_min (i,ai, acc as SOME(mini, minai)) = 674 if mlibOmegaint.<(ai, minai) andalso not(ai == mlibOmegaint.zero) then SOME(i,ai) else acc 675in 676 #1 (valOf (Array.foldli find_min NONE sums)) 677end 678 679(* ---------------------------------------------------------------------- 680 generate_row (pt0, em, i, up, lows, next, kont) 681 682 Types: 683 pt0 :ptree 684 em :elimmode 685 i :int (a variable index) 686 up :dfactoid 687 lows :dfactoid list 688 next :ptree -> (result -> 'a) -> 'a 689 kont :result -> 'a 690 691 "Resolves" dfactoid against all the factoids in lows, producing new 692 factoids, which get added to the ptree. If a factoid is directly 693 contradictory, then return it immediately, using kont. If a factoid 694 is vacuous, drop it. 695 ---------------------------------------------------------------------- *) 696 697fun generate_row pt0 em i up lows next kont = 698 case lows of 699 [] => next pt0 kont 700 | low::lowtl => let 701 val (df as (f,d)) = 702 gcd_check_dfactoid (combine_dfactoids em i (low, up)) 703 fun after_add pt k = generate_row pt em i up lowtl next k 704 in 705 if true_factoid f then after_add pt0 kont 706 else if false_factoid f then kont (CONTR d) 707 else add_check_factoid pt0 df after_add kont 708 end 709 710(* ---------------------------------------------------------------------- 711 generate_crossproduct (pt0, em, i, ups, lows, next, kont) 712 713 Types: 714 pt0 :ptree 715 em :elimmode 716 i :integer (a variable index) 717 ups :dfactoid list 718 lows :dfactoid list 719 next :ptree -> (result -> 'a) -> 'a 720 kont :result -> 'a 721 ---------------------------------------------------------------------- *) 722 723fun generate_crossproduct (pt0, em, i, ups, lows, next, kont) = 724 case ups of 725 [] => next pt0 kont 726 | u::us => let 727 fun after pt k = generate_crossproduct(pt, em, i, us, lows, next, k) 728 in 729 generate_row pt0 em i u lows after kont 730 end 731 732(* ---------------------------------------------------------------------- 733 extend_vmap ptree i vmap 734 735 vmap provides values for all of the variables present in ptree's 736 factoids except i. Use it to evaluate all of the factoids, except 737 at variable i and to then return vmap extended with a value for 738 variable i that respects all of the factoids. 739 740 Note definite parallels with code in one_var_analysis. 741 ---------------------------------------------------------------------- *) 742 743fun extend_vmap ptree i vmap = let 744 fun categorise ((f, _), (acc as (lower, upper))) = let 745 val c0 = eval_factoid_except vmap i f 746 val (fk as k) = factoid_key f 747 open mlibOmegaint 748 val coeff = sub_key(k, i) 749 in 750 case compare(coeff, zero) of 751 LESS => let (* upper *) 752 val c = c0 div ~coeff 753 in 754 case upper of 755 NONE => (lower, SOME c) 756 | SOME c' => if c < c' then (lower, SOME c) else acc 757 end 758 | EQUAL => acc 759 | GREATER => let (* lower *) 760 val c = ~(c0 div coeff) 761 in 762 case lower of 763 NONE => (SOME c, upper) 764 | SOME c' => if c > c' then (SOME c, upper) else acc 765 end 766 end 767 val (lower, upper) = dbfold categorise (NONE, NONE) ptree 768 open mlibOmegaint 769 val _ = valOf lower <= valOf upper orelse raise Fail "urk in extend map" 770in 771 mlibPatricia.add i (valOf lower) vmap 772end 773 774(* ---------------------------------------------------------------------- 775 zero_upto n 776 777 Returns an integer map that maps the integers from 0..n (inclusive) 778 to mlibOmegaint.zero. 779 ---------------------------------------------------------------------- *) 780 781fun zero_upto n = 782 if n < 0 then raise Fail "zero_upto called with negative argument" 783 else if n = 0 then mlibPatricia.add 0 mlibOmegaint.zero mlibPatricia.empty 784 else mlibPatricia.add n mlibOmegaint.zero (zero_upto (n - 1)) 785 786 787(* ---------------------------------------------------------------------- 788 one_step ptree em nextstage kont 789 790 Types: 791 em :elimmode 792 nextstage :ptree -> elimmode -> (result -> 'a) -> 'a 793 kont :result -> 'a 794 795 Assume that ptree doesn't contain anything directly contradictory, 796 and that there aren't any redundant constraints around (these have 797 been thrown away by throwaway_redundant_factoids). 798 799 Perform one step of the shadow calculation algorithm: 800 (a) if there is only one variable left in all of the constraints 801 in ptree do one_var_analysis to return some sort of result, 802 and pass this result to kont. 803 (b) otherwise, decide on a variable to eliminate. If there's a 804 variable that allows for an exact shadow (has coefficients of 805 one in all of its lower-bound or upper-bound occurrences), 806 pick this. Otherwise, take the variable whose coefficients' 807 absolute values sum to the least amount. 808 (c) the initial tree of the result contains all of those 809 constraints from the original which do not mention the 810 variable to be eliminated. 811 (d) divide the constraints that do mention the variable to be 812 eliminated into upper and lower bound sets. 813 Then (in generate_crossproduct): 814 (e) work through all of the possible combinations in 815 upper x lower 816 using combine_dfactoids. Each new dfactoid has to be added 817 to the accumulating tree. This may produce a direct 818 contradiction, in which case, stop and return this (again, 819 using the kont function). Otherwise keep going. 820 (f) pass the completed new ptree to nextstage, with an augmented 821 continuation that copes with satisfiable results (though 822 satisfiable results won't come back if the mode is REAL 823 of course). 824 ---------------------------------------------------------------------- *) 825 826fun one_step ptree em nextstage kont = let 827 val (var_to_elim, mode) = 828 case exact_var ptree of 829 SOME i => (i, em) 830 | NONE => (least_coeff_var ptree, inexactify em) 831 fun categorise (df, (notmentioned, uppers, lowers)) = let 832 open mlibOmegaint 833 val fk = dfactoid_key df 834 in 835 case compare(sub_key(fk, var_to_elim), zero) of 836 LESS => (notmentioned, df::uppers, lowers) 837 | EQUAL => (dbinsert df notmentioned, uppers, lowers) 838 | GREATER => (notmentioned, uppers, df::lowers) 839 end 840 val (newtree0, uppers, lowers) = 841 dbfold categorise (dbempty (dbwidth ptree), [], []) ptree 842 fun drop_contr (CONTR _) = NO_CONCL 843 | drop_contr x = x 844 fun extend_satisfiable r = 845 case r of 846 SATISFIABLE vmap => 847 SATISFIABLE (extend_vmap ptree var_to_elim vmap) 848 | _ => r 849 fun newkont r = 850 case (em, mode) of 851 (EXACT, EXACT) => kont (extend_satisfiable r) 852 | (EXACT, REAL) => let 853 in 854 case r of 855 CONTR _ => kont r 856 | _ => one_step ptree EDARK nextstage kont 857 end 858 | (REAL, REAL) => kont r 859 | (EDARK, EDARK) => kont (drop_contr (extend_satisfiable r)) 860 | (EDARK, DARK) => kont (drop_contr (extend_satisfiable r)) 861 | (DARK, DARK) => kont (drop_contr (extend_satisfiable r)) 862 | _ => raise Fail "Can't happen - in newkont calculation" 863 fun newnextstage pt k = nextstage pt mode k 864in 865 generate_crossproduct(newtree0, mode, var_to_elim, uppers, lowers, 866 newnextstage, newkont) 867end 868 869(* ---------------------------------------------------------------------- 870 toplevel ptree em kont 871 872 ---------------------------------------------------------------------- *) 873 874fun toplevel ptree em kont = let 875 fun after_throwaway pt k = 876 if dbsize pt = 0 then 877 k (SATISFIABLE (zero_upto (dbwidth pt - 2))) 878 else if has_one_var ptree then 879 k (mode_result em (one_var_analysis ptree em)) 880 else 881 one_step pt em toplevel k 882in 883 throwaway_redundant_factoids ptree after_throwaway kont 884end 885 886(* ---------------------------------------------------------------------- 887 work ptree kont 888 889 ---------------------------------------------------------------------- *) 890 891fun work ptree k = toplevel ptree EXACT k 892 893(* ---------------------------------------------------------------------- 894 simple tests 895 896 open mlibOmegaMLShadow 897 fun failkont r = raise Fail "kont" 898 fun next1 pt em k = pt 899 900 fun fromList l = (mlibOmegaMLShadow.fromList l, ASM l) 901 902 datatype 'a prettyr = SAT of bool * (string * mlibOmegaint.int) list 903 | PCONTR of 'a derivation 904 | PNOCONC 905 fun v2list vs = #1 (List.foldr (fn(v,(acc,n)) => 906 (("X"^Int.toString n, v)::acc, n - 1)) 907 ([], length vs - 1) vs) 908 909 fun display_result orig r = 910 case r of 911 SATISFIABLE v => let 912 val limit = length (hd orig) 913 val satlist = List.tabulate(limit, 914 (fn i => if i = limit - 1 then 915 mlibOmegaint.one 916 else mlibPatricia.find i v)) 917 val vslistlist = map (fn c => ListPair.zip(c, satlist)) orig 918 open mlibOmegaint 919 val sumprod = List.foldl (fn((x,y),acc) => acc + fromInt x * y) zero 920 in 921 SAT (List.all (fn l => zero <= sumprod l) vslistlist, 922 v2list satlist) 923 end 924 | CONTR x => PCONTR x 925 | NO_CONCL => PNOCONC; 926 927 fun dbinsert db f = 928 add_check_factoid db f (fn db => fn g => db) failkont 929 fun gentree l = 930 List.foldl (fn (df,t) => dbinsert t (fromList df)) 931 (dbempty (length (hd l))) l; 932 933 val seed = Random.newgen() 934 fun gen_int() = Random.range(~15,16) seed 935 fun gen_constraint n = List.tabulate(n, fn n => gen_int()) 936 fun gen_test m n = List.tabulate(m, fn _ => gen_constraint n) 937 938 939 fun test l = work (gentree l) (display_result l); 940 941 fun test csts = let 942 fun add_csts db normalc exc csts = 943 case csts of 944 [] => normalc db 945 | (c::cs) => 946 add_check_factoid db (fromList c) 947 (fn db => fn k => add_csts db normalc k cs) 948 exc 949 in 950 add_csts (dbempty (length (hd csts))) work display_result 951 end 952 953 test [[2,3,6],[~1,~4,7]]; (* exact elimination *) 954 test [[2,3,4],[~3,~4,7]]; (* dark shadow test *) 955 test [[2,3,4],[~3,~4,7],[4,5,~10]]; (* another dark shadow *) 956 test [[2,3,4],[~3,~4,7],[4,~5,~10]]; (* also satisfiable *) 957 test [[~3,~1,6],[2,1,~5],[3,~1,~3]]; (* a contradiction *) 958 test [[11,13,~27], [~11,~13,45],[7,~9,10],[~7,9,4]]; (* no conclusion *) 959 960 (* satisfiable after throwing away everything because of var 1 *) 961 test [[1,2,3,4],[2,1,4,3],[5,6,7,~8],[~3,2,~1,6]]; 962 963 test [[1,2,3,4],[2,~2,3,~10],[2,3,~5,6],[~3,~2,1,7]]; 964 965 test [[~9, ~11, ~8, 9, 11], [15, 0, 8, ~7, 8], [4, 3, 11, ~2, ~13], 966 [~13, ~8, ~14, 15, 8], [~10, 9, 15, ~13, 9], [~15, ~14, ~3, 2, 5], 967 [2, ~1, ~5, 14, ~7]]; 968 969 970 (* very poor performance exhibited here: *) 971 test [[13, ~5, ~12, 11, 1, ~5, ~3, ~6], [14, ~6, ~8, ~1, 4, ~10, 15, 6], 972 [9, ~5, 1, 10, 2, ~1, ~2, 0], [~1, ~13, 14, ~3, 11, ~9, 14, 9], 973 [~13, 1, 5, ~14, ~6, ~3, 14, 12], [11, 8, 9, 12, 1, ~2, ~6, 8], 974 [~14, ~8, ~4, ~2, ~3, 13, ~7, ~10], [~3, 3, ~14, ~3, ~7, 4, ~6, ~6]]; 975 976 fun print_test (l:int list list) = ignore (printVal l) 977 978 val current_goal = ref ([] : int list list) 979 val slow_goals = ref ([] : int list list list) 980 fun sattest i timelimit numcsts numvars = if i <= 0 then () 981 else let 982 val t = gen_test numcsts numvars 983 val _ = current_goal := t (* OK *) 984 val ctimer = Timer.startCPUTimer() 985 val result = test t 986 val timetaken = Timer.checkCPUTimer ctimer 987 val _ = if Time.>(#usr timetaken, 988 Time.fromSeconds timelimit) then 989 slow_goals := t :: (!slow_goals) (* OK *) 990 else () 991 in 992 case result of 993 SAT (b, _) => if b then print "SAT - OK\n" 994 else (print "SAT - FAILED\n"; 995 print_test t; 996 print "SAT - FAILED\n\n") 997 | PCONTR _ => print "CONTR\n" 998 | PNOCONC => print "PNOCONC\n"; 999 sattest (i - 1) timelimit numcsts numvars 1000 end 1001 1002 1003 fun can_throwaway ptree = let 1004 val next = throwaway_redundant_factoids ptree (fn p => fn k => p) failkont 1005 in 1006 dbsize next < dbsize ptree 1007 end 1008 1009 fun dbonestep ptree em = 1010 throwaway_redundant_factoids ptree 1011 (fn p => 1012 fn k => 1013 one_step ptree em 1014 (fn pt => 1015 fn em => 1016 fn k => (pt, em)) 1017 k) failkont 1018 1019 ---------------------------------------------------------------------- *) 1020 1021 1022end (* struct *) 1023