Use the proper dependency for lablgtk3-sourceview3.
[why3.git] / examples / multiprecision / lineardecision.mlw
blobafc81239b0afe694210747e4143df4a2c1d9849b
1 module LinearEquationsCoeffs
3 type a
4 function (+) a a : a
5 function (*) a a : a
6 function (-_) a : a
7 function azero: a
8 function aone: a
9 predicate ale a a
11 clone algebra.OrderedUnitaryCommutativeRing as A with
12   type t = a, function (+) = (+), function (*) = (*),
13   function (-_) = (-_), constant zero = azero,
14   constant one=aone, predicate (<=) = ale,
15   axiom .
17 function (-) a a : a
19 axiom sub_def: forall a1 a2. a1 - a2 = a1 + (- a2)
21 type t
22 type vars = int -> a
23 type cvars
24 exception Unknown
26 function interp t cvars : a
28 val constant czero : t
29 val constant cone : t
31 axiom zero_def: forall y. interp czero y = azero
32 axiom one_def: forall y. interp cone y = aone
34 lemma neg_mul:
35   forall x y: a. (-x) * y = - (x*y)
37 val add (a b: t) : t
38   ensures { forall v: cvars. interp result v = interp a v + interp b v }
39   raises  { Unknown -> true }
41 val mul (a b: t) : t
42   ensures { forall v: cvars. interp result v = interp a v * interp b v }
43   raises  { Unknown -> true }
45 val opp (a:t) : t
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 }
51 val inv (a:t) : t
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 }
57 end
59 module LinearEquationsDecision
61 use int.Int
62 type coeff
64 clone LinearEquationsCoeffs as C with type t = coeff, axiom .
65 type vars = C.vars
67 type expr = Term coeff int | Add expr expr | Cst coeff
69 let rec predicate valid_expr (e:expr)
70   variant { e }
71 = match e with
72   | Term _ i -> 0 <= i
73   | Cst _ -> true
74   | Add e1 e2 -> valid_expr e1 && valid_expr e2
75   end
77 let rec predicate expr_bound (e:expr) (b:int)
78   variant { e }
79 = match e with
80   | Term _ i -> 0 <= i <= b
81   | Cst _ -> true
82   | Add e1 e2 -> expr_bound e1 b && expr_bound e2 b
83   end
85 function interp (e:expr) (y:vars) (z:C.cvars) : C.a
86 = match e with
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
90   end
92 use bool.Bool
93 use list.List
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 }
114   variant  { e }
115 = match e with
116   | Add e1 e2 -> expr_bound_w e1 b1 b2; expr_bound_w e2 b1 b2
117   | Cst _ -> ()
118   | Term _ _ -> ()
119   end
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 }
127   variant  { l }
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
134 = match l with
135   | Nil -> interp_eq g y z
136   | Cons h t -> (interp_eq h y z) -> (interp_ctx t g y z)
137   end
139 use mach.int.Int63
140 use seq.Seq
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]);
152     done
153   done;
154   r
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]);
164     done
165   done;
166   r
168 use ref.Ref
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]);
176   done;
177   !r
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)
198     done;
199     set r i m.columns v[i]
200   done;
201   r
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 }
212     r[i] <- v[i]
213   done;
214   r
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
221   else
222     let r = ref true in
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;
226     done;
227     !r
229 use int.MinMax
230 use list.Length
232 let rec function max_var (e:expr) : int
233   variant { e }
234   requires { valid_expr e }
235   ensures { 0 <= result }
236   ensures { expr_bound e result }
237 = match e with
238   | Term _ i -> i
239   | Cst _ -> 0
240   | Add e1 e2 -> max (max_var e1) (max_var e2)
241   end
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
250   variant { l }
251   requires { valid_ctx l }
252   ensures { 0 <= result }
253   ensures { ctx_bound l result }
254 = match l with
255   | Nil -> 0
256   | Cons e t -> max (max_var_e e) (max_var_ctx t)
257   end
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 }
263   variant { e }
264 = match e with
265   | Cst c -> Cst (C.opp c)
266   | Term c j ->
267     let oc = C.opp c in
268     let r = Term oc j in
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) };
273     r
274   | Add e1 e2 ->
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) };
283       Add e1' e2'
284   end
286 predicate atom (e:expr)
287 = match e with
288   | Add _ _ -> false | _ -> true
289   end
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 }
300   variant { ex }
301 = match ex with
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
305                  norm_eq_aux e2 ae ac
306   end
308 use debug.Debug
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 }
315 = match e with
316   | (e1, e2) ->
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
320       (e, c) ->
321         let ec = C.opp c in
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 };
328         e, ec
329     end
330   end
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 }
336   variant  { ctx }
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 }
341   variant  { ctx }
342 = match ctx with Nil -> () | Cons _ t -> interp_ctx_valid t g end
344 use list.Append
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 }
348   variant { ctx }
349 = match ctx with
350   | Nil -> ()
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 }
355   variant { l }
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 }
362   variant { e }
363   raises  { C.Unknown -> true }
364 = if C.eq c C.czero then Cst C.czero
365   else match e with
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)
369   end
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) }
374   variant { e2 }
375   raises  { C.Unknown -> true }
377   let term_or_cst c i
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)
381     requires { atom a }
382     returns { r,_ -> forall y z. interp r y z
383                      = C.(+) (interp e y z) (interp a y z) }
384     variant { e }
385     raises  { C.Unknown -> true }
386   = match (e,a) with
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)
401     | Add e1 e2, _ ->
402       let r, b = add_atom e1 a in
403       if b
404       then
405         match r with
406           | Cst c ->
407             if C.eq c C.czero
408             then begin
409               assert { forall y z. C.(+) (interp e1 y z) (interp a y z) = C.azero };
410               e2, True end
411             else Add r e2, True
412           | _ -> Add r e2, True
413         end
414       else
415         let r,b = add_atom e2 a in
416         match r with
417           | Cst c ->
418             if C.eq c C.czero
419             then begin
420               assert { forall y z. C.(+) (interp e2 y z) (interp a y z) = C.azero };
421               e1, True end
422             else Add e1 r, b
423           | _ -> Add e1 r, b
424         end
425     | _, Add _ _ -> absurd
426     end
427   in
428   match e2 with
429     | Add e1' e2' -> add_expr (add_expr e1 e1') e2'
430     | _ -> let r,_= add_atom e1 e2 in r
431   end
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
446   let r = (a,b) 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 }
450     variant { l }
451   = match l with Nil -> () | Cons _ t -> aux t end in
452   r
453   end
455 let rec zero_expr (e:expr) : bool
456   ensures { result -> forall y z. interp e y z = C.azero }
457   variant { e }
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 }
462     variant { e }
463     = match e with
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
467     end
468   in
469   let e' = add_expr (Cst C.czero) e in (* simplifies expr *)
470   all_zero e'
472 let sub_expr (e1 e2:expr)
473   ensures { forall y z. C.(+) (interp result y z) (interp e2 y z)
474                         = interp e1 y z }
475   raises  { C.Unknown -> true }
476 = let r = add_expr e1 (mul_expr e2 (C.opp C.cone)) in
477   assert { forall y z.
478            let v1 = interp e1 y z in
479            let v2 = interp e2 y z in
480            let vr = interp r y z in
481            C.(+) vr v2 = v1
482            by C.(*) v2 (C.(-_) C.aone) = C.(-_) v2
483            so C.(+) vr v2
484            = C.(+) (C.(+) v1 (C.(*) v2 (C.(-_) C.aone))) v2
485            = C.(+) (C.(+) v1 (C.(-_) v2)) v2 = v1 };
486   r
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)
497 use option.Option
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 }
502   variant { l }
503 = match l with
504   | Nil -> Nil
505   | Cons h t ->
506     let ex, c = norm_eq h in
507     Cons (ex, Cst c) (norm_context t)
508   end
510 let rec print_lc ctx v : unit variant { ctx }
511 = match ctx, v with
512   | Nil, Nil -> ()
513   | Cons l t, Cons v t2 ->
514    (if C.eq C.czero v then ()
515     else (print l; print v));
516     print_lc t t2
517   | _ -> ()
518   end
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 }
532     variant { l }
533   = match (l, v) with
534     | Nil, Nil -> Some s
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
538       else begin
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
545     | _ -> None
546     end
547   in
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)
550   | None -> false
551   end
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)
560     done
561   done;
562   r
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);
569     set m i2 j c
570   done
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))
579   done
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)))
587   done
589 use ref.Ref
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 }
597   let n = a.rows in
598   let m = a.columns in
599   (* print n; print m; *)
600   let rec find_nonz (i j:int63)
601     requires { 0 <= i <= n }
602     requires { 0 <= j < m }
603     variant { n-i }
604     ensures { i <= result <= n }
605     ensures { result < n -> not (C.eq (a.elts result j) C.czero) }
606     = if i >= n then n
607     else
608       if C.eq (get a i j) C.czero
609       then find_nonz (i+1) j
610       else i in
611   let pivots = Array63.make n 0 in
612   let r = ref (-1) in
613   for j = 0 to m-2 do
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 }
618     label Start in
619     let k = find_nonz (!r+1) j in
620     if k < n
621     then begin
622       r := !r + 1;
623       pivots[!r] <- j;
624       mul_row a k (C.inv(get a k j));
625       if k <> !r then swap_rows a k !r;
626       for i = 0 to n-1 do
627         if i <> !r
628         then addmul_row a !r i (C.opp(get a i j))
629       done;
630     end;
631     assert { forall i1 i2: int. 0 <= i1 < i2 <= !r -> pivots[i1] < pivots[i2]
632              by pivots[i1] = pivots[i1] at Start
633              so [@case_split]
634                 ((i2 < !r so pivots[i2] = pivots[i2] at Start)
635                 \/ (i2 = !r so pivots[i1] < j(* = pivots[i2])*))) };
636   done;
637   if !r < 0 then None (* matrix is all zeroes *)
638   else begin
639     let v = Array63.make m(*(m-1)*) C.czero in
640     for i = 0 to !r do
641       v[pivots[i]] <- get a i (m-1)
642     done;
643     Some v (*pivots[!r] < m-1*)  (*pivot on last column, no solution*)
644   end
646 let rec function to_list (a: array 'a) (l u: int63) : list 'a
647   requires { l >= 0 /\ u <= Array63.length a }
648   variant  { u - l }
649 = if u <= l then Nil else Cons a[l] (to_list a (l+1) u)
651 exception Failure
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
663   end;
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
670     variant { ex }
671     raises  { C.Unknown -> true }
672     requires { 0 <= i < length l }
673     requires { expr_bound ex nv }
674     raises  { Failure -> true }
675   = match ex with
676     | Cst c -> if C.eq c C.czero then () else raise Failure
677     | Term c j ->
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
681     end in
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 }
688   = match ctx with
689     | Nil -> ()
690     | Cons e t ->
691       assert { i < length l };
692       try
693         let ex, c = norm_eq e in
694         if (not (C.eq c C.czero)) then b[i] <- C.add b[i] c;
695         fill_expr ex i;
696       with C.Unknown -> () (* some equalities are in the context but cannot be normalized, typically they are useless, ignore them *)
697       end;
698       fill_ctx t (i+1)
699     end in
700   let rec fill_goal (ex:expr) : unit
701     requires { expr_bound ex nv }
702     variant { ex }
703     raises { C.Unknown -> true }
704     raises  { Failure -> true }
705   = match ex with
706     | Cst c -> if C.eq c C.czero then () else raise Failure
707     | Term c j ->
708       let j = Int63.of_int j in
709       v[j] <- C.add v[j] c
710     | Add e1 e2 -> fill_goal e1; fill_goal e2
711     end in
712   fill_ctx l 0;
713   let (ex, d) = norm_eq g in
714   fill_goal ex;
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
719     | Some r ->
720       check_combination l g (to_list r 0 ll)
721     | None -> false
722   end
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
730 = match e with
731   | C c -> C.interp c z
732   | Times e1 e2 -> C.(*) (interp_c e1 y z) (interp_c e2 y z)
733   end
735 function interp' (e:expr') (y:vars) (z:C.cvars) : C.a
736 = match e with
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)
741   | Var n -> y n
742   | Coeff c -> C.interp c z
743   end
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
754 = match l with
755   | Nil -> interp_eq' g y z
756   | Cons h t -> (interp_eq' h y z) -> (interp_ctx' t g y z)
757   end
759 let rec predicate valid_expr' (e:expr')
760   variant { e }
761 = match e with
762   | Var i -> 0 <= i
763   | Sum e1 e2 | Diff e1 e2 -> valid_expr' e1 && valid_expr' e2
764   | Coeff _ -> true
765   | ProdL e _ | ProdR _ e -> valid_expr' e
766   end
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 }
778   variant { e }
780   let rec simp_c (e:cprod) : coeff
781     ensures { forall y z. C.interp result z = interp_c e y z }
782     variant { e }
783     raises  { C.Unknown -> true }
784   =
785     match e with
786     | C c -> c
787     | Times c1 c2 -> C.mul (simp_c c1) (simp_c c2)
788     end
789   in
790   match e with
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
794   | Coeff c -> Cst c
795   | ProdL e c | ProdR c e ->
796     mul_expr (simp e) (simp_c c)
797   end
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 }
811   variant { ctx }
812 = match ctx with
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
816   end
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
831 use int.Int
832 use real.RealInfix
833 use real.FromInt
834 use int.Abs
836 type t = (int, int)
837 type rvars = int -> real
839 exception QError
841 let constant rzero = (0,1)
842 let constant rone = (1,1)
844 function rinterp (t:t) (v:rvars) : real
845 = match t with
846   | (n,d) ->  from_int n /. from_int d
847   end
849 let lemma prod_compat_eq (a b c:real)
850   requires { c <> 0.0 }
851   requires { a *. c = b *. c }
852   ensures  { a = b }
853 = ()
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
882 use ref.Ref
883 use number.Gcd
885 let gcd (x:int) (y:int)
886   requires { x > 0 /\ y > 0 }
887   ensures { result = gcd x y }
888   ensures { result > 0 }
889   =
890   let ghost ox = x in
891   let x = ref x in let y = ref y in
892   label Pre in
893   while (!y > 0) do
894      invariant { !x >= 0 /\ !y >= 0 }
895      invariant { gcd !x !y = gcd (!x at Pre) (!y at Pre) }
896      variant { !y }
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 };
900      x := !y; y := r;
901   done;
902   !x
904 let simp (t:t) : t
905   ensures { forall v:rvars. rinterp result v = rinterp t v }
906 = match t with
907   | (n,d) ->
908     if d = 0 then t
909     else if n = 0 then rzero
910     else
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' };
915     (n', d')
916   end
918 let radd (a b:t)
919   ensures { forall y. rinterp result y = rinterp a y +. rinterp b y }
920   raises  { QError -> true }
921 = match (a,b) with
922   | (n1,d1), (n2,d2) ->
923   if d1 = 0 || d2 = 0 then raise QError
924   else begin
925     let r = (n1*d2 + n2*d1, d1*d2) in
926     let ghost d = from_int d1 *. from_int d2 in
927     assert { forall y.
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 };
934     simp r end
935  end
937 let rmul (a b:t)
938   ensures { forall y. rinterp result y = rinterp a y *. rinterp b y }
939   raises  { QError -> true }
940 = match (a,b) with
941   | (n1,d1), (n2, d2) ->
942     if d1 = 0 || d2 = 0 then raise QError
943     else begin
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 };
950       r
951     end
952   end
954 let ropp (a:t)
955   ensures { forall y. rinterp result y = -. rinterp a y }
956 = match a with
957   | (n,d) -> (-n, d)
958   end
960 let predicate req (a b:t)
961   ensures { result -> forall y. rinterp a y = rinterp b y }
962 = match (a,b) with
963   | (n1,d1), (n2,d2) -> n1 = n2 && d1 = d2 || (d1 <> 0 && d2 <> 0 && n1 * d2 = n2 * d1)
964   end
966 let rinv (a:t)
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 }
971 = match a with
972   | (n,d) -> if n = 0 || d = 0 then raise QError else (d,n)
973   end
975 let is_zero (a:t)
976   ensures { result <-> req a rzero }
977 = match a with
978   | (n,d) -> n = 0 && d <> 0
979   end
983 module LinearDecisionRational
985 use RationalCoeffs
986 use real.RealInfix
987 use real.FromInt
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
995 use int.Int
997 type t' =  IC int | Error
999 function interp_id (t:t') (v:int -> int) : int
1000 = match t with
1001   | IC i -> i
1002   | Error -> 0 (* never created *)
1003   end
1005 let constant izero = IC 0
1007 let constant ione = IC 1
1009 let predicate ieq (a b:t') = false
1011 exception NError
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 }
1016 = raise NError
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 }
1021 = raise NError
1023 let iopp (a:t') : t'
1024   ensures { forall z. interp_id result z = - interp_id a z }
1025   raises  { NError -> true }
1026 = raise NError
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 }
1032 = raise NError
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 .
1037 use real.FromInt
1038 use RationalCoeffs
1039 use LinearDecisionRational as R
1040 use list.List
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 }
1049 = match t with
1050   | IC x -> (x,1)
1051   | _ -> raise NError
1052   end
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 }
1058   variant { e }
1059 = match e with
1060   | C c -> R.C (m c)
1061   | Times c1 c2 -> R.Times (m_cprod c1) (m_cprod c2)
1062   end
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 }
1069   variant { e }
1070 = match e with
1071   | Var i -> R.Var i
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)
1077   end
1079 use list.Length
1080 use debug.Debug
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 }
1096   variant { ctx }
1097 = match ctx with
1098   | Nil -> Nil, m_eq g
1099   | Cons h t ->
1100     let c',g' = m_ctx t g in
1101     (Cons (m_eq h) c',g')
1102     end
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
1111   R.decision l' g'
1116 module Test
1118 use RationalCoeffs
1119 use LinearDecisionRational
1120 use int.Int
1121 use real.RealInfix
1122 use real.FromInt
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)
1133 module TestInt
1135 use LinearDecisionInt
1136 use int.Int
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 ->
1143      x = 5
1147 module MP64Coeffs
1149 use mach.int.UInt64 as M
1150 use real.RealInfix
1151 use real.FromInt
1152 use real.PowerReal
1153 use RationalCoeffs as Q
1154 use int.Int
1156 use debug.Debug
1158 type evars = int -> int
1161 type exp = Lit int | Var int | Plus exp exp | Minus exp | Sub exp exp
1162 type t = (Q.t, 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
1175 = match e with
1176   | Lit n -> n
1177   | Var v -> y v
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)
1181   end
1183 function minterp (t:t) (y:evars) : real
1184 = match t with
1185   (q,e) ->
1186   qinterp q *. pow rradix (from_int (interp_exp e y))
1187   end
1189 exception MPError
1191 let rec opp_exp (e:exp)
1192   ensures { forall y. interp_exp result y = - interp_exp e y }
1193   variant { e }
1194 = match e with
1195   | Lit n -> Lit (-n)
1196   | Minus e' -> e'
1197   | Plus e1 e2 -> Plus (opp_exp e1) (opp_exp e2)
1198   | Sub e1 e2 -> Sub e2 e1
1199   | Var _ -> Minus e
1200   end
1202 let rec add_sub_exp (e1 e2:exp) (s:bool) : exp
1203   ensures { forall y.
1204             if s
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 }
1208   variant { e2, e1 }
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 }
1215     variant { e }
1216   = match (e,a) with
1217     | Lit n1, Lit n2 -> (if s then Lit (n1+n2) else Lit (n1-n2)), True
1218     | Lit n, Var i
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
1221     | Var i, Lit n
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
1234                       else
1235                         if i = j then Lit 0, True
1236                         else Sub e a, False
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
1241     | Plus e1 e2, _ ->
1242       let r, b = add_atom e1 a s in
1243       if b then
1244         match r with
1245         | Lit n -> if n = 0 then e2, True else Plus r e2, True
1246         | _ -> Plus r e2, True
1247         end
1248       else let r, b = add_atom e2 a s in Plus e1 r, b
1249     | Sub e1 e2, _ ->
1250       let r, b = add_atom e1 a s in
1251       if b then
1252         match r with
1253         | Lit n -> if n = 0 then opp_exp e2, True else Sub r e2, True
1254         | _ -> Sub r e2, True
1255         end
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
1261     end
1262   in
1263   match e2 with
1264    | Plus e1' e2' ->
1265      let r = add_sub_exp e1 e1' s in
1266      match r with
1267      | Lit n -> if n = 0
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
1271      end
1272    | Sub e1' e2' ->
1273      let r = add_sub_exp e1 e1' s in
1274      match r with
1275      | Lit n -> if n = 0
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)
1279      end
1280    | _ -> let r, _ = add_atom e1 e2 s in r
1281   end
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 }
1291   variant { e }
1292   raises  { MPError -> true }
1294   let rec all_zero (e:exp) : bool
1295     ensures { result -> forall y. interp_exp e y = 0 }
1296     variant { e }
1297   = match e with
1298     | Lit n -> n = 0
1299     | Var _ -> false
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
1303     end
1304   in
1305   let e' = add_exp (Lit 0) e in (* simplifies exp *)
1306   all_zero e'
1308 let rec same_exp (e1 e2: exp)
1309   ensures { result -> forall y. interp_exp e1 y = interp_exp e2 y }
1310   variant { e1, e2 }
1311   raises  { MPError -> true }
1312 = match e1, e2 with
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))
1317   end
1319 let madd (a b:t)
1320   ensures { forall y. minterp result y = minterp a y +. minterp b y }
1321   raises  { MPError -> true }
1322   raises  { Q.QError -> true }
1323 = match a, b with
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
1328     then begin
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 };
1336       (q,e1) end
1337     else raise MPError
1338   end
1340 let mmul (a b:t)
1341   ensures { forall y. minterp result y = minterp a y *. minterp b y }
1342   raises  { Q.QError -> true }
1343   raises  { MPError -> true }
1344 = match a, b with
1345   | (q1,e1), (q2,e2) ->
1346     let q = Q.rmul q1 q2 in
1347     if Q.is_zero q then mzero
1348     else begin
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
1355                   so p = p1 *. p2
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 };
1360       (q,e)
1361     end
1362   end
1364 let mopp (a:t)
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 }
1370   variant { e1, e2 }
1371 = match e1, e2 with
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)
1378   | _ -> false
1379   end
1381 let predicate meq (a b:t)
1382   ensures { result -> forall y. minterp a y = minterp b y }
1383 = match (a,b) with
1384   | (q1,e1), (q2,e2) -> (Q.req q1 q2 && pure_same_exp e1 e2) || (Q.is_zero q1 && Q.is_zero q2)
1385   end
1387 let minv (a:t)
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 }
1392 = match a with
1393   | (q,e) -> (Q.rinv q, opp_exp e)
1394   end
1398 module LinearDecisionRationalMP
1400 use MP64Coeffs
1401 use real.RealInfix
1403 type coeff = t
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
1410 use int.Int
1411 use int.Power
1412 use MP64Coeffs
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
1420 = match t with
1421   | I n -> n
1422   | E e -> power M.radix (interp_exp e y)
1423   | R -> M.radix
1424   end
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 }
1431 = raise MPError
1433 let mpmul (a b:t) : t
1434   ensures { forall y. mpinterp result y = mpinterp a y * mpinterp b y }
1435   raises  { MPError -> true }
1436 = raise MPError
1438 let mpopp (a:t) : t
1439   ensures { forall y. mpinterp result y = - mpinterp a y }
1440   raises  { MPError -> true }
1441 = raise MPError
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)
1447   end*)
1449 let mpinv (a:t) : t
1450   ensures { not mpeq result mpzero }
1451   raises  { MPError -> true }
1452 = raise MPError
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
1458 use real.FromInt
1459 use real.PowerReal
1460 use real.RealInfix
1461 use int.Int
1463 use list.List
1465 predicate pos_exp (t:t) (y:evars)
1466 = match t with
1467   | E e -> 0 <= interp_exp e y
1468   | I _ | R -> true
1469   end
1471 predicate pos_cprod (e:cprod) (y:evars)
1472 = match e with
1473   | C c -> pos_exp c y
1474   | Times c1 c2 -> pos_cprod c1 y && pos_cprod c2 y
1475   end
1477 predicate pos_expr' (e:expr') (y:evars)
1478 = match e with
1479   | Coeff c -> pos_exp c y
1480   | Var _ -> true
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
1484   end
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) }
1495 = match t with
1496   | I n -> ((n,1), Lit 0)
1497   | E e -> ((1,1), e)
1498   | R -> ((1,1), Lit 1) (* or ((radix, 1), Lit 0) ? *)
1499   end
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) }
1508 = match e with
1509   | C c -> R.C (m c)
1510   | Times c1 c2 -> R.Times (m_cprod c1) (m_cprod c2)
1511   end
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}
1517 = match e with
1518   | Var i -> R.Var i
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)
1524   end
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
1532 use list.Length
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 }
1540   variant { ctx }
1541 = match ctx with
1542   | Nil -> Nil
1543   | Cons h t ->
1544     let r = Cons (m_eq h) (m_ctx t) in
1545     r
1546     end
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)
1560 module EqPropMP
1562 use int.Int
1563 use LinearDecisionIntMP
1564 use array.Array
1565 use int.MinMax
1566 use option.Option
1567 use list.List
1568 use list.Append
1571 use MP64Coeffs as E
1573 let rec predicate expr_bound' (e:expr') (b:int)
1574   variant { e }
1575 = match e with
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
1579   | Coeff _ -> true
1580   end
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 }
1592   variant  { e }
1593 = match e with
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
1597   | _ -> ()
1598   end
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 }
1606   variant  { l }
1607 = match l with Nil -> () | Cons _ t -> ctx_bound_w' t b1 b2 end
1609 let rec function max_var' (e:expr') : int
1610   variant { e }
1611   requires { valid_expr' e }
1612   ensures { 0 <= result }
1613   ensures { expr_bound' e result }
1614 = match e with
1615   | Var i -> i
1616   | Coeff _ -> 0
1617   | Sum e1 e2 | Diff e1 e2 -> max (max_var' e1) (max_var' e2)
1618   | ProdL e _ | ProdR _ e -> max_var' e
1619   end
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
1628   variant { l }
1629   requires { valid_ctx' l }
1630   ensures { 0 <= result }
1631   ensures { ctx_bound' l result }
1632 = match l with
1633   | Nil -> 0
1634   | Cons e t -> max (max_var_e' e) (max_var_ctx' t)
1635   end
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 }
1639   variant  { ctx }
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 }
1644   variant { ctx }
1645 = match ctx with
1646   | Nil -> ()
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 }
1651   variant { l }
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 }
1658 = ()
1660 predicate ctx_impl_ctx' (c1 c2: context')
1661 = match c2 with
1662   | Nil -> true
1663   | Cons eq t -> ctx_impl_ctx' c1 t /\ forall y z. y=z -> interp_ctx' c1 eq y z
1664   end
1666 predicate ctx_holds' (c: context') (y z:vars)
1667 = match c with
1668   | Nil -> true
1669   | Cons h t -> interp_eq' h y z /\ ctx_holds' t y z
1670   end
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 }
1675   variant { l }
1676 = match l with
1677   | Nil -> ()
1678   | Cons h t ->
1679     if interp_eq' h y z then holds_interp_ctx' t g y z
1680   end
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 }
1686   variant  { l }
1687 = match l with
1688   | Nil -> ()
1689   | Cons _ t -> interp_holds' t g y z
1690   end
1692 let rec lemma impl_holds' (c1 c2: context') (y z: vars)
1693   requires { ctx_impl_ctx' c1 c2 }
1694   requires { y=z }
1695   requires { ctx_holds' c1 y z }
1696   ensures  { ctx_holds' c2 y z }
1697   variant  { c2 }
1698 = match c2 with
1699   | Nil -> ()
1700   | Cons h t -> interp_holds' c1 h y z; impl_holds' c1 t y z
1701   end
1703 let rec lemma ctx_impl' (c1 c2: context') (g:equality') (y z: vars)
1704   requires { ctx_impl_ctx' c1 c2 }
1705   requires { y=z }
1706   requires { interp_ctx' c2 g y z }
1707   ensures  { interp_ctx' c1 g y z }
1708   variant  { c2 }
1709 = if ctx_holds' c1 y z
1710   then begin
1711        impl_holds' c1 c2 y z;
1712        interp_holds' c2 g y z;
1713        holds_interp_ctx' c1 g y z
1714        end
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 }
1719   variant  { ctx }
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) }
1726 = ()
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 }
1731   variant  { c2 }
1732 = match c2 with
1733   | Nil -> ()
1734   | Cons h t -> interp_ctx_wl' c1 (Cons e Nil) h; impl_wl' c1 t e
1735   end
1737 let rec lemma impl_self (c:context')
1738   ensures { ctx_impl_ctx' c c }
1739   variant { c }
1740 = match c with
1741   | Nil -> ()
1742   | Cons h t -> (impl_self t; impl_wl' c t h)
1743   end
1745 predicate is_eq_tbl (a:array (option E.exp)) (l:context')
1746 = forall i. 0 <= i < length a ->
1747   match a[i] with
1748   | None -> true
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
1751   end
1752 use int.NumOf
1753 use array.NumOfEq
1754 use list.Length
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 }
1773     variant { e }
1774   = match e with
1775     | Var i -> Some (E.Var i)
1776     | Sum e1 e2 ->
1777       let r1,r2 = (exp_of_expr' e1, exp_of_expr' e2) in
1778       match r1,r2 with
1779       | Some ex1, Some ex2 -> Some (E.Plus ex1 ex2)
1780       | _ -> None
1781       end
1782     | Diff e1 e2 ->
1783       let r1,r2 = (exp_of_expr' e1, exp_of_expr' e2) in
1784       match r1,r2 with
1785       | Some ex1, Some ex2 -> Some (E.Sub ex1 ex2)
1786       | _ -> None
1787       end
1788     | Coeff (I n) -> Some (E.Lit n)
1789     | _ -> None
1790     end
1791   in
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) }
1796   = match eq with
1797     | Var i, e | e, Var i ->
1798       let r = exp_of_expr' e in
1799       match r with
1800         | None -> ()
1801         | Some ex ->
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 };
1805           a[i] <- Some ex
1806       end
1807     | _ -> ()
1808     end
1809   in
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 }
1814     variant { l }
1815   = match l with
1816     | Nil -> ()
1817     | Cons eq l -> fill_tbl_ctx l; fill_tbl_eq eq
1818     end
1819   in
1820   fill_tbl_ctx l;
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] }
1832   =
1833     label Start in
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) }
1842       variant { e }
1843       = match e with
1844         | E.Lit _ -> e
1845         | E.Var j ->
1846           if (not (defensive_get seen j)) then propagate_in_tbl j;
1847           match (defensive_get a j) with
1848           | None -> e
1849           | Some e' -> e'
1850           end
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)
1854         end
1855     in
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] };
1860     match a[i] with
1861     | None -> ()
1862     | Some e -> a[i] <- Some (prop e)
1863     end;
1864   in
1865   for i = 0 to m do
1866     invariant { is_eq_tbl a l }
1867     if not seen[i] then propagate_in_tbl i;
1868   done;
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 }
1872     variant { e }
1873     raises { OutOfBounds -> true }
1874   = match e with
1875     | E.Lit _ -> e
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)
1880     end
1881   in
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 }
1888   = match c with
1889     | I _ -> c
1890     | E e -> E (propagate_exp e)
1891     | R -> R
1892     end
1893   in
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 }
1897     variant { c }
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 }
1901   = match c with
1902     | C c -> C (propagate_coeff c)
1903     | Times c1 c2 -> Times (propagate_c c1) (propagate_c c2)
1904     end
1905   in
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 }
1910     variant { e }
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 }
1916   = match e with
1917     | Var _ -> e
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)
1923     end
1924   in
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
1935   in
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 }
1942     variant { 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 }
1948   = match acc with
1949      | Nil -> Nil
1950      | Cons h t ->
1951        let h' = propagate_eq h in
1952        let t' = propagate t in
1953        Cons h' t'
1954     end
1955   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
1969     mp_decision l' g'
1973 module TestMP
1975 use EqPropMP
1976 use mach.int.UInt64
1977 use int.Int
1978 use int.Power
1980 meta "compute_max_steps" 0x10000
1982 goal g: forall i x c r: int.
1983   0 <= i ->
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.
1988   0 <= i -> 0 <= j ->
1989   (power radix i) * a = b ->
1990   i+1 = j ->
1991   (power radix j) * a = radix*b
1993 goal g'': forall r r' i c x x' y l: int.
1994   0 <= i ->
1995   c = 0 ->
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 ? *)
2006 module Test2
2008 use int.Int
2009 use LinearDecisionInt
2011 meta "compute_max_steps" 0x10000
2013 goal g: forall x y z: int.
2014   x + y = 0 ->
2015   y - z = 0 ->
2016   x = 0
2020 module Fmla
2022 use map.Map
2023 use int.Int
2025 type value
2026 constant dummy : value
2027 predicate foo 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 =
2034   match t with
2035   | Val n -> b n
2036   | Add t1 t2 -> add (interp_term t1 b) (interp_term t2 b)
2037   end
2039 function interp_fmla (f:fmla) (l:int) (b:int->value) : bool =
2040   match f with
2041   | Foo t -> foo (interp_term t b)
2042   | Forall f -> forall v. interp_fmla f (l-1) b[l <- v]
2043   end
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 }
2050 = false
2053 module TestFmla
2055 use Fmla
2057 meta "compute_max_steps" 0x10000
2059 goal g:
2060      forall a: value.
2061      ((forall x. forall y. foo (add x (add (add a dummy) y))) = True)