1 module LinearEquationsCoeffs
11 clone algebra.OrderedUnitaryCommutativeRing as A with
12 type t = a, function (+) = (+), function (*) = (*),
13 function (-_) = (-_), constant zero = azero,
14 constant one=aone, predicate (<=) = ale,
19 axiom sub_def: forall a1 a2. a1 - a2 = a1 + (- a2)
26 function interp t cvars : a
28 val constant czero : t
31 axiom zero_def: forall y. interp czero y = azero
32 axiom one_def: forall y. interp cone y = aone
35 forall x y: a. (-x) * y = - (x*y)
38 ensures { forall v: cvars. interp result v = interp a v + interp b v }
39 raises { Unknown -> true }
42 ensures { forall v: cvars. interp result v = interp a v * interp b v }
43 raises { Unknown -> true }
46 ensures { forall v: cvars. interp result v = - (interp a v) }
48 val predicate eq (a b:t)
49 ensures { result -> forall y:cvars. interp a y = interp b y }
52 requires { not (eq a czero) }
53 (* ensures { forall v: cvars. interp result v * interp a v = aone } no proof needed, but had better be true *)
54 ensures { not (eq result czero) }
55 raises { Unknown -> true }
59 module LinearEquationsDecision
64 clone LinearEquationsCoeffs as C with type t = coeff, axiom .
67 type expr = Term coeff int | Add expr expr | Cst coeff
69 let rec predicate valid_expr (e:expr)
74 | Add e1 e2 -> valid_expr e1 && valid_expr e2
77 let rec predicate expr_bound (e:expr) (b:int)
80 | Term _ i -> 0 <= i <= b
82 | Add e1 e2 -> expr_bound e1 b && expr_bound e2 b
85 function interp (e:expr) (y:vars) (z:C.cvars) : C.a
87 | Term c v -> C.(*) (C.interp c z) (y v)
88 | Add e1 e2 -> C.(+) (interp e1 y z) (interp e2 y z)
89 | Cst c -> C.interp c z
95 type equality = (expr, expr)
96 type context = list equality
98 let predicate valid_eq (eq:equality)
99 = match eq with (e1,e2) -> valid_expr e1 && valid_expr e2 end
101 let predicate eq_bound (eq:equality) (b:int)
102 = match eq with (e1,e2) -> expr_bound e1 b && expr_bound e2 b end
104 let rec predicate valid_ctx (ctx:context)
105 = match ctx with Nil -> true | Cons eq t -> valid_eq eq && valid_ctx t end
107 let rec predicate ctx_bound (ctx:context) (b:int)
108 = match ctx with Nil -> true | Cons eq t -> eq_bound eq b && ctx_bound t b end
110 let rec lemma expr_bound_w (e:expr) (b1 b2:int)
111 requires { b1 <= b2 }
112 requires { expr_bound e b1 }
113 ensures { expr_bound e b2 }
116 | Add e1 e2 -> expr_bound_w e1 b1 b2; expr_bound_w e2 b1 b2
121 lemma eq_bound_w: forall e:equality, b1 b2:int. eq_bound e b1 -> b1 <= b2 -> eq_bound e b2
123 let rec lemma ctx_bound_w (l:context) (b1 b2:int)
124 requires { ctx_bound l b1 }
125 requires { b1 <= b2 }
126 ensures { ctx_bound l b2 }
128 = match l with Nil -> () | Cons _ t -> ctx_bound_w t b1 b2 end
130 function interp_eq (g:equality) (y:vars) (z:C.cvars) : bool
131 = match g with (g1, g2) -> interp g1 y z = interp g2 y z end
133 function interp_ctx (l: context) (g: equality) (y: vars) (z:C.cvars) : bool
135 | Nil -> interp_eq g y z
136 | Cons h t -> (interp_eq h y z) -> (interp_ctx t g y z)
141 use mach.array.Array63
142 use mach.matrix.Matrix63
144 let apply_r (m: matrix coeff) (v: array coeff) : array coeff
145 requires { v.length = m.columns }
146 ensures { result.length = m.rows }
147 raises { C.Unknown -> true }
148 = let r = Array63.make m.rows C.czero in
149 for i = 0 to m.rows - 1 do
150 for j = 0 to m.columns - 1 do
151 r[i] <- C.add r[i] (C.mul (get m i j) v[j]);
156 let apply_l (v: array coeff) (m: matrix coeff) : array coeff
157 requires { v.length = m.rows }
158 ensures { result.length = m.columns }
159 raises { C.Unknown -> true }
160 = let r = Array63.make m.columns C.czero in
161 for j = 0 to m.columns - 1 do
162 for i = 0 to m.rows - 1 do
163 r[j] <- C.add r[j] (C.mul (get m i j) v[i]);
170 let sprod (a b: array coeff) : coeff
171 requires { a.length = b.length }
172 raises { C.Unknown -> true }
173 = let r = ref C.czero in
174 for i = 0 to a.length - 1 do
175 r := C.add !r (C.mul a[i] b[i]);
179 let m_append (m: matrix coeff) (v:array coeff) : matrix coeff
180 requires { m.rows = v.length }
181 requires { m.columns < int63'maxInt }
182 ensures { result.rows = m.rows }
183 ensures { result.columns = m.columns + 1 }
184 ensures { forall i j. 0 <= i < m.rows -> 0 <= j < m.columns ->
185 result.elts i j = m.elts i j }
186 ensures { forall i. 0 <= i < m.rows -> result.elts i m.columns = v[i] }
187 = let r = Matrix63.make m.rows (m.columns + 1) C.czero in
188 for i = 0 to m.rows - 1 do
189 invariant { forall k j. 0 <= k < i -> 0 <= j < m.columns ->
190 r.elts k j = m.elts k j }
191 invariant { forall k. 0 <= k < i -> r.elts k m.columns = v[k] }
192 for j = 0 to m.columns - 1 do
193 invariant { forall k j. 0 <= k < i -> 0 <= j < m.columns ->
194 r.elts k j = m.elts k j }
195 invariant { forall k. 0 <= k < i -> r.elts k m.columns = v[k] }
196 invariant { forall l. 0 <= l < j -> r.elts i l = m.elts i l }
197 set r i j (get m i j)
199 set r i m.columns v[i]
203 let v_append (v: array coeff) (c: coeff) : array coeff
204 requires { length v < int63'maxInt }
205 ensures { length result = length v + 1 }
206 ensures { forall k. 0 <= k < v.length -> result[k] = v[k] }
207 ensures { result[v.length] = c }
208 = let r = Array63.make (v.length + 1) c in
209 for i = 0 to v.length - 1 do
210 invariant { forall k. 0 <= k < i -> r[k] = v[k] }
211 invariant { r[v.length] = c }
216 let predicate (==) (a b: array coeff)
217 ensures { result = true -> length a = length b /\
218 forall i. 0 <= i < length a -> C.eq a[i] b[i] }
220 if length a <> length b then false
223 for i = 0 to length a - 1 do
224 invariant { !r = true -> forall j. 0 <= j < i -> C.eq a[j] b[j] }
225 if not (C.eq a[i] b[i]) then r := false;
232 let rec function max_var (e:expr) : int
234 requires { valid_expr e }
235 ensures { 0 <= result }
236 ensures { expr_bound e result }
240 | Add e1 e2 -> max (max_var e1) (max_var e2)
243 let function max_var_e (e:equality) : int
244 requires { valid_eq e }
245 ensures { 0 <= result }
246 ensures { eq_bound e result }
247 = match e with (e1,e2) -> max (max_var e1) (max_var e2) end
249 let rec function max_var_ctx (l:context) : int
251 requires { valid_ctx l }
252 ensures { 0 <= result }
253 ensures { ctx_bound l result }
256 | Cons e t -> max (max_var_e e) (max_var_ctx t)
259 let rec opp_expr (e:expr) : expr
260 ensures { forall y z. interp result y z = C.(-_) (interp e y z) }
261 ensures { valid_expr e -> valid_expr result }
262 ensures { forall b. expr_bound e b -> expr_bound result b }
265 | Cst c -> Cst (C.opp c)
269 assert { forall y z. interp r y z = C.(*) (C.interp oc z) (y j)
270 = C.(*) (C.(-_) (C.interp c z)) (y j)
271 = C.(-_) (C.(*) (C.interp c z) (y j))
272 = C.(-_) (interp e y z) };
275 let e1' = opp_expr e1 in
276 let e2' = opp_expr e2 in
277 assert { forall a1 a2. C.(+) (C.(-_) a1) (C.(-_) a2) = C.(-_) (C.(+) a1 a2) };
278 assert { forall y z. interp (Add e1' e2') y z = C.(-_) (interp e y z) by
279 interp (Add e1' e2') y z = C.(+) (interp e1' y z) (interp e2' y z)
280 = C.(+) (C.(-_) (interp e1 y z)) (C.(-_) (interp e2 y z))
281 = C.(-_) (C.(+) (interp e1 y z) (interp e2 y z))
282 = C.(-_) (interp e y z) };
286 predicate atom (e:expr)
288 | Add _ _ -> false | _ -> true
291 (*TODO put this back in norm_eq*)
292 let rec norm_eq_aux (ex acc_e:expr) (acc_c:coeff) : (expr, coeff)
293 returns { (rex, rc) -> forall y z.
294 C.(+) (interp rex y z) (interp (Cst rc) y z)
295 = C.(+) (interp ex y z)
296 (C.(+) (interp acc_e y z) (interp (Cst acc_c) y z)) }
297 returns { (rex, _) -> forall b:int. expr_bound ex b /\ expr_bound acc_e b
298 -> expr_bound rex b }
299 raises { C.Unknown -> true }
302 | Cst c -> acc_e, (C.add c acc_c)
303 | Term _ _ -> (Add acc_e ex, acc_c)
304 | Add e1 e2 -> let ae, ac = norm_eq_aux e1 acc_e acc_c in
310 let norm_eq (e:equality) : (expr, coeff)
311 returns { (ex, c) -> forall y z.
312 interp_eq e y z <-> interp_eq (ex, Cst c) y z }
313 returns { (ex, _) -> forall b:int. eq_bound e b -> expr_bound ex b }
314 raises { C.Unknown -> true }
317 let s = Add e1 (opp_expr e2) in
318 assert { forall b. eq_bound e b -> expr_bound s b };
319 match norm_eq_aux s (Cst C.czero) C.czero with
322 assert { forall a1 a2. C.(+) a1 a2 = C.azero -> a1 = C.(-_) a2 };
323 assert { forall y z. interp_eq (e1,e2) y z -> interp_eq (e, Cst ec) y z
324 by interp_eq (s, Cst C.czero) y z so interp s y z = C.azero
325 so C.(+) (interp e y z) (interp (Cst c) y z) = C.azero
326 so interp e y z = C.(-_) (interp (Cst c) y z)
327 = interp (Cst ec) y z };
333 let rec lemma interp_ctx_impl (ctx: context) (g1 g2:equality)
334 requires { forall y z. interp_eq g1 y z -> interp_eq g2 y z }
335 ensures { forall y z. interp_ctx ctx g1 y z -> interp_ctx ctx g2 y z }
337 = match ctx with Nil -> () | Cons _ t -> interp_ctx_impl t g1 g2 end
339 let rec lemma interp_ctx_valid (ctx:context) (g:equality)
340 ensures { forall y z. interp_eq g y z -> interp_ctx ctx g y z }
342 = match ctx with Nil -> () | Cons _ t -> interp_ctx_valid t g end
346 let rec lemma interp_ctx_wr (ctx l:context) (g:equality)
347 ensures { forall y z. interp_ctx ctx g y z -> interp_ctx (ctx ++ l) g y z }
351 | Cons _ t -> interp_ctx_wr t l g end
353 let rec lemma interp_ctx_wl (ctx l: context) (g:equality)
354 ensures { forall y z. interp_ctx ctx g y z -> interp_ctx (l ++ ctx) g y z }
356 = match l with Nil -> () | Cons _ t -> interp_ctx_wl ctx t g end
358 let rec mul_expr (e:expr) (c:coeff) : expr
359 ensures { forall y z. interp result y z
360 = C.(*) (C.interp c z) (interp e y z) }
361 ensures { valid_expr e -> valid_expr result }
363 raises { C.Unknown -> true }
364 = if C.eq c C.czero then Cst C.czero
366 | Cst c1 -> Cst (C.mul c c1)
367 | Term c1 v -> Term (C.mul c c1) v
368 | Add e1 e2 -> Add (mul_expr e1 c) (mul_expr e2 c)
371 let rec add_expr (e1 e2: expr) : expr
372 ensures { forall y z. interp result y z
373 = C.(+) (interp e1 y z) (interp e2 y z) }
375 raises { C.Unknown -> true }
378 ensures { forall y z. interp result y z = interp (Term c i) y z }
379 = if C.eq C.czero c then Cst C.czero else Term c i in
380 let rec add_atom (e a:expr) : (expr, bool)
382 returns { r,_ -> forall y z. interp r y z
383 = C.(+) (interp e y z) (interp a y z) }
385 raises { C.Unknown -> true }
387 | Term ce ie, Term ca ia ->
388 if ie = ia then (term_or_cst (C.add ce ca) ie, True)
389 else if C.eq ce C.czero then (term_or_cst ca ia, True)
390 else if C.eq ca C.czero then (e,True)
391 else (Add e a, False)
392 | Cst ce, Cst ca -> Cst (C.add ce ca), True
393 | Cst ce, Term ca _ ->
394 if C.eq ca C.czero then (e, True)
395 else if C.eq ce C.czero then (a, True)
396 else (Add e a, False)
397 | Term ce _, Cst ca ->
398 if C.eq ce C.czero then (a, True)
399 else if C.eq ca C.czero then (e, True)
400 else (Add e a, False)
402 let r, b = add_atom e1 a in
409 assert { forall y z. C.(+) (interp e1 y z) (interp a y z) = C.azero };
412 | _ -> Add r e2, True
415 let r,b = add_atom e2 a in
420 assert { forall y z. C.(+) (interp e2 y z) (interp a y z) = C.azero };
425 | _, Add _ _ -> absurd
429 | Add e1' e2' -> add_expr (add_expr e1 e1') e2'
430 | _ -> let r,_= add_atom e1 e2 in r
433 let mul_eq (eq:equality) (c:coeff)
434 ensures { forall y z. interp_eq eq y z -> interp_eq result y z }
435 raises { C.Unknown -> true }
436 = match eq with (e1,e2) -> (mul_expr e1 c, mul_expr e2 c) end
438 let add_eq (eq1 eq2:equality)
439 ensures { forall y z. interp_eq eq1 y z -> interp_eq eq2 y z
440 -> interp_eq result y z }
441 ensures { forall y z ctx. interp_ctx ctx eq1 y z -> interp_ctx ctx eq2 y z
442 -> interp_ctx ctx result y z }
443 raises { C.Unknown -> true }
444 = match eq1, eq2 with ((a1,b1), (a2,b2)) ->
445 let a = add_expr a1 a2 in let b = add_expr b1 b2 in
447 let rec lemma aux (l:context)
448 ensures { forall y z. interp_ctx l eq1 y z -> interp_ctx l eq2 y z
449 -> interp_ctx l r y z }
451 = match l with Nil -> () | Cons _ t -> aux t end in
455 let rec zero_expr (e:expr) : bool
456 ensures { result -> forall y z. interp e y z = C.azero }
458 raises { C.Unknown -> true }
460 let rec all_zero (e:expr) : bool
461 ensures { result -> forall y z. interp e y z = C.azero }
464 | Cst c -> C.eq c C.czero
465 | Term c _ -> C.eq c C.czero
466 | Add e1 e2 -> all_zero e1 && all_zero e2
469 let e' = add_expr (Cst C.czero) e in (* simplifies expr *)
472 let sub_expr (e1 e2:expr)
473 ensures { forall y z. C.(+) (interp result y z) (interp e2 y z)
475 raises { C.Unknown -> true }
476 = let r = add_expr e1 (mul_expr e2 (C.opp C.cone)) in
478 let v1 = interp e1 y z in
479 let v2 = interp e2 y z in
480 let vr = interp r y z in
482 by C.(*) v2 (C.(-_) C.aone) = C.(-_) v2
484 = C.(+) (C.(+) v1 (C.(*) v2 (C.(-_) C.aone))) v2
485 = C.(+) (C.(+) v1 (C.(-_) v2)) v2 = v1 };
488 let rec same_eq (eq1 eq2: equality) : bool
489 ensures { result -> forall y z. interp_eq eq1 y z -> interp_eq eq2 y z }
490 raises { C.Unknown -> true }
491 = let (e1,c1) = norm_eq eq1 in
492 let (e2,c2) = norm_eq eq2 in
493 let e = sub_expr e1 e2 in
494 if zero_expr e && C.eq c1 c2 then true
495 else (print (add_expr (Cst C.czero) e); print c1; print c2; false)
499 let rec norm_context (l:context) : context
500 ensures { forall g y z. interp_ctx result g y z -> interp_ctx l g y z }
501 raises { C.Unknown -> true }
506 let ex, c = norm_eq h in
507 Cons (ex, Cst c) (norm_context t)
510 let rec print_lc ctx v : unit variant { ctx }
513 | Cons l t, Cons v t2 ->
514 (if C.eq C.czero v then ()
515 else (print l; print v));
520 let check_combination (ctx:context) (g:equality) (v:list coeff) : bool
521 ensures { result = true -> forall y z. interp_ctx ctx g y z}
522 raises { C.Unknown -> true }
524 (*let ctx = norm_context ctx in
525 let (g,c) = norm_eq g in*)
526 (* normalize before for fewer Unknown exceptions in computations ? *)
527 let rec aux (l:context) (ghost acc: context) (s:equality) (v:list coeff) : option equality
528 requires { forall y z. interp_ctx acc s y z }
529 requires { ctx = acc ++ l }
530 returns { Some r -> forall y z. interp_ctx ctx r y z | None -> true }
531 raises { C.Unknown -> true }
535 | Cons eq te, Cons c tc ->
536 let ghost nacc = acc ++ (Cons eq Nil) in
537 if C.eq c C.czero then aux te nacc s tc
539 let ns = (add_eq s (mul_eq eq c)) in
540 interp_ctx_wr ctx (Cons eq Nil) s;
541 interp_ctx_wl ctx (Cons eq Nil) eq;
542 assert { forall y z. interp_ctx nacc ns y z
543 by interp_ctx nacc s y z /\ interp_ctx nacc eq y z };
544 aux te nacc ns tc end
548 match aux ctx Nil (Cst C.czero, Cst C.czero) v with
549 | Some sum -> if same_eq sum g then true else (print_lc ctx v; false)
553 let transpose (m:matrix coeff) : matrix coeff
554 ensures { result.rows = m.columns /\ result.columns = m.rows }
556 let r = Matrix63.make m.columns m.rows C.czero in
557 for i = 0 to m.rows - 1 do
558 for j = 0 to m.columns - 1 do
559 set r j i (get m i j)
564 let swap_rows (m:matrix coeff) (i1 i2: int63) : unit
565 requires { 0 <= i1 < m.rows /\ 0 <= i2 < m.rows }
566 = for j = 0 to m.columns - 1 do
567 let c = get m i1 j in
568 set m i1 j (get m i2 j);
572 let mul_row (m:matrix coeff) (i: int63) (c: coeff) : unit
573 requires { 0 <= i < m.rows }
574 requires { not (C.eq c C.czero) }
575 raises { C.Unknown -> true }
576 = if C.eq c C.cone then () else
577 for j = 0 to m.columns - 1 do
578 set m i j (C.mul c (get m i j))
581 let addmul_row (m:matrix coeff) (src dst: int63) (c: coeff) : unit
582 requires { 0 <= src < m.rows /\ 0 <= dst < m.rows }
583 raises { C.Unknown -> true }
584 = if C.eq c C.czero then () else
585 for j = 0 to m.columns - 1 do
586 set m dst j (C.add (get m dst j) (C.mul c (get m src j)))
591 let gauss_jordan (a: matrix coeff) : option (array coeff)
592 (*AX=B, a=(A|B), result=X*)
593 returns { Some r -> Array63.length r = a.columns | None -> true }
594 requires { 1 <= a.rows /\ 1 <= a.columns }
595 raises { C.Unknown -> true }
599 (* print n; print m; *)
600 let rec find_nonz (i j:int63)
601 requires { 0 <= i <= n }
602 requires { 0 <= j < m }
604 ensures { i <= result <= n }
605 ensures { result < n -> not (C.eq (a.elts result j) C.czero) }
608 if C.eq (get a i j) C.czero
609 then find_nonz (i+1) j
611 let pivots = Array63.make n 0 in
614 invariant { -1 <= !r < n }
615 invariant { forall i. 0 <= i <= !r -> 0 <= pivots[i] }
616 invariant { forall i1 i2: int. 0 <= i1 < i2 <= !r -> pivots[i1] < pivots[i2] }
617 invariant { !r >= 0 -> pivots[!r] < j }
619 let k = find_nonz (!r+1) j in
624 mul_row a k (C.inv(get a k j));
625 if k <> !r then swap_rows a k !r;
628 then addmul_row a !r i (C.opp(get a i j))
631 assert { forall i1 i2: int. 0 <= i1 < i2 <= !r -> pivots[i1] < pivots[i2]
632 by pivots[i1] = pivots[i1] at Start
634 ((i2 < !r so pivots[i2] = pivots[i2] at Start)
635 \/ (i2 = !r so pivots[i1] < j(* = pivots[i2])*))) };
637 if !r < 0 then None (* matrix is all zeroes *)
639 let v = Array63.make m(*(m-1)*) C.czero in
641 v[pivots[i]] <- get a i (m-1)
643 Some v (*pivots[!r] < m-1*) (*pivot on last column, no solution*)
646 let rec function to_list (a: array 'a) (l u: int63) : list 'a
647 requires { l >= 0 /\ u <= Array63.length a }
649 = if u <= l then Nil else Cons a[l] (to_list a (l+1) u)
653 let linear_decision (l: context) (g: equality) : bool
654 requires { valid_ctx l }
655 requires { valid_eq g }
656 requires { length l < 100000 } (* integer overflows *)
657 ensures { forall y z. result -> interp_ctx l g y z }
658 raises { C.Unknown -> true | Failure -> true }
660 let nv = (max (max_var_e g) (max_var_ctx l)) in
661 begin ensures { nv < 100000 }
662 if nv >= 100000 then raise Failure
664 let nv = Int63.of_int nv in
665 let ll = Int63.of_int (length l) in
666 let a = Matrix63.make ll (nv+1) C.czero in
667 let b = Array63.make ll C.czero in (* ax = b *)
668 let v = Array63.make (nv+1) C.czero in (* goal *)
669 let rec fill_expr (ex: expr) (i:int63): unit
671 raises { C.Unknown -> true }
672 requires { 0 <= i < length l }
673 requires { expr_bound ex nv }
674 raises { Failure -> true }
676 | Cst c -> if C.eq c C.czero then () else raise Failure
678 let j = Int63.of_int j in
679 set a i j (C.add (get a i j) c)
680 | Add e1 e2 -> fill_expr e1 i; fill_expr e2 i
682 let rec fill_ctx (ctx:context) (i:int63) : unit
683 requires { ctx_bound ctx nv }
684 variant { length l - i }
685 requires { length l - i = length ctx }
686 requires { 0 <= i <= length l }
687 raises { Failure -> true }
691 assert { i < length l };
693 let ex, c = norm_eq e in
694 if (not (C.eq c C.czero)) then b[i] <- C.add b[i] c;
696 with C.Unknown -> () (* some equalities are in the context but cannot be normalized, typically they are useless, ignore them *)
700 let rec fill_goal (ex:expr) : unit
701 requires { expr_bound ex nv }
703 raises { C.Unknown -> true }
704 raises { Failure -> true }
706 | Cst c -> if C.eq c C.czero then () else raise Failure
708 let j = Int63.of_int j in
710 | Add e1 e2 -> fill_goal e1; fill_goal e2
713 let (ex, d) = norm_eq g in
715 let ab = m_append a b in
716 let cd = v_append v d in
717 let ab' = transpose ab in
718 match gauss_jordan (m_append ab' cd) with
720 check_combination l g (to_list r 0 ll)
724 type expr' = | Sum expr' expr' | ProdL expr' cprod | ProdR cprod expr' | Diff expr' expr'
725 | Var int | Coeff coeff
727 with cprod = | C coeff | Times cprod cprod
729 function interp_c (e:cprod) (y:vars) (z:C.cvars) : C.a
731 | C c -> C.interp c z
732 | Times e1 e2 -> C.(*) (interp_c e1 y z) (interp_c e2 y z)
735 function interp' (e:expr') (y:vars) (z:C.cvars) : C.a
737 | Sum e1 e2 -> C.(+) (interp' e1 y z) (interp' e2 y z)
738 | ProdL e c -> C.(*) (interp' e y z) (interp_c c y z)
739 | ProdR c e -> C.(*) (interp_c c y z) (interp' e y z)
740 | Diff e1 e2 -> C.(-) (interp' e1 y z) (interp' e2 y z)
742 | Coeff c -> C.interp c z
745 (*exception NonLinear*)
747 type equality' = (expr', expr')
748 type context' = list equality'
750 function interp_eq' (g:equality') (y:vars) (z:C.cvars) : bool
751 = match g with (g1, g2) -> interp' g1 y z = interp' g2 y z end
753 function interp_ctx' (l: context') (g: equality') (y: vars) (z:C.cvars) : bool
755 | Nil -> interp_eq' g y z
756 | Cons h t -> (interp_eq' h y z) -> (interp_ctx' t g y z)
759 let rec predicate valid_expr' (e:expr')
763 | Sum e1 e2 | Diff e1 e2 -> valid_expr' e1 && valid_expr' e2
765 | ProdL e _ | ProdR _ e -> valid_expr' e
768 let predicate valid_eq' (eq:equality')
769 = match eq with (e1,e2) -> valid_expr' e1 && valid_expr' e2 end
771 let rec predicate valid_ctx' (ctx:context')
772 = match ctx with Nil -> true | Cons eq t -> valid_eq' eq && valid_ctx' t end
774 let rec simp (e:expr') : expr
775 ensures { forall y z. interp result y z = interp' e y z }
776 ensures { valid_expr' e -> valid_expr result }
777 raises { C.Unknown -> true }
780 let rec simp_c (e:cprod) : coeff
781 ensures { forall y z. C.interp result z = interp_c e y z }
783 raises { C.Unknown -> true }
787 | Times c1 c2 -> C.mul (simp_c c1) (simp_c c2)
791 | Sum e1 e2 -> Add (simp e1) (simp e2)
792 | Diff e1 e2 -> Add (simp e1) (opp_expr (simp e2))
793 | Var n -> Term C.cone n
795 | ProdL e c | ProdR c e ->
796 mul_expr (simp e) (simp_c c)
799 let simp_eq (eq:equality') : equality
800 ensures { forall y z. interp_eq result y z = interp_eq' eq y z }
801 ensures { valid_eq' eq -> valid_eq result }
802 raises { (*NonLinear -> true | *)C.Unknown -> true }
803 = match eq with (g1, g2) -> (simp g1, simp g2) end
805 let rec simp_ctx (ctx: context') (g:equality') : (context, equality)
806 returns { (rc, rg) ->
807 (valid_ctx' ctx -> valid_eq' g -> valid_ctx rc /\ valid_eq rg) /\
808 length rc = length ctx /\
809 forall y z. interp_ctx rc rg y z = interp_ctx' ctx g y z }
810 raises { (*NonLinear -> true | *) C.Unknown -> true }
813 | Nil -> Nil, simp_eq g
814 | Cons eq t -> let rt, rg = simp_ctx t g in
815 Cons (simp_eq eq) rt, rg
818 let decision (l:context') (g:equality')
819 requires { valid_ctx' l }
820 requires { valid_eq' g }
821 requires { length l < 100000 }
822 ensures { forall y z. result -> interp_ctx' l g y z }
823 raises { (* NonLinear -> true | *) C.Unknown -> true | Failure -> true }
824 = let sl, sg = simp_ctx l g in
825 linear_decision sl sg
829 module RationalCoeffs
837 type rvars = int -> real
841 let constant rzero = (0,1)
842 let constant rone = (1,1)
844 function rinterp (t:t) (v:rvars) : real
846 | (n,d) -> from_int n /. from_int d
849 let lemma prod_compat_eq (a b c:real)
850 requires { c <> 0.0 }
851 requires { a *. c = b *. c }
855 let lemma cross_d (n1 d1 n2 d2:int)
856 requires { d1 <> 0 /\ d2 <> 0 }
857 requires { n1 * d2 = n2 * d1 }
858 ensures { forall v. rinterp (n1,d1) v = rinterp (n2,d2) v }
859 = let d = from_int (d1 * d2) in
860 assert { forall v. rinterp (n1, d1) v = rinterp (n2, d2) v
861 by rinterp (n1, d1) v *. d = rinterp (n2,d2) v *. d }
863 let lemma cross_ind (n1 d1 n2 d2:int)
864 requires { d1 <> 0 /\ d2 <> 0 }
865 requires { forall v. rinterp (n1,d1) v = rinterp (n2,d2) v }
866 ensures { n1 * d2 = n2 * d1 }
867 = assert { from_int d1 <> 0.0 /\ from_int d2 <> 0.0 };
868 assert { from_int n1 /. from_int d1 = from_int n2 /. from_int d2 };
869 assert { from_int n1 *. from_int d2 = from_int n2 *. from_int d1
870 by from_int n1 *. from_int d2
871 = (from_int n1 /. from_int d1) *. from_int d1 *. from_int d2
872 = (from_int n2 /. from_int d2) *. from_int d1 *. from_int d2
873 = from_int n2 *. from_int d1 };
874 assert { from_int (n1*d2) = from_int (n2 * d1) }
877 lemma cross: forall n1 d1 n2 d2: int. d1 <> 0 -> d2 <> 0 ->
878 n1 * d2 = n2 * d1 <->
879 forall v. rinterp (n1,d1) v = rinterp (n2,d2) v
881 use int.ComputerDivision
885 let gcd (x:int) (y:int)
886 requires { x > 0 /\ y > 0 }
887 ensures { result = gcd x y }
888 ensures { result > 0 }
891 let x = ref x in let y = ref y in
894 invariant { !x >= 0 /\ !y >= 0 }
895 invariant { gcd !x !y = gcd (!x at Pre) (!y at Pre) }
897 invariant { ox > 0 -> !x > 0 }
898 let r = mod !x !y in let ghost q = div !x !y in
899 assert { r = !x - q * !y };
905 ensures { forall v:rvars. rinterp result v = rinterp t v }
909 else if n = 0 then rzero
911 let g = gcd (abs n) (abs d) in
912 let n', d' = (div n g, div d g) in
913 assert { n = g * n' /\ d = g * d' };
914 assert { n' * d = n * d' };
919 ensures { forall y. rinterp result y = rinterp a y +. rinterp b y }
920 raises { QError -> true }
922 | (n1,d1), (n2,d2) ->
923 if d1 = 0 || d2 = 0 then raise QError
925 let r = (n1*d2 + n2*d1, d1*d2) in
926 let ghost d = from_int d1 *. from_int d2 in
928 rinterp a y +. rinterp b y = rinterp r y
929 by rinterp a y *. d = from_int n1 *. from_int d2
930 so rinterp b y *. d = from_int n2 *. from_int d1
931 so (rinterp a y +. rinterp b y) *. d
932 = from_int (n1*d2 + n2*d1)
933 = rinterp r y *. d };
938 ensures { forall y. rinterp result y = rinterp a y *. rinterp b y }
939 raises { QError -> true }
941 | (n1,d1), (n2, d2) ->
942 if d1 = 0 || d2 = 0 then raise QError
944 let r = (n1*n2, d1*d2) in
945 assert { forall y. rinterp r y = rinterp a y *. rinterp b y
946 by rinterp r y = from_int (n1*n2) /. from_int(d1*d2)
947 = (from_int n1 *. from_int n2) /. (from_int d1 *. from_int d2)
948 = (from_int n1 /. from_int d1) *. (from_int n2 /. from_int d2)
949 = rinterp a y *. rinterp b y };
955 ensures { forall y. rinterp result y = -. rinterp a y }
960 let predicate req (a b:t)
961 ensures { result -> forall y. rinterp a y = rinterp b y }
963 | (n1,d1), (n2,d2) -> n1 = n2 && d1 = d2 || (d1 <> 0 && d2 <> 0 && n1 * d2 = n2 * d1)
967 requires { not req a rzero }
968 ensures { not req result rzero }
969 ensures { forall y. rinterp result y *. rinterp a y = 1.0 }
970 raises { QError -> true }
972 | (n,d) -> if n = 0 || d = 0 then raise QError else (d,n)
976 ensures { result <-> req a rzero }
978 | (n,d) -> n = 0 && d <> 0
983 module LinearDecisionRational
989 clone export LinearEquationsDecision with type C.a = real, function C.(+) = (+.), function C.( * ) = ( *. ), function C.(-_) = (-._), function C.(-) = (-.), type coeff = t, type C.cvars=int -> real, function C.interp=rinterp, exception C.Unknown = QError, constant C.azero = Real.zero, constant C.aone = Real.one, predicate C.ale = (<=.), val C.czero=rzero, val C.cone=rone, lemma C.sub_def, lemma C.zero_def, lemma C.one_def, val C.add=radd, val C.mul=rmul, val C.opp=ropp, val C.eq=req, val C.inv=rinv, goal .
993 module LinearDecisionInt
997 type t' = IC int | Error
999 function interp_id (t:t') (v:int -> int) : int
1002 | Error -> 0 (* never created *)
1005 let constant izero = IC 0
1007 let constant ione = IC 1
1009 let predicate ieq (a b:t') = false
1013 let iadd (a b:t') : t'
1014 ensures { forall z. interp_id result z = interp_id a z + interp_id b z }
1015 raises { NError -> true }
1018 let imul (a b:t') : t'
1019 ensures { forall z. interp_id result z = interp_id a z * interp_id b z }
1020 raises { NError -> true }
1023 let iopp (a:t') : t'
1024 ensures { forall z. interp_id result z = - interp_id a z }
1025 raises { NError -> true }
1028 let iinv (t:t') : t'
1029 (*ensures { forall v: int -> int. id result v * id t v = one }*)
1030 ensures { not (ieq result izero) }
1031 raises { NError -> true }
1034 clone export LinearEquationsDecision with type C.a = int, function C.(+)=(+), function C.(*) = (*), function C.(-_) = (-_), function C.(-) = (-), type coeff = t', type C.cvars = int->int,function C.interp = interp_id, constant C.azero = zero, constant C.aone = one, predicate C.ale= (<=), val C.czero = izero, val C.cone = ione, lemma C.sub_def, lemma C.zero_def, lemma C.one_def, val C.add = iadd, val C.mul = imul, val C.opp = iopp, val C.eq = ieq, val C.inv = iinv, goal .
1039 use LinearDecisionRational as R
1042 let ghost function m_y (y:int -> int): (int -> real)
1043 ensures { forall i. result i = from_int (y i) }
1044 = fun i -> from_int (y i)
1046 let m (t:t') : (int, int)
1047 ensures { forall z. rinterp result (m_y z) = from_int (interp_id t z) }
1048 raises { NError -> true }
1054 let rec m_cprod (e:cprod) : R.cprod
1055 ensures { forall y z. R.interp_c result (m_y y) (m_y z)
1056 = from_int (interp_c e y z) }
1057 raises { NError -> true }
1061 | Times c1 c2 -> R.Times (m_cprod c1) (m_cprod c2)
1064 let rec m_expr (e:expr') : R.expr'
1065 ensures { forall y z. R.interp' result (m_y y) (m_y z)
1066 = from_int (interp' e y z) }
1067 ensures { valid_expr' e -> R.valid_expr' result }
1068 raises { NError -> true }
1072 | Coeff c -> R.Coeff (m c)
1073 | Sum e1 e2 -> R.Sum (m_expr e1) (m_expr e2)
1074 | Diff e1 e2 -> R.Diff (m_expr e1) (m_expr e2)
1075 | ProdL e c -> R.ProdL (m_expr e) (m_cprod c)
1076 | ProdR c e -> R.ProdR (m_cprod c) (m_expr e)
1082 let m_eq (eq:equality') : R.equality'
1083 ensures { forall y z. R.interp_eq' result (m_y y) (m_y z)
1084 <-> interp_eq' eq y z }
1085 ensures { valid_eq' eq -> R.valid_eq' result }
1086 raises { NError -> true }
1087 = match eq with (e1,e2) -> (m_expr e1, m_expr e2) end
1089 let rec m_ctx (ctx:context') (g:equality') : (R.context', R.equality')
1090 returns { c',g' -> forall y z. R.interp_ctx' c' g' (m_y y) (m_y z) <->
1091 interp_ctx' ctx g y z }
1092 returns { c', _ -> valid_ctx' ctx -> R.valid_ctx' c' }
1093 returns { c', _ -> length c' = length ctx }
1094 returns { _, g' -> valid_eq' g -> R.valid_eq' g' }
1095 raises { NError -> true }
1098 | Nil -> Nil, m_eq g
1100 let c',g' = m_ctx t g in
1101 (Cons (m_eq h) c',g')
1104 let int_decision (l: context') (g: equality') : bool
1105 requires { valid_ctx' l }
1106 requires { valid_eq' g }
1107 requires { length l < 100000 }
1108 ensures { forall y z. result -> interp_ctx' l g y z }
1109 raises { R.Failure -> true | QError -> true | NError -> true }
1110 = let l',g' = m_ctx l g in
1119 use LinearDecisionRational
1124 meta "compute_max_steps" 0x10000
1125 meta coercion function from_int
1127 goal g: forall x y: real.
1128 (from_int 3 /. from_int 1) *. x +. (from_int 2/. from_int 1) *. y = (from_int 21/. from_int 1) ->
1129 (from_int 7 /. from_int 1) *. x +. (from_int 4/. from_int 1) *. y = (from_int 47/. from_int 1) ->
1130 x = (from_int 5 /. from_int 1)
1135 use LinearDecisionInt
1138 meta "compute_max_steps" 0x10000
1140 goal g: forall x y:int.
1141 3 * x + 2 * y = 21 ->
1142 7 * x + 4 * y = 47 ->
1149 use mach.int.UInt64 as M
1153 use RationalCoeffs as Q
1158 type evars = int -> int
1161 type exp = Lit int | Var int | Plus exp exp | Minus exp | Sub exp exp
1164 let constant mzero = (Q.rzero, Lit 0)
1165 let constant mone = (Q.rone, Lit 0)
1167 constant rradix: real = from_int (M.radix)
1169 function qinterp (q:Q.t) : real
1170 = match q with (n,d) -> from_int n /. from_int d end
1172 lemma qinterp_def: forall q v. qinterp q = Q.rinterp q v
1174 function interp_exp (e:exp) (y:evars) : int
1178 | Plus e1 e2 -> interp_exp e1 y + interp_exp e2 y
1179 | Sub e1 e2 -> interp_exp e1 y - interp_exp e2 y
1180 | Minus e' -> - (interp_exp e' y)
1183 function minterp (t:t) (y:evars) : real
1186 qinterp q *. pow rradix (from_int (interp_exp e y))
1191 let rec opp_exp (e:exp)
1192 ensures { forall y. interp_exp result y = - interp_exp e y }
1197 | Plus e1 e2 -> Plus (opp_exp e1) (opp_exp e2)
1198 | Sub e1 e2 -> Sub e2 e1
1202 let rec add_sub_exp (e1 e2:exp) (s:bool) : exp
1205 then interp_exp result y = interp_exp e1 y + interp_exp e2 y
1206 else interp_exp result y = interp_exp e1 y - interp_exp e2 y }
1207 raises { MPError -> true }
1210 let rec add_atom (e a:exp) (s:bool) : (exp, bool)
1211 returns { r, _ -> forall y.
1212 if s then interp_exp r y = interp_exp e y + interp_exp a y
1213 else interp_exp r y = interp_exp e y - interp_exp a y }
1214 raises { MPError -> true }
1217 | Lit n1, Lit n2 -> (if s then Lit (n1+n2) else Lit (n1-n2)), True
1219 -> if n = 0 then (if s then Var i else Minus (Var i)), True
1220 else (if s then Plus e a else Sub e a), False
1222 -> if n = 0 then Var i, true
1223 else (if s then Plus e a else Sub e a), False
1224 | Lit n, Minus e' ->
1225 if n = 0 then (if s then Minus e' else e'), True
1226 else (if s then Plus e a else Sub e a), False
1227 | Minus e', Lit n ->
1228 if n = 0 then Minus e', True
1229 else (if s then Plus e a else Sub e a), False
1230 | Var i, Minus (Var j) | Minus (Var j), Var i ->
1231 if s && (i = j) then (Lit 0, true)
1232 else (if s then Plus e a else Sub e a), False
1233 | Var i, Var j -> if s then Plus e a, False
1235 if i = j then Lit 0, True
1237 | Minus (Var i), Minus (Var j) ->
1238 if (not s) && (i=j) then Lit 0, true
1239 else (if s then Plus e a else Sub e a), False
1240 | Minus _, Minus _ -> (if s then Plus e a else Sub e a), False
1242 let r, b = add_atom e1 a s in
1245 | Lit n -> if n = 0 then e2, True else Plus r e2, True
1246 | _ -> Plus r e2, True
1248 else let r, b = add_atom e2 a s in Plus e1 r, b
1250 let r, b = add_atom e1 a s in
1253 | Lit n -> if n = 0 then opp_exp e2, True else Sub r e2, True
1254 | _ -> Sub r e2, True
1256 else let r, b = add_atom e2 a (not s) in
1257 if b then Sub e1 r, True
1258 else if s then Sub (Plus e1 a) e2, False
1259 else Sub e1 (Plus e2 a), False
1260 | _ -> raise MPError
1265 let r = add_sub_exp e1 e1' s in
1268 then (if s then e2' else opp_exp e2')
1269 else add_sub_exp r e2' s
1270 | _ -> add_sub_exp r e2' s
1273 let r = add_sub_exp e1 e1' s in
1276 then (if s then opp_exp e2' else e2')
1277 else add_sub_exp r e2' (not s)
1278 | _ -> add_sub_exp r e2' (not s)
1280 | _ -> let r, _ = add_atom e1 e2 s in r
1283 let add_exp (e1 e2:exp) : exp
1284 ensures { forall y. interp_exp result y = interp_exp e1 y + interp_exp e2 y }
1285 raises { MPError -> True }
1286 = add_sub_exp e1 e2 True
1289 let rec zero_exp (e:exp) : bool
1290 ensures { result -> forall y. interp_exp e y = 0 }
1292 raises { MPError -> true }
1294 let rec all_zero (e:exp) : bool
1295 ensures { result -> forall y. interp_exp e y = 0 }
1300 | Minus e -> all_zero e
1301 | Plus e1 e2 -> all_zero e1 && all_zero e2
1302 | Sub e1 e2 -> all_zero e1 && all_zero e2
1305 let e' = add_exp (Lit 0) e in (* simplifies exp *)
1308 let rec same_exp (e1 e2: exp)
1309 ensures { result -> forall y. interp_exp e1 y = interp_exp e2 y }
1311 raises { MPError -> true }
1313 | Lit n1, Lit n2 -> n1 = n2
1314 | Var v1, Var v2 -> v1 = v2
1315 | Minus e1', Minus e2' -> same_exp e1' e2'
1316 | _ -> zero_exp (add_exp e1 (opp_exp e2))
1320 ensures { forall y. minterp result y = minterp a y +. minterp b y }
1321 raises { MPError -> true }
1322 raises { Q.QError -> true }
1324 | (q1, e1), (q2, e2) ->
1325 if Q.is_zero q1 then b
1326 else if Q.is_zero q2 then a
1327 else if same_exp e1 e2
1329 let q = Q.radd q1 q2 in
1330 assert { forall y. minterp (q, e1) y = minterp a y +. minterp b y
1331 by let p = pow rradix (from_int (interp_exp e1 y)) in
1332 minterp (q, e1) y = (qinterp q) *. p
1333 = (qinterp q1 +. qinterp q2) *. p
1334 = qinterp q1 *. p +. qinterp q2 *. p
1335 = minterp a y +. minterp b y };
1341 ensures { forall y. minterp result y = minterp a y *. minterp b y }
1342 raises { Q.QError -> true }
1343 raises { MPError -> true }
1345 | (q1,e1), (q2,e2) ->
1346 let q = Q.rmul q1 q2 in
1347 if Q.is_zero q then mzero
1349 let e = add_exp e1 e2 in
1350 assert { forall y. minterp (q,e) y = minterp a y *. minterp b y
1351 by let p1 = pow rradix (from_int (interp_exp e1 y)) in
1352 let p2 = pow rradix (from_int (interp_exp e2 y)) in
1353 let p = pow rradix (from_int (interp_exp e y)) in
1354 interp_exp e y = interp_exp e1 y + interp_exp e2 y
1356 so minterp (q,e) y = qinterp q *. p
1357 = (qinterp q1 *. qinterp q2) *. p
1358 = (qinterp q1 *. qinterp q2) *. p1 *. p2
1359 = minterp a y *. minterp b y };
1365 ensures { forall y. minterp result y = -. minterp a y }
1366 = match a with (q,e) -> (Q.ropp q, e) end
1368 let rec predicate pure_same_exp (e1 e2: exp)
1369 ensures { result -> forall y. interp_exp e1 y = interp_exp e2 y }
1372 | Lit n1, Lit n2 -> n1 = n2
1373 | Var v1, Var v2 -> v1 = v2
1374 | Minus e1', Minus e2' -> pure_same_exp e1' e2'
1375 | Plus a1 a2, Plus b1 b2 ->
1376 (pure_same_exp a1 b1 && pure_same_exp a2 b2) ||
1377 (pure_same_exp a1 b2 && pure_same_exp a2 b1)
1381 let predicate meq (a b:t)
1382 ensures { result -> forall y. minterp a y = minterp b y }
1384 | (q1,e1), (q2,e2) -> (Q.req q1 q2 && pure_same_exp e1 e2) || (Q.is_zero q1 && Q.is_zero q2)
1388 requires { not meq a mzero }
1389 ensures { not meq result mzero }
1390 (* ensures { forall y. minterp result y *. minterp a y = 1.0 } no need to prove this*)
1391 raises { Q.QError -> true }
1393 | (q,e) -> (Q.rinv q, opp_exp e)
1398 module LinearDecisionRationalMP
1405 clone export LinearEquationsDecision with type C.a = real, function C.(+) = (+.), function C.(*) = ( *.), function C.(-_) = (-._), function C.(-) = (-.), type coeff = t, type C.cvars=evars, function C.interp=minterp, exception C.Unknown = MPError, constant C.azero = Real.zero, constant C.aone = Real.one, predicate C.ale = (<=.), val C.czero=mzero, val C.cone=mone, lemma C.sub_def, lemma C.zero_def, lemma C.one_def, val C.add=madd, val C.mul=mmul, val C.opp=mopp, val C.eq=meq, val C.inv=minv, goal .
1408 module LinearDecisionIntMP
1414 type t = | I int | E exp | R
1416 let constant mpzero: t = I 0
1417 let constant mpone: t = I 1
1419 function mpinterp (t:t) (y:evars) : int
1422 | E e -> power M.radix (interp_exp e y)
1426 (* TODO restructure stuff so that expr, eq, ctx, valid_ can be imported without having to implement these *)
1428 let mpadd (a b:t) : t
1429 ensures { forall y. mpinterp result y = mpinterp a y + mpinterp b y }
1430 raises { MPError -> true }
1433 let mpmul (a b:t) : t
1434 ensures { forall y. mpinterp result y = mpinterp a y * mpinterp b y }
1435 raises { MPError -> true }
1439 ensures { forall y. mpinterp result y = - mpinterp a y }
1440 raises { MPError -> true }
1443 let predicate mpeq (a b:t)
1444 ensures { result -> forall y. mpinterp a y = mpinterp b y }
1445 = false (*match a, b with
1446 (n1, e1), (n2, e2) -> n1=n2 && (n1 = 0 || same_exp e1 e2)
1450 ensures { not mpeq result mpzero }
1451 raises { MPError -> true }
1455 clone export LinearEquationsDecision with type C.a = int, function C.(+) = (+), function C.(*) = (*), function C.(-_) = (-_), function C.(-) = (-), type coeff = t, type C.cvars = int->int, function C.interp = mpinterp, constant C.azero = zero, constant C.aone = one, val C.czero = mpzero, val C.cone = mpone, predicate C.ale = (<=), lemma C.sub_def, lemma C.zero_def, lemma C.one_def, val C.add = mpadd, val C.mul = mpmul, val C.opp = mpopp, val C.eq = mpeq, val C.inv = mpinv, goal .
1457 use LinearDecisionRationalMP as R
1465 predicate pos_exp (t:t) (y:evars)
1467 | E e -> 0 <= interp_exp e y
1471 predicate pos_cprod (e:cprod) (y:evars)
1473 | C c -> pos_exp c y
1474 | Times c1 c2 -> pos_cprod c1 y && pos_cprod c2 y
1477 predicate pos_expr' (e:expr') (y:evars)
1479 | Coeff c -> pos_exp c y
1481 | Sum e1 e2 | Diff e1 e2
1482 -> pos_expr' e1 y /\ pos_expr' e2 y
1483 | ProdL e c | ProdR c e -> pos_expr' e y && pos_cprod c y
1486 predicate pos_eq' (eq:equality') (y:evars)
1487 = match eq with (e1, e2) -> pos_expr' e1 y /\ pos_expr' e2 y end
1489 predicate pos_ctx' (l:context') (y:evars)
1490 = match l with Nil -> true | Cons h t -> pos_eq' h y /\ pos_ctx' t y end
1492 let rec function m (t:t) : R.coeff
1493 ensures { forall y. pos_exp t y -> minterp result y
1494 = from_int (mpinterp t y) }
1496 | I n -> ((n,1), Lit 0)
1498 | R -> ((1,1), Lit 1) (* or ((radix, 1), Lit 0) ? *)
1501 let ghost function m_y (y:int->int): (int -> real)
1502 ensures { forall i. result i = from_int (y i) }
1503 = fun i -> from_int (y i)
1505 let rec function m_cprod (e:cprod) : R.cprod
1506 ensures { forall y z. pos_cprod e z -> R.interp_c result (m_y y) z
1507 = from_int (interp_c e y z) }
1510 | Times c1 c2 -> R.Times (m_cprod c1) (m_cprod c2)
1513 let rec function m_expr (e:expr') : R.expr'
1514 ensures { forall y z. pos_expr' e z -> R.interp' result (m_y y) z
1515 = from_int (interp' e y z) }
1516 ensures { valid_expr' e -> R.valid_expr' result}
1519 | Coeff c -> R.Coeff (m c)
1520 | Sum e1 e2 -> R.Sum (m_expr e1) (m_expr e2)
1521 | Diff e1 e2 -> R.Diff (m_expr e1) (m_expr e2)
1522 | ProdL e c -> R.ProdL (m_expr e) (m_cprod c)
1523 | ProdR c e -> R.ProdR (m_cprod c) (m_expr e)
1526 let function m_eq (eq:equality') : R.equality'
1527 ensures { forall y z. pos_eq' eq z -> (R.interp_eq' result (m_y y) z
1528 <-> interp_eq' eq y z) }
1529 ensures { valid_eq' eq -> R.valid_eq' result }
1530 = match eq with (e1,e2) -> (m_expr e1, m_expr e2) end
1534 let rec function m_ctx (ctx:context') : R.context'
1535 ensures { forall y z g. pos_ctx' ctx z -> pos_eq' g z ->
1536 (R.interp_ctx' result (m_eq g) (m_y y) z
1537 <-> interp_ctx' ctx g y z) }
1538 ensures { length result = length ctx }
1539 ensures { valid_ctx' ctx -> R.valid_ctx' result }
1544 let r = Cons (m_eq h) (m_ctx t) in
1548 let mp_decision (l: context') (g: equality') : bool
1549 requires { valid_ctx' l }
1550 requires { valid_eq' g }
1551 requires { length l < 100000 }
1552 ensures { forall y z. result -> pos_ctx' l z -> pos_eq' g z
1553 -> interp_ctx' l g y z }
1554 raises { R.Failure -> true | MPError -> true | Q.QError -> true }
1556 R.decision (m_ctx l) (m_eq g)
1563 use LinearDecisionIntMP
1573 let rec predicate expr_bound' (e:expr') (b:int)
1576 | Sum e1 e2 | Diff e1 e2 -> expr_bound' e1 b && expr_bound' e2 b
1577 | ProdL e _ | ProdR _ e -> expr_bound' e b
1578 | Var n -> 0 <= n <= b
1582 let predicate eq_bound' (eq:equality') (b:int)
1583 = match eq with (e1,e2) -> expr_bound' e1 b && expr_bound' e2 b end
1585 let rec predicate ctx_bound' (ctx:context') (b:int)
1586 = match ctx with Nil -> true | Cons eq t -> eq_bound' eq b && ctx_bound' t b end
1588 let rec lemma expr_bound_w' (e:expr') (b1 b2:int)
1589 requires { b1 <= b2 }
1590 requires { expr_bound' e b1 }
1591 ensures { expr_bound' e b2 }
1594 | Sum e1 e2 | Diff e1 e2 ->
1595 expr_bound_w' e1 b1 b2; expr_bound_w' e2 b1 b2
1596 | ProdL e _ | ProdR _ e -> expr_bound_w' e b1 b2
1600 lemma eq_bound_w': forall e:equality', b1 b2:int. eq_bound' e b1 -> b1 <= b2 -> eq_bound' e b2
1602 let rec lemma ctx_bound_w' (l:context') (b1 b2:int)
1603 requires { ctx_bound' l b1 }
1604 requires { b1 <= b2 }
1605 ensures { ctx_bound' l b2 }
1607 = match l with Nil -> () | Cons _ t -> ctx_bound_w' t b1 b2 end
1609 let rec function max_var' (e:expr') : int
1611 requires { valid_expr' e }
1612 ensures { 0 <= result }
1613 ensures { expr_bound' e result }
1617 | Sum e1 e2 | Diff e1 e2 -> max (max_var' e1) (max_var' e2)
1618 | ProdL e _ | ProdR _ e -> max_var' e
1621 let function max_var_e' (e:equality') : int
1622 requires { valid_eq' e }
1623 ensures { 0 <= result }
1624 ensures { eq_bound' e result }
1625 = match e with (e1,e2) -> max (max_var' e1) (max_var' e2) end
1627 let rec function max_var_ctx' (l:context') : int
1629 requires { valid_ctx' l }
1630 ensures { 0 <= result }
1631 ensures { ctx_bound' l result }
1634 | Cons e t -> max (max_var_e' e) (max_var_ctx' t)
1637 let rec lemma interp_ctx_valid' (ctx:context') (g:equality')
1638 ensures { forall y z. interp_eq' g y z -> interp_ctx' ctx g y z }
1640 = match ctx with Nil -> () | Cons _ t -> interp_ctx_valid' t g end
1642 let rec lemma interp_ctx_wr' (ctx l:context') (g:equality')
1643 ensures { forall y z. interp_ctx' ctx g y z -> interp_ctx' (ctx ++ l) g y z }
1647 | Cons _ t -> interp_ctx_wr' t l g end
1649 let rec lemma interp_ctx_wl' (ctx l: context') (g:equality')
1650 ensures { forall y z. interp_ctx' ctx g y z -> interp_ctx' (l ++ ctx) g y z }
1652 = match l with Nil -> () | Cons _ t -> interp_ctx_wl' ctx t g end
1655 let lemma interp_ctx_cons' (e:equality') (l:context') (g:equality')
1656 requires { forall y z. (interp_eq' e y z -> interp_ctx' l g y z) }
1657 ensures { forall y z. interp_ctx' (Cons e l) g y z }
1660 predicate ctx_impl_ctx' (c1 c2: context')
1663 | Cons eq t -> ctx_impl_ctx' c1 t /\ forall y z. y=z -> interp_ctx' c1 eq y z
1666 predicate ctx_holds' (c: context') (y z:vars)
1669 | Cons h t -> interp_eq' h y z /\ ctx_holds' t y z
1672 let rec lemma holds_interp_ctx' (l:context') (g:equality') (y z:vars)
1673 requires { ctx_holds' l y z -> interp_eq' g y z }
1674 ensures { interp_ctx' l g y z }
1679 if interp_eq' h y z then holds_interp_ctx' t g y z
1682 let rec lemma interp_holds' (l:context') (g:equality') (y z:vars)
1683 requires { interp_ctx' l g y z }
1684 requires { ctx_holds' l y z }
1685 ensures { interp_eq' g y z }
1689 | Cons _ t -> interp_holds' t g y z
1692 let rec lemma impl_holds' (c1 c2: context') (y z: vars)
1693 requires { ctx_impl_ctx' c1 c2 }
1695 requires { ctx_holds' c1 y z }
1696 ensures { ctx_holds' c2 y z }
1700 | Cons h t -> interp_holds' c1 h y z; impl_holds' c1 t y z
1703 let rec lemma ctx_impl' (c1 c2: context') (g:equality') (y z: vars)
1704 requires { ctx_impl_ctx' c1 c2 }
1706 requires { interp_ctx' c2 g y z }
1707 ensures { interp_ctx' c1 g y z }
1709 = if ctx_holds' c1 y z
1711 impl_holds' c1 c2 y z;
1712 interp_holds' c2 g y z;
1713 holds_interp_ctx' c1 g y z
1716 let rec lemma interp_ctx_impl' (ctx: context') (g1 g2:equality')
1717 requires { forall y z. interp_eq' g1 y z -> interp_eq' g2 y z }
1718 ensures { forall y z. interp_ctx' ctx g1 y z -> interp_ctx' ctx g2 y z }
1720 = match ctx with Nil -> () | Cons _ t -> interp_ctx_impl' t g1 g2 end
1722 let lemma impl_cons (c1 c2:context') (e:equality') (y z:vars)
1723 requires { ctx_impl_ctx' c1 c2 }
1724 requires { forall y z. interp_ctx' c1 e y z }
1725 ensures { ctx_impl_ctx' c1 (Cons e c2) }
1728 let rec lemma impl_wl' (c1 c2:context') (e:equality')
1729 requires { ctx_impl_ctx' c1 c2 }
1730 ensures { ctx_impl_ctx' (Cons e c1) c2 }
1734 | Cons h t -> interp_ctx_wl' c1 (Cons e Nil) h; impl_wl' c1 t e
1737 let rec lemma impl_self (c:context')
1738 ensures { ctx_impl_ctx' c c }
1742 | Cons h t -> (impl_self t; impl_wl' c t h)
1745 predicate is_eq_tbl (a:array (option E.exp)) (l:context')
1746 = forall i. 0 <= i < length a ->
1749 | Some e -> forall y z. y=z -> ctx_holds' l y z
1750 -> E.interp_exp (E.Var i) z = E.interp_exp e z
1756 let prop_ctx (l:context') (g:equality') : (context', equality')
1757 requires { valid_ctx' l }
1758 requires { valid_eq' g }
1759 returns { rl, _ -> length rl = length l }
1760 returns { rl, rg -> valid_ctx' rl /\ valid_eq' rg
1761 /\ forall y z. y=z -> interp_ctx' rl rg y z
1762 -> interp_ctx' l g y z }
1763 returns { rl, rg -> forall y z. y=z -> ctx_holds' l y z
1764 -> pos_ctx' l z -> pos_eq' g z
1765 -> (pos_ctx' rl z /\ pos_eq' rg z) }
1766 raises { OutOfBounds -> true }
1768 let m = max (max_var_ctx' l) (max_var_e' g) in
1769 let a : array (option E.exp) = Array.make (m+1) None in
1770 let rec exp_of_expr' (e:expr') : option E.exp
1771 returns { | None -> true
1772 | Some ex -> forall y z. y=z -> interp' e y z = E.interp_exp ex z }
1775 | Var i -> Some (E.Var i)
1777 let r1,r2 = (exp_of_expr' e1, exp_of_expr' e2) in
1779 | Some ex1, Some ex2 -> Some (E.Plus ex1 ex2)
1783 let r1,r2 = (exp_of_expr' e1, exp_of_expr' e2) in
1785 | Some ex1, Some ex2 -> Some (E.Sub ex1 ex2)
1788 | Coeff (I n) -> Some (E.Lit n)
1792 let fill_tbl_eq (eq:equality') : unit
1793 requires { eq_bound' eq m }
1794 ensures { forall l. is_eq_tbl (old a) l ->
1795 is_eq_tbl a (Cons eq l) }
1797 | Var i, e | e, Var i ->
1798 let r = exp_of_expr' e in
1802 assert { forall l y z. y=z -> ctx_holds' (Cons eq l) y z ->
1803 E.interp_exp ex z = interp' e y z
1804 = interp' (Var i) y z = y i };
1810 let rec fill_tbl_ctx (l:context') : unit
1811 requires { is_eq_tbl a Nil }
1812 ensures { is_eq_tbl a l }
1813 requires { ctx_bound' l m }
1817 | Cons eq l -> fill_tbl_ctx l; fill_tbl_eq eq
1821 (* a contains equalities, let us propagate them so that
1822 only a single pass on the context is needed *)
1823 let seen = Array.make (m+1) false in
1824 let rec propagate_in_tbl (i:int) : unit
1825 requires { is_eq_tbl a l }
1826 ensures { is_eq_tbl a l }
1827 raises { OutOfBounds -> true }
1828 variant { numof seen false 0 (m+1) }
1829 requires { seen[i] = false }
1830 ensures { seen[i] = true }
1831 ensures { forall j. old seen[j] -> seen[j] }
1834 let rec prop (e:E.exp) : E.exp
1835 requires { is_eq_tbl a l }
1836 ensures { is_eq_tbl a l }
1837 ensures { forall y z. y=z -> ctx_holds' l y z ->
1838 E.interp_exp e z = E.interp_exp result z }
1839 ensures { forall j. old seen[j] -> seen[j] }
1840 raises { OutOfBounds -> true }
1841 requires { numof seen false 0 (m+1) < numof (seen at Start) false 0 (m+1) }
1846 if (not (defensive_get seen j)) then propagate_in_tbl j;
1847 match (defensive_get a j) with
1851 | E.Plus e1 e2 -> E.Plus (prop e1) (prop e2)
1852 | E.Sub e1 e2 -> E.Sub (prop e1) (prop e2)
1853 | E.Minus e -> E.Minus (prop e)
1856 defensive_set seen i true;
1857 assert { numof seen false 0 (m+1) < numof (old seen) false 0 (m+1)
1858 by forall j. 0 <= j < m+1 -> (old seen)[j] -> seen[j]
1859 so not (old seen)[i] so seen[i] };
1862 | Some e -> a[i] <- Some (prop e)
1866 invariant { is_eq_tbl a l }
1867 if not seen[i] then propagate_in_tbl i;
1869 let rec propagate_exp (e:E.exp)
1870 ensures { forall y z. y=z -> ctx_holds' l y z ->
1871 E.interp_exp e z = E.interp_exp result z }
1873 raises { OutOfBounds -> true }
1876 | E.Var i -> match (defensive_get a i) with Some e' -> e' | None -> e end
1877 | E.Plus e1 e2 -> E.Plus (propagate_exp e1) (propagate_exp e2)
1878 | E.Sub e1 e2 -> E.Sub (propagate_exp e1) (propagate_exp e2)
1879 | E.Minus e -> E.Minus (propagate_exp e)
1882 let propagate_coeff (c:t)
1883 ensures { forall y z. y=z -> ctx_holds' l y z ->
1884 interp_eq' (Coeff c, Coeff result) y z }
1885 ensures { forall y z. y = z -> ctx_holds' l y z ->
1886 pos_exp c z -> pos_exp result z }
1887 raises { OutOfBounds -> true }
1890 | E e -> E (propagate_exp e)
1894 let rec propagate_c (c:cprod)
1895 ensures { forall y z. y=z -> ctx_holds' l y z ->
1896 interp_c c y z = interp_c result y z }
1898 raises { OutOfBounds -> true }
1899 ensures { forall y z. y = z -> ctx_holds' l y z ->
1900 pos_cprod c z -> pos_cprod result z }
1902 | C c -> C (propagate_coeff c)
1903 | Times c1 c2 -> Times (propagate_c c1) (propagate_c c2)
1906 let rec propagate_e (e:expr')
1907 requires { expr_bound' e m }
1908 ensures { expr_bound' result m }
1909 ensures { forall y z. y=z -> ctx_holds' l y z -> interp_eq' (e,result) y z }
1911 raises { OutOfBounds -> true }
1912 requires { valid_expr' e }
1913 ensures { valid_expr' result }
1914 ensures { forall y z. y = z -> ctx_holds' l y z
1915 -> pos_expr' e z -> pos_expr' result z }
1918 | ProdL e c -> ProdL (propagate_e e) (propagate_c c)
1919 | ProdR c e -> ProdR (propagate_c c) (propagate_e e)
1920 | Sum e1 e2 -> Sum (propagate_e e1) (propagate_e e2)
1921 | Diff e1 e2 -> Diff (propagate_e e1) (propagate_e e2)
1922 | Coeff c -> Coeff (propagate_coeff c)
1925 let rec propagate_eq (eq:equality')
1926 requires { eq_bound' eq m }
1927 ensures { eq_bound' result m }
1928 ensures { forall y z. y=z -> interp_ctx' l eq y z <-> interp_ctx' l result y z }
1929 raises { OutOfBounds -> true }
1930 requires { valid_eq' eq }
1931 ensures { valid_eq' result }
1932 ensures { forall y z. y = z -> ctx_holds' l y z
1933 -> pos_eq' eq z -> pos_eq' result z }
1934 = match eq with (a,b) -> (propagate_e a, propagate_e b) end
1936 let rec propagate (acc:context') : context'
1937 requires { ctx_bound' acc m }
1938 ensures { ctx_bound' result m }
1939 requires { ctx_impl_ctx' l acc }
1940 ensures { ctx_impl_ctx' l result }
1941 ensures { length result = length acc }
1943 requires { valid_ctx' acc }
1944 ensures { valid_ctx' result }
1945 ensures { forall y z. y = z -> ctx_holds' l y z
1946 -> pos_ctx' acc z -> pos_ctx' result z }
1947 raises { OutOfBounds -> true }
1951 let h' = propagate_eq h in
1952 let t' = propagate t in
1956 propagate l, propagate_eq g
1958 use LinearDecisionRationalMP as R
1960 let prop_mp_decision (l:context') (g:equality') : bool
1961 requires { valid_ctx' l }
1962 requires { valid_eq' g }
1963 requires { length l < 100000 }
1964 ensures { forall y z. result -> pos_ctx' l z -> pos_eq' g z
1965 -> y = z -> interp_ctx' l g y z }
1966 raises { | OutOfBounds -> true | E.MPError -> true
1967 | E.Q.QError -> true | R.Failure -> true}
1968 = let l', g' = prop_ctx l g in
1980 meta "compute_max_steps" 0x10000
1982 goal g: forall i x c r: int.
1984 x + (2 * (power radix i) * c) = r ->
1985 radix * x + (2 * (power radix (i+1)) * c) = radix * r
1987 goal g': forall a b i j: int.
1989 (power radix i) * a = b ->
1991 (power radix j) * a = radix*b
1993 goal g'': forall r r' i c x x' y l: int.
1996 r + power radix i * c = x + y ->
1997 r' = r + power radix i * l ->
1998 x' = x + power radix i * l ->
1999 r' + power radix (i+1) * c = x' + y
2001 (*tries to add power radix i and power radix (i+1), fails
2002 -> cst propagation ? *)
2009 use LinearDecisionInt
2011 meta "compute_max_steps" 0x10000
2013 goal g: forall x y z: int.
2026 constant dummy : value
2028 function add value value : value
2030 type term = Val int | Add term term
2031 type fmla = Forall fmla | Foo term
2033 function interp_term (t:term) (b:int->value) : value =
2036 | Add t1 t2 -> add (interp_term t1 b) (interp_term t2 b)
2039 function interp_fmla (f:fmla) (l:int) (b:int->value) : bool =
2041 | Foo t -> foo (interp_term t b)
2042 | Forall f -> forall v. interp_fmla f (l-1) b[l <- v]
2045 function interp (f:fmla) (b: int -> value) : bool =
2046 interp_fmla f (-1) b
2048 let f (f:fmla) : bool
2049 ensures { result -> forall b. interp f b }
2057 meta "compute_max_steps" 0x10000
2061 ((forall x. forall y. foo (add x (add (add a dummy) y))) = True)