Commit 993d4916 by Jean-Christophe Filliâtre

### bresenham: simplified spec and proof

parent d7473542
 ... @@ -17,39 +17,40 @@ module M ... @@ -17,39 +17,40 @@ module M constant y2: int constant y2: int axiom first_octant: 0 <= y2 <= x2 axiom first_octant: 0 <= y2 <= x2 (* The code. (* [best x y] expresses that the point [(x,y)] is the best [(best x y)] expresses that the point [(x,y)] is the best possible point i.e. the closest to the real line possible point i.e. the closest to the real line (see the Coq file). i.e. for all y', we have |y - x*y2/x2| <= |y' - x*y2/x2| The invariant relates [x], [y], and [e] and We stay in type [int] by multiplying everything by [x2]. *) gives lower and upper bound for [e] (see the Coq file). *) use import int.Abs use import int.Abs predicate best (x y: int) = predicate best (x y: int) = forall y': int. abs (x2 * y - x * y2) <= abs (x2 * y' - x * y2) forall y': int. abs (x2 * y - x * y2) <= abs (x2 * y' - x * y2) predicate invariant_ (x y e: int) = (** Key lemma for Bresenham's proof: if [b] is at distance less or equal e = 2 * (x + 1) * y2 - (2 * y + 1) * x2 /\ than [1/2] from the rational [c/a], then it is the closest such integer. 2 * (y2 - x2) <= e <= 2 * y2 We express this property using integers by multiplying everything by [2a]. *) lemma invariant_is_ok: forall x y e: int. invariant_ x y e -> best x y lemma closest : forall a b c: int. 0 < a -> abs (2 * a * b - 2 * c) <= a -> forall b': int. abs (a * b - c) <= abs (a * b' - c) let bresenham () = let bresenham () = let x = ref 0 in let y = ref 0 in let y = ref 0 in let e = ref (2 * y2 - x2) in let e = ref (2 * y2 - x2) in while !x <= x2 do for x = 0 to x2 do invariant { 0 <= !x <= x2 + 1 /\ invariant_ !x !y !e } invariant { !e = 2 * (x + 1) * y2 - (2 * !y + 1) * x2 } variant { x2 + 1 - !x } invariant { 2 * (y2 - x2) <= !e <= 2 * y2 } (* here we would plot (x, y) *) (* here we would plot (x, y), assert { best !x !y }; so we assert this is the best possible row y for column x *) assert { best x !y }; if !e < 0 then if !e < 0 then e := !e + 2 * y2 e := !e + 2 * y2 else begin else begin y := !y + 1; y := !y + 1; e := !e + 2 * (y2 - x2) e := !e + 2 * (y2 - x2) end; end x := !x + 1 done done end end
 (* This file is generated by Why3's Coq driver *) (* This file is generated by Why3's Coq 8.4 driver *) (* Beware! Only edit allowed sections below *) (* Beware! Only edit allowed sections below *) Require Import ZArith. Require Import BuiltIn. Require Import Rbase. Require BuiltIn. Definition unit := unit. Require int.Int. Require int.Abs. Parameter mark : Type. (* Why3 assumption *) Definition unit := unit. Parameter at1: forall (a:Type), a -> mark -> a. (* Why3 assumption *) Inductive ref (a:Type) {a_WT:WhyType a} := Implicit Arguments at1. Parameter old: forall (a:Type), a -> a. Implicit Arguments old. Inductive ref (a:Type) := | mk_ref : a -> ref a. | mk_ref : a -> ref a. Implicit Arguments mk_ref. Axiom ref_WhyType : forall (a:Type) {a_WT:WhyType a}, WhyType (ref a). Existing Instance ref_WhyType. Definition contents (a:Type)(u:(ref a)): a := Implicit Arguments mk_ref [[a] [a_WT]]. match u with | mk_ref contents1 => contents1 (* Why3 assumption *) Definition contents {a:Type} {a_WT:WhyType a} (v:(@ref a a_WT)): a := match v with | (mk_ref x) => x end. end. Implicit Arguments contents. Parameter x2: Z. Parameter y2: Z. Axiom first_octant : (0%Z <= (y2 ))%Z /\ ((y2 ) <= (x2 ))%Z. Parameter x2: Z. Axiom Abs_pos : forall (x:Z), (0%Z <= (Zabs x))%Z. Parameter y2: Z. Definition best(x:Z) (y:Z): Prop := forall (yqt:Z), Axiom first_octant : (0%Z <= y2)%Z /\ (y2 <= x2)%Z. ((Zabs (((x2 ) * y)%Z - (x * (y2 ))%Z)%Z) <= (Zabs (((x2 ) * yqt)%Z - (x * (y2 ))%Z)%Z))%Z. Definition invariant_(x:Z) (y:Z) (e:Z): Prop := (* Why3 assumption *) (e = (((2%Z * (x + 1%Z)%Z)%Z * (y2 ))%Z - (((2%Z * y)%Z + 1%Z)%Z * (x2 ))%Z)%Z) /\ Definition best (x:Z) (y:Z): Prop := forall (y':Z), (((2%Z * ((y2 ) - (x2 ))%Z)%Z <= e)%Z /\ (e <= (2%Z * (y2 ))%Z)%Z). ((Zabs ((x2 * y)%Z - (x * y2)%Z)%Z) <= (Zabs ((x2 * y')%Z - (x * y2)%Z)%Z))%Z. (* YOU MAY EDIT THE CONTEXT BELOW *) (*s First a tactic [Case_Zabs] to do case split over [(Zabs x)]: (*s First a tactic [Case_Zabs] to do case split over [(Zabs x)]: introduces two subgoals, one where [x] is assumed to be non negative introduces two subgoals, one where [x] is assumed to be non negative and thus where [Zabs x] is replaced by [x]; and another where and thus where [Zabs x] is replaced by [x]; and another where ... @@ -86,18 +75,17 @@ Ltac ZCompare x y H := ... @@ -86,18 +75,17 @@ Ltac ZCompare x y H := Ltac RingSimpl x y := replace x with y; [ idtac | ring ]. Ltac RingSimpl x y := replace x with y; [ idtac | ring ]. (*s Key lemma for Bresenham's proof: if [b] is at distance less or equal Require Import Why3. than [1/2] from the rational [c/a], then it is the closest such integer. Ltac ae := why3 "Alt-Ergo,0.95.1," timelimit 3. We express this property in [Z], thus multiplying everything by [2a]. *) Lemma closest : (* Why3 goal *) forall a b c:Z, Theorem closest : forall (a:Z) (b:Z) (c:Z), (0%Z < a)%Z -> (0 <= a)%Z -> (((Zabs (((2%Z * a)%Z * b)%Z - (2%Z * c)%Z)%Z) <= a)%Z -> forall (b':Z), (Zabs (2 * a * b - 2 * c) <= a)%Z -> ((Zabs ((a * b)%Z - c)%Z) <= (Zabs ((a * b')%Z - c)%Z))%Z). forall b':Z, (Zabs (a * b - c) <= Zabs (a * b' - c))%Z. (* Why3 intros a b c h1 h2 b'. *) Proof. intros a b c Ha Hmin. intros a b c Ha Hmin. generalize (proj2 (Zabs_le (2 * a * b - 2 * c) a Ha) Hmin). assert (Ha': (0 <= a)%Z) by omega. generalize (proj2 (Zabs_le (2 * a * b - 2 * c) a Ha') Hmin). intros Hmin' b'. intros Hmin' b'. elim (Z_le_gt_dec (2 * a * b) (2 * c)); intro Habc. elim (Z_le_gt_dec (2 * a * b) (2 * c)); intro Habc. (* 2ab <= 2c *) (* 2ab <= 2c *) ... @@ -105,22 +93,8 @@ rewrite (Zabs_non_eq (a * b - c)). ... @@ -105,22 +93,8 @@ rewrite (Zabs_non_eq (a * b - c)). ZCompare b b' Hbb'. ZCompare b b' Hbb'. (* b > b' *) (* b > b' *) rewrite (Zabs_non_eq (a * b' - c)). rewrite (Zabs_non_eq (a * b' - c)). apply Zle_left_rev. ae. RingSimpl (Zopp (a * b' - c) + Zopp (Zopp (a * b - c)))%Z ae. (a * (b - b'))%Z. apply Zmult_le_0_compat; omega. apply Zge_le. apply Zge_trans with (m := (a * b - c)%Z). apply Zmult_ge_reg_r with (p := 2%Z). omega. RingSimpl (0 * 2)%Z 0%Z. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. RingSimpl (a * b' - c)%Z (a * b' + Zopp c)%Z. RingSimpl (a * b - c)%Z (a * b + Zopp c)%Z. apply Zle_ge. apply Zplus_le_compat_r. apply Zmult_le_compat_l; omega. (* b < b' *) (* b < b' *) rewrite (Zabs_eq (a * b' - c)). rewrite (Zabs_eq (a * b' - c)). apply Zmult_le_reg_r with (p := 2%Z). apply Zmult_le_reg_r with (p := 2%Z). ... @@ -135,7 +109,7 @@ ZCompare b b' Hbb'. ... @@ -135,7 +109,7 @@ ZCompare b b' Hbb'. apply Zplus_le_compat. apply Zplus_le_compat. RingSimpl (2 * a)%Z (2 * a * 1)%Z. RingSimpl (2 * a)%Z (2 * a * 1)%Z. RingSimpl (2 * (a * b' - a * b))%Z (2 * a * (b' - b))%Z. RingSimpl (2 * (a * b' - a * b))%Z (2 * a * (b' - b))%Z. apply Zmult_le_compat_l; omega. ae. RingSimpl (2 * (a * b - c))%Z (2 * a * b - 2 * c)%Z. RingSimpl (2 * (a * b - c))%Z (2 * a * b - 2 * c)%Z. omega. omega. (* 0 <= ab'-c *) (* 0 <= ab'-c *) ... @@ -144,29 +118,14 @@ ZCompare b b' Hbb'. ... @@ -144,29 +118,14 @@ ZCompare b b' Hbb'. apply Zplus_le_compat. apply Zplus_le_compat. RingSimpl a (a * 1)%Z. RingSimpl a (a * 1)%Z. RingSimpl (a * 1 * b' - a * 1 * b)%Z (a * (b' - b))%Z. RingSimpl (a * 1 * b' - a * 1 * b)%Z (a * (b' - b))%Z. apply Zmult_le_compat_l; omega. ae. apply Zmult_le_reg_r with (p := 2%Z). omega. apply Zle_trans with (Zopp a). apply Zle_trans with (Zopp a). omega. omega. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. ae. (* b = b' *) (* b = b' *) rewrite <- Hbb'. ae. rewrite (Zabs_non_eq (a * b - c)). ae. omega. apply Zge_le. apply Zmult_ge_reg_r with (p := 2%Z). omega. RingSimpl (0 * 2)%Z 0%Z. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. apply Zge_le. apply Zmult_ge_reg_r with (p := 2%Z). omega. RingSimpl (0 * 2)%Z 0%Z. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. (* 2ab > 2c *) (* 2ab > 2c *) rewrite (Zabs_eq (a * b - c)). rewrite (Zabs_eq (a * b - c)). ... @@ -178,89 +137,32 @@ ZCompare b b' Hbb'. ... @@ -178,89 +137,32 @@ ZCompare b b' Hbb'. RingSimpl (Zopp (a * b' - c) * 2)%Z RingSimpl (Zopp (a * b' - c) * 2)%Z (2 * (c - a * b) + 2 * (a * b - a * b'))%Z. (2 * (c - a * b) + 2 * (a * b - a * b'))%Z. apply Zle_trans with a. apply Zle_trans with a. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. ae. omega. apply Zle_trans with (Zopp a + 2 * a)%Z. apply Zle_trans with (Zopp a + 2 * a)%Z. omega. omega. apply Zplus_le_compat. apply Zplus_le_compat. RingSimpl (2 * (c - a * b))%Z (2 * c - 2 * a * b)%Z. ae. omega. RingSimpl (2 * a)%Z (2 * a * 1)%Z. RingSimpl (2 * a)%Z (2 * a * 1)%Z. RingSimpl (2 * (a * b - a * b'))%Z (2 * a * (b - b'))%Z. RingSimpl (2 * (a * b - a * b'))%Z (2 * a * (b - b'))%Z. apply Zmult_le_compat_l; omega. ae. (* 0 >= ab'-c *) (* 0 >= ab'-c *) RingSimpl (a * b' - c)%Z (a * b' - a * b + (a * b - c))%Z. RingSimpl (a * b' - c)%Z (a * b' - a * b + (a * b - c))%Z. RingSimpl 0%Z (Zopp a + a)%Z. RingSimpl 0%Z (Zopp a + a)%Z. apply Zplus_le_compat. apply Zplus_le_compat. RingSimpl (Zopp a) (a * (-1))%Z. RingSimpl (Zopp a) (a * (-1))%Z. RingSimpl (a * b' - a * b)%Z (a * (b' - b))%Z. RingSimpl (a * b' - a * b)%Z (a * (b' - b))%Z. apply Zmult_le_compat_l; omega. ae. apply Zmult_le_reg_r with (p := 2%Z). ae. omega. apply Zle_trans with a. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. omega. (* b < b' *) (* b < b' *) rewrite (Zabs_eq (a * b' - c)). rewrite (Zabs_eq (a * b' - c)). apply Zle_left_rev. apply Zle_left_rev. RingSimpl (a * b' - c + Zopp (a * b - c))%Z (a * (b' - b))%Z. RingSimpl (a * b' - c + Zopp (a * b - c))%Z (a * (b' - b))%Z. apply Zmult_le_0_compat; omega. ae. apply Zle_trans with (m := (a * b - c)%Z). apply Zle_trans with (m := (a * b - c)%Z). apply Zmult_le_reg_r with (p := 2%Z). ae. omega. ae. RingSimpl (0 * 2)%Z 0%Z. ae. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. ae. omega. RingSimpl (a * b' - c)%Z (a * b' + Zopp c)%Z. RingSimpl (a * b - c)%Z (a * b + Zopp c)%Z. apply Zplus_le_compat_r. apply Zmult_le_compat_l; omega. (* b = b' *) rewrite <- Hbb'. rewrite (Zabs_eq (a * b - c)). omega. apply Zmult_le_reg_r with (p := 2%Z). omega. RingSimpl (0 * 2)%Z 0%Z. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. apply Zmult_le_reg_r with (p := 2%Z). omega. RingSimpl (0 * 2)%Z 0%Z. RingSimpl ((a * b - c) * 2)%Z (2 * a * b - 2 * c)%Z. omega. Qed. (* DO NOT EDIT BELOW *) Theorem invariant_is_ok : forall (x:Z) (y:Z) (e:Z), (invariant_ x y e) -> (best x y). (* YOU MAY EDIT THE PROOF BELOW *) Proof. intros x y e. unfold invariant_; unfold best; intros [E I'] y'. cut (0 <= x2)%Z; [ intro Hx2 | idtac ]. apply closest. assumption. apply (proj1 (Zabs_le (2 * x2 * y - 2 * (x * y2)) x2 Hx2)). rewrite E in I'. split. (* 0 <= x2 *) generalize (proj2 I'). RingSimpl (2 * (x + 1) * y2 - (2 * y + 1) * x2)%Z (2 * x * y2 - 2 * x2 * y + 2 * y2 - x2)%Z. intro. RingSimpl (2 * (x * y2))%Z (2 * x * y2)%Z. omega. (* 0 <= x2 *) generalize (proj1 I'). RingSimpl (2 * (x + 1) * y2 - (2 * y + 1) * x2)%Z (2 * x * y2 - 2 * x2 * y + 2 * y2 - x2)%Z. RingSimpl (2 * (y2 - x2))%Z (2 * y2 - 2 * x2)%Z. RingSimpl (2 * (x * y2))%Z (2 * x * y2)%Z. omega. omega. Qed. Qed. (* DO NOT EDIT BELOW *)
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!