(****************************************************************************)
(*                 The Calculus of Inductive Constructions                  *)
(*                                                                          *)
(*                                Projet Coq                                *)
(*                                                                          *)
(*                     INRIA        LRI-CNRS        ENS-CNRS                *)
(*              Rocquencourt         Orsay          Lyon                    *)
(*                                                                          *)
(*                                 Coq V6.2                                 *)
(*                               May 1st 1998                               *)
(*                                                                          *)
(****************************************************************************)
(* Certification of Imperative Programs                                     *)
(* Jean-Christophe Filliātre                                                *)
(****************************************************************************)
(*                                  exp_int.v                               *)
(****************************************************************************)

(* Efficient computation of X^n using
 * 
 *    X^(2n)   =     (X^n) ^ 2
 *    X^(2n+1) = X . (X^n) ^ 2
 *
 * Proofs of both fonctional and imperative programs.
 *)

Require Zpower.
Require Zcomplements.

Require Programs.
Require Ring.
Require Omega.

Definition Zdouble := [n:Z]`2*n`.

Definition Zsquare := [n:Z](Zmult n n).

(* Some auxiliary lemmas about Zdiv2 are necessary *)

Lemma Zdiv2_ge_0 : (x:Z) `x >= 0` -> `(Zdiv2 x) >= 0`.
Proof.
Destruct x; Auto.
(Destruct p; Auto).
Simpl. Omega.
Intros. (Absurd `(NEG p) >= 0`; Red; Auto).
Save.

Lemma Zdiv2_lt : (x:Z) `x > 0` -> `(Zdiv2 x) < x`.
Proof.
Destruct x.
Intro. Absurd `0 > 0`; [ Omega | Assumption ].
(Destruct p; Auto).

Simpl.
Intro p.
Replace (POS (xI p0)) with `2*(POS p0)+1`.
Omega.
Simpl. Auto.

Intro p.
Simpl.
Replace (POS (xO p0)) with `2*(POS p0)`.
Omega.
Simpl. Auto.

Simpl. Omega.

Intros. 
(Absurd `(NEG p) > 0`; Red; Auto).
(Elim p; Auto).
Omega.
Save.

(* A property of Zpower:  x^(2*n) = (x^2)^n *)

Lemma Zpower_2n : 
  (x,n:Z)`n >= 0` -> (Zpower x (Zdouble n))=(Zpower (Zsquare x) n).
Proof.
Unfold Zdouble.
Intros x n Hn.
Replace `2*n` with `n+n`.
Rewrite Zpower_exp.
Pattern n.
Apply natlike_ind.

Simpl. Auto.

Intros.
Unfold Zs.
Rewrite Zpower_exp.
Rewrite Zpower_exp.
Replace (Zpower x `1`) with x.
Replace (Zpower (Zsquare x) `1`) with (Zsquare x).
Rewrite <- H0.
Unfold Zsquare.
Ring.

Unfold Zpower; Unfold Zpower_pos; Simpl. Omega.

Unfold Zpower; Unfold Zpower_pos; Simpl. Omega.

Omega.
Omega.
Omega.
Omega.
Omega.
Assumption.
Assumption.
Omega.
Save.


(* The program *)

Correctness i_exp
  fun (x:Z)(n:Z) ->
    { `n >= 0` }
   (let y = ref 1 in
    let m = ref x in
    let e = ref n in
    begin
      while !e > 0 do
        { invariant (Zpower x n)=(Zmult y (Zpower m e)) /\ `e>=0` as I
          variant e }
        (if not (Zeven_odd_dec !e) then y := (Zmult !y !m))
          { (Zpower x n) = (Zmult y (Zpower m (Zdouble (Zdiv2 e)))) };
        m := (Zsquare !m);
        e := (Zdiv2 !e)
      done;
      !y
    end)
    { result=(Zpower x n) }
.
Proof.
(* Zodd *)
Decompose [and] I.
Rewrite (Zodd_div2 e0 H1 Test1) in H0. Rewrite H0.
Rewrite Zpower_exp.
Unfold Zdouble.
Replace (Zpower m0 `1`) with m0.
Ring.
(Unfold Zpower; Unfold Zpower_pos; Simpl; Ring).
(Generalize (Zdiv2_ge_0 e0); Omega).
Omega.
(* Zeven *)
Decompose [and] I.
Rewrite (Zeven_div2 e0 Test1) in H0. Rewrite H0.
Auto.
(* Zwf /\ I *)
Split.
(* Zwf *)
Unfold Zwf.
Repeat Split.
Generalize (Zdiv2_ge_0 e0); Omega.
Omega.
Exact (Zdiv2_lt e0 Test2).
(* invariant *)
Split.
Rewrite Post3. Unfold Zdouble. Unfold Zsquare.
Rewrite (Zpower_2n).
Trivial.
Generalize (Zdiv2_ge_0 e0); Omega.
Generalize (Zdiv2_ge_0 e0); Omega.
Split; [ Ring | Assumption ].
(* exit fo loop *)
Decompose [and] I.
Cut `e0 = 0`. Intro.
Rewrite H. Rewrite H0.
Simpl; Ring.
Omega.
Save.


(* Recursive version. *)

Correctness r_exp
  let rec exp (x:Z) (n:Z) : Z { variant n } =
    { `n >= 0` }
    (if n = 0 then
       1
     else
       let y = (exp x (Zdiv2 n)) in
       (if (Zeven_odd_dec n) then
       	  (Zmult y y)
        else
          (Zmult x (Zmult y y))) { result=(Zpower x n) }
    ) 
    { result=(Zpower x n) }
.
Proof.
Rewrite Test2. Auto.
(* w.f. *)
Unfold Zwf.
Repeat Split.
Generalize (Zdiv2_ge_0 n0) ; Omega.
Omega.
Generalize (Zdiv2_lt n0) ; Omega.
(* rec. call *)
Generalize (Zdiv2_ge_0 n0) ; Omega.
(* invariant: case even *)
Generalize (Zeven_div2 n0 Test1).
Intro Heq. Rewrite Heq.
Rewrite Post5.
Replace `2*(Zdiv2 n0)` with `(Zdiv2 n0)+(Zdiv2 n0)`.
Rewrite Zpower_exp.
Auto.
Generalize (Zdiv2_ge_0 n0) ; Omega.
Generalize (Zdiv2_ge_0 n0) ; Omega.
Omega.
(* invariant: cas odd *)
Generalize (Zodd_div2 n0 Pre1 Test1).
Intro Heq. Rewrite Heq.
Rewrite Post5.
Rewrite Zpower_exp.
Replace `2*(Zdiv2 n0)` with `(Zdiv2 n0)+(Zdiv2 n0)`.
Rewrite Zpower_exp.
Replace `(Zpower x0 1)` with x0.
Ring.
Unfold Zpower; Unfold Zpower_pos; Simpl. Omega.
Generalize (Zdiv2_ge_0 n0) ; Omega.
Generalize (Zdiv2_ge_0 n0) ; Omega.
Omega.
Generalize (Zdiv2_ge_0 n0) ; Omega.
Omega.
Save.
