LAProof.C.verif_cholesky

Require Import VST.floyd.proofauto.
Require Import vcfloat.FPStdLib.
Require Import vcfloat.FPStdCompCert.
Require Import VSTlib.spec_math.
Require Import LAProof.C.floatlib.
Require Import LAProof.C.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 Onil | 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 kf (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 jf (Z.of_nat j)) (iota 0 (Z.to_nat n)).

Definition lists_of_fun (n: Z) (f: Z Z val) :=
 map (fun ilist_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.