LAProof.C.verif_cholesky

Require Import VST.floyd.proofauto.
From vcfloat Require Import FPStdLib FPStdCompCert.
Require Import VSTlib.spec_math.
From LAProof.C Require Import floatlib cholesky.
Import FPCore FPCompCert.

#[export] Instance CompSpecs : compspecs.
Definition Vprog : varspecs.
Definition Gprog: funspecs := [sqrt_spec].


Open Scope logic.

Section WithSTD.
Context {NAN: FPCore.Nans} {t : type} {STD: is_standard t}.

Definition neg_zero := ftype_of_float (Binary.B754_zero (fprec t) (femax t) true).

Definition sumF := fold_right BPLUS neg_zero.

Fixpoint iota (i j : nat) :=
match j with O ⇒ nil | S j' ⇒ i :: iota (S i) j' end.

Definition subtract_loop A R (i j k: Z) :=
  fold_left BMINUS
    (map (fun k' ⇒ BMULT (R k' i) (R k' j)) (map Z.of_nat (iota 0 (Z.to_nat k))))
     (A i j).

Definition sum_upto (n: Z) (f: Z → ftype t) :=
  sumF (map (fun k ⇒ f (Z.of_nat k)) (iota 0 (Z.to_nat n))).

Definition cholesky_jik_ij (n: Z) (A R: Z → Z → ftype t) (i j: Z) : Prop :=
   (0 ≤ j < n) →
     (0 ≤ i < j → R i j = BDIV (subtract_loop A R i j i) (R i i))
   ∧ (i=j → R i j = BSQRT _ _ (subtract_loop A R i j i)).

Definition cholesky_jik_spec (n: Z) (A R: Z → Z → ftype t) : Prop :=
  ∀ i j, cholesky_jik_ij n A R i j.

Definition cholesky_jik_upto (n imax jmax : Z) (A R : Z → Z → ftype t) : Prop :=
  ∀ i j,
    (j<jmax → cholesky_jik_ij n A R i j)
   ∧ (j=jmax → i<imax → cholesky_jik_ij n A R i j).

Definition done_to_ij (n i j: Z) (R: Z → Z → ftype Tdouble) (i' j': Z) : val :=
  if zlt i' 0 then Vundef
  else if zlt j' 0 then Vundef
  else if zlt j' i' then Vundef
  else if zlt (j'+1) j
         then Vfloat (float_of_ftype (R i' j'))
  else if zeq (j'+1) j
       then if zlt i' i
           then Vfloat (float_of_ftype (R i' j'))
           else Vundef
  else Vundef.

Definition done_to_n (n: Z) := done_to_ij n n n.

End WithSTD.

Definition N : Z := proj1_sig (opaque_constant 1000).
Definition N_eq : N=1000 := proj2_sig (opaque_constant _).
Hint Rewrite N_eq : rep_lia.

Definition matrix := tarray (tarray tdouble N) N.

Definition list_of_fun (n: Z) (f: Z → val) : list val :=
 map (fun j ⇒ f (Z.of_nat j)) (iota 0 (Z.to_nat n)).

Definition lists_of_fun (n: Z) (f: Z → Z → val) :=
 map (fun i ⇒ list_of_fun n (f (Z.of_nat i))) (iota 0 (Z.to_nat n)).

Definition cholesky_spec {NAN: Nans}:=
 DECLARE _cholesky
 WITH rsh: share, sh: share, n: Z, A: (Z → Z → ftype Tdouble), a: val, r: val
 PRE [ tuint , tptr (tarray tdouble N), tptr (tarray tdouble N) ]
    PROP (readable_share rsh; writable_share sh; 0 ≤ n ≤ N)
    PARAMS ( Vint (Int.repr n); a; r)
    SEP (data_at rsh matrix (lists_of_fun N (done_to_n n A)) a;
         data_at_ sh matrix r)
 POST [ tvoid ]
   EX R: Z → Z → ftype Tdouble,
    PROP (cholesky_jik_spec n A R)
    RETURN ()
    SEP (data_at rsh matrix (lists_of_fun N (done_to_n n A)) a;
         data_at sh matrix (lists_of_fun N (done_to_n n R)) r).

Lemma Znth_i_iota:
  ∀ lo i hi,
   0 ≤ i < Z.of_nat hi → Znth i (iota lo hi) = (lo+Z.to_nat i)%nat.

Lemma Znth_i_list_of_fun:
  ∀ d f i n, 0 ≤ i < n → @Znth _ d i (list_of_fun n f) = f i.

Lemma Zlength_iota:
  ∀ lo n, Zlength (iota lo n) = Z.of_nat n.

Lemma length_iota:
  ∀ lo n, length (iota lo n) = n.

Lemma Znth_lists_done: ∀ N n A d d' i j imax jmax,
   n ≤ N →
   imax ≤ n → jmax ≤ n →
   0 ≤ i →
   0 ≤ j < jmax →
   i ≤ j →
   (j+1=jmax → i<imax) →
   @Znth _ d j (@Znth _ d' i (lists_of_fun N (done_to_ij n imax jmax A))) =
   Vfloat (A i j).

Lemma iota_plus1:
  ∀ lo n, iota lo (n + 1)%nat = iota lo n ++ [(lo+n)%nat].

Definition updij {T} (R: Z → Z → T) i j x :=
  fun i' j' ⇒ if zeq i i' then if zeq j j' then x else R i' j' else R i' j'.

Lemma upto_iota:
 ∀ n, upto n= map Z.of_nat (iota O n).

Lemma iota_range: ∀ k n, In k (iota 0 n) → (k<n)%nat.

Lemma upd_Znth_lists_of_fun:
  ∀ d N n R i j x,
   0 ≤ i ≤ j → j < N →
  
   upd_Znth i (lists_of_fun N (done_to_ij n i (j + 1) R))
     (upd_Znth j (@Znth _ d i (lists_of_fun N (done_to_ij n i (j + 1) R)))
        (Vfloat x))
    = lists_of_fun N (done_to_ij n (i+1) (j+1) (updij R i j x)).

Lemma update_i_lt_j:
  ∀ {t: type} {STD: is_standard t} n i j (A R: Z → Z → ftype t),
   0 ≤ i < j → j < n →
   cholesky_jik_upto n i j A R →
   let rij := BDIV (subtract_loop A R i j i) (R i i) in
    cholesky_jik_upto n (i + 1) j A (updij R i j rij).

Lemma subtract_another:
  ∀
   {t: type} {STD: is_standard t} i j k (A R: Z → Z → ftype t),
    0 ≤ i ≤ j → 0 ≤ k < j →
    subtract_loop A R i j (k+1) =
     BMINUS (subtract_loop A R i j k) (BMULT (R k i) (R k j)).

Lemma body_cholesky :
  semax_body Vprog Gprog f_cholesky cholesky_spec.