Library alpha

Require Import Dekker.
Require Import discriminant3.
Require Export Caduceus.
Require Export WhyFloatsStrictLegacy.
Require Import Wf_nat.

Hint Resolve pdGreaterThanOne pdGivesBound EvenClosestRoundedModeP RND_EvenClosest_correct RND_EvenClosest_canonic (FcanonicBound radix).

Local Coercion FtoRradix: Float.float >-> R.

Lemma Rplus_eq_compat: forall a1 b1 a2 b2:R,
    a1=a2 -> b1=b2 -> (a1+b1=a2+b2)%R.

Lemma Rabs_le_trans: forall a b c d:R,
  (Rabs (a-c)+Rabs (c-b) <= d)%R ->
  (Rabs (a-b) <= d)%R.

Section Calculaaux.

Variable x:R.
Variable f:float.

Hypothesis feq: (f= RND_EvenClosest bdouble radix 53 x)%R.

Hypothesis xGe: (powerRZ radix (-1010) <= Rabs x )%R.

Lemma RoudGeOv2: (FtoRradix f <> 0)%R
  -> (Rabs x/2 <= Rabs f)%R.

Lemma RoundGe1010: (powerRZ radix (-1010) <= Rabs f )%R.

Lemma f1010normal: Fnormal radix bdouble f.

Lemma RoundErrLe1:
  (Rabs (f-x) <= (powerRZ radix (-53) * Rabs f))%R.

Lemma RoundErrLe2:
  (Rabs (f-x) <= (powerRZ radix (-53)/(1-powerRZ radix (-53)) * Rabs x))%R.

End Calculaaux.

Section CalculaEquatlities.

Lemma powerRZ_53_eq: (1 / 9007199254740992 = powerRZ radix (-53))%R.

Lemma powerRZ_51_eq: (1 / 2251799813685248 = powerRZ radix (-51))%R.

Lemma powerRZ_50_eq: (1 / 1125899906842624 = powerRZ radix (-50))%R.

Lemma powerRZ_1000_eq: ((1 * 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376)
   = powerRZ radix 1000)%R.

End CalculaEquatlities.
Section Calcula1.

Variable dx dt v t a1:double.

Hypothesis dxPos: (0 < double_value dx)%R.
Hypothesis dtPos: (0 < double_value dt)%R.
Hypothesis vPos: (0 < double_value v)%R.

Hypothesis dxEPos: (0 < double_exact dx)%R.

Hypothesis Reldx: (Rabs (double_exact dx -double_value dx)/double_value dx <= powerRZ radix (-53))%R.
Hypothesis Reldt: (Rabs (double_exact dt -double_value dt)/double_value dt <= powerRZ radix (-51))%R.
Hypothesis vexact: double_exact v=double_value v.

Hypothesis dxLe: (double_value dx <= /2)%R.
Hypothesis dtGe: (powerRZ radix (-1011) <= double_value dt)%R.
Hypothesis ValLe: (double_value dt/double_value dx*double_value v <= 1)%R.
Hypothesis ValGe: (powerRZ radix (-503) <= double_value dt/double_value dx*double_value v)%R.

Hypothesis teq: (div_double_post nearest_even dt dx t).
Hypothesis a1eq: (mul_double_post nearest_even t v a1).

Lemma dtdxLe: (powerRZ radix (-1010) <= Rabs (double_value dt/double_value dx) )%R.

Lemma tvGe: (powerRZ radix (-504) <= Rabs (double_value t * double_value v))%R.

Lemma a1Ge0: (0 <= double_value a1)%R.

Lemma a1Ge: (powerRZ radix (-504) <= double_value a1)%R.

Lemma a1_err:
 (double_round_error a1 <= powerRZ radix (-50) - powerRZ radix (-54))%R.

End Calcula1.

Section Calcula.

Variable dx dt v t a1 a:double.

Hypothesis dxPos: (0 < double_value dx)%R.
Hypothesis dtPos: (0 < double_value dt)%R.
Hypothesis vPos: (0 < double_value v)%R.
Hypothesis dxEPos: (0 < double_exact dx)%R.
Hypothesis dtEPos: (0 < double_exact dt)%R.

Hypothesis Reldx: (Rabs (double_exact dx -double_value dx)/double_value dx <= powerRZ radix (-53))%R.
Hypothesis Reldt: (Rabs (double_exact dt -double_value dt)/double_value dt <= powerRZ radix (-51))%R.
Hypothesis vexact: double_exact v=double_value v.

Hypothesis dxLe: (double_value dx <= /2)%R.
Hypothesis dtGe: (powerRZ radix (-1000) <= double_exact dt)%R.
Hypothesis ValLe: (double_exact dt/double_exact dx*double_exact v <= 1 - powerRZ radix (-50))%R.
Hypothesis ValGe: (powerRZ radix (-500) <= double_exact dt/double_exact dx*double_exact v)%R.

Theorem dtGe2: (powerRZ radix (-1011) <= double_value dt)%R.

Lemma ValGe2: (powerRZ radix (-503) <= double_value dt/double_value dx*double_value v)%R.

Lemma ValLe2: (double_value dt/double_value dx*double_value v <= 1)%R.

Hypothesis teq: (div_double_post nearest_even dt dx t).
Hypothesis a1eq: (mul_double_post nearest_even t v a1).
Hypothesis aeq: (mul_double_post nearest_even a1 a1 a).

Lemma aGe0: (0 <= double_value a)%R.

Lemma a1Lt1: (double_value a1 < 1)%R.

Lemma aLe1: (double_value a <= 1)%R.

Lemma a_err:
 (double_round_error a <= powerRZ radix (-49))%R.

End Calcula.

Section Sum_Z_Z.

Fixpoint sum_f_z_aux (f:Z->R) (a:Z) (b:Z) (n:nat) {struct n } :R
   := match n with
        O => f a
      | S m => (f a + sum_f_z_aux f (a+1) b m)%R
    end.

Definition sum_f_z (f:Z->R) (a:Z) (b:Z) : R :=
     sum_f_z_aux f a b (Zabs_nat (b-a)).

Lemma sum_f_z_1: forall (i:Z) (f:Z->R), sum_f_z f i i = f i.

Lemma sum_f_z_down: forall (i j:Z) (f:Z->R),
   (i < j)%Z ->
   (sum_f_z f i j = f i + sum_f_z f (i+1) j)%R.

Lemma sum_f_z_up_aux: forall (n:nat) (i j:Z) (f:Z->R),
    (i < j)%Z -> (n=Zabs_nat (j -1- i)) ->
   (sum_f_z_aux f i j (S n)= f j + sum_f_z_aux f i (j-1) n)%R.

Lemma sum_f_z_up: forall (i j:Z) (f:Z->R),
    (i < j)%Z ->
   (sum_f_z f i j = f j + sum_f_z f i (j-1))%R.

Lemma sum_f_z_split_aux: forall (n:nat) (i j k:Z) (f:Z->R),
    (n=Zabs_nat (k -1- j)) ->
    (i <= j < k)%Z ->
   (sum_f_z_aux f i k (Zabs_nat (k - i))= sum_f_z_aux f i j (Zabs_nat (j-i))
         + sum_f_z_aux f (j+1) k n)%R.

Lemma sum_f_z_split: forall (i j k:Z) (f:Z->R),
    (i <= j < k)%Z ->
   (sum_f_z f i k = sum_f_z f i j+ sum_f_z f (j+1) k)%R.

Lemma sum_f_z_eq: forall (i j:Z) (f1 f2:Z->R),
    (i <= j)%Z ->
    (forall (k:Z), (i<=k)%Z -> (k<=j)%Z -> f1 k = f2 k)
    -> sum_f_z f1 i j=sum_f_z f2 i j.

Lemma sum_f_z_le: forall (i j:Z) (f1 f2:Z->R),
     (i <= j)%Z ->
    (forall (k:Z), (i<=k)%Z -> (k<=j)%Z -> (f1 k <= f2 k)%R)
    -> (sum_f_z f1 i j <= sum_f_z f2 i j)%R.

Lemma sum_f_z_zero:forall (i j:Z),
   (i <= j)%Z -> (sum_f_z (fun _ => 0) i j =0)%R.

Lemma sum_f_z_pos: forall (i j:Z) (f:Z->R),
     (i <= j)%Z ->
    (forall (k:Z), (i<=k)%Z -> (k<=j)%Z -> (0 <= f k)%R)
    -> (0 <= sum_f_z f i j)%R.

Lemma sum_f_z_cte:forall (i j:Z) (C:R),
   (i <= j)%Z -> (sum_f_z (fun _ => C) i j =(j-i+1)*C)%R.

Lemma sum_f_z_plus: forall (i j:Z) (f1 f2:Z->R),
   (i <= j)%Z ->
   (sum_f_z (fun i => (f1 i+f2 i)%R) i j=
      sum_f_z f1 i j + sum_f_z f2 i j)%R.

Lemma sum_f_z_mult: forall (i j:Z) (f:Z->R) (C:R),
   (i <= j)%Z ->
   (sum_f_z (fun i => (C*f i)%R) i j=
      C* sum_f_z f i j)%R.

Lemma sum_f_z_shift1: forall (i j:Z) (f:Z->R),
   (i <=j)%Z ->
   (sum_f_z f i j=sum_f_z (fun i=> f (i+1)%Z) (i-1) (j-1))%R.

Lemma sum_f_z_shift2: forall (i j:Z) (f:Z->R),
   (i <= j)%Z ->
   sum_f_z f i j=sum_f_z (fun i=> f (i-1)%Z) (i+1) (j+1).

Lemma sum_f_z_shift: forall (i j k:Z) (f:Z->R),
   (i <=j)%Z ->
   (sum_f_z f (i-k) (j-k)=sum_f_z (fun i=> f (i-k)%Z) i j)%R.

Lemma sum_f_z_Rabs: forall (i j:Z) (f:Z->R),
   (i <= j)%Z ->
   (Rabs (sum_f_z f i j) <= sum_f_z (fun i=> (Rabs (f i))) i j)%R.

Lemma sum_f_z_switch: forall (i j:Z) (f:Z->R),
   (i <=j)%Z ->
   (sum_f_z f i j=sum_f_z (fun k=> f (-k)%Z) (-j) (-i))%R.

End Sum_Z_Z.

Section AlphaDefinition.
Variable a:R.
Hypothesis aGt0: (0 < a)%R.

Fixpoint alpha (i:Z) (k:nat) {struct k} : R := match k with
     O => match i with 0%Z => 1%R
                                  | _ => 0%R
              end
   | S k' => match k' with
              O => match i with 0%Z => (2*(1-a))%R
                                  | (Zpos xH) => a
                                  | (Zneg xH) => a
                                  | _ => 0%R
                        end
           | (S k'') => (a*(alpha (i-1) k'+ alpha (i+1) k') +
                         2*(1-a)* alpha i k' - alpha i k'')%R
             end
     end.

Lemma alpha_rew: forall (k:nat) (i:Z),
 alpha i (S (S k)) = (a*(alpha (i-1) (S k)+ alpha (i+1) (S k)) +
                         2*(1-a)* alpha i (S k) - alpha i k)%R.

Lemma alpha_sym: forall (k:nat) ( i:Z), alpha i k = alpha (-i) k.

Lemma alpha_zero1: forall (k:nat) (i:Z),
    (i < - Z_of_nat k)%Z -> alpha i k =0%R.

Lemma alpha_zero2: forall (i:Z) (k:nat), (Z_of_nat k < i)%Z
    -> alpha i k =0%R.

Lemma alpha_ak: forall (k:nat),
   alpha (Z_of_nat k) k =Rpower a (INR k).

End AlphaDefinition.

Section alphaSum.
Variable a:R.
Hypothesis aGt0: (0 < a)%R.

Lemma alpha_sum: forall (k:nat),
   (sum_f_z (fun i => alpha a i k) (- Z_of_nat k) (Z_of_nat k)
     = INR k+1)%R.

Axiom alpha_pos: forall (k:nat) ( i:Z), (0 <= alpha a i k)%R.

End alphaSum.
Section alphaOK.

Variable a:R.
Variable eps:Z-> Z-> R.

Hypothesis aGt0: (0 < a)%R.

Definition Err := fun i k:Z =>
      (sum_f_z
         (fun l => sum_f_z
                  (fun j => (alpha a j (Zabs_nat l)
                                    * (eps (i+j)%Z (k-l)%Z))%R)
          (-l) l)
    0 k).

Lemma alpha_err_ok_aux1: forall i:Z, forall k:nat,
  (0 < k)%Z ->
  (Err i (k+1) = eps i (k+1)+
                      ( a*(eps (i-1) k+eps (i+1) k)+2*(1-a)*eps i k)
                     + sum_f_z
                      (fun l => sum_f_z
                             (fun j => (alpha a j (Zabs_nat l+1)
                                        * (eps (i+j)%Z (k-l)%Z))%R)
                            (-l-1) (l+1))
                      1 k)%R.

Lemma alpha_err_ok_aux2: forall i:Z, forall k:nat,
   (0 < k)%Z ->
  (Err i (k-1) = sum_f_z
                       (fun l => sum_f_z
                             (fun j => (alpha a j ((Zabs_nat l)-1)
                                        * (eps (i+j)%Z (k-l)%Z))%R)
                            (-l+1) (l-1))
                      1 k)%R.

Lemma alpha_err_ok_aux3: forall i:Z, forall k:nat,
  (0 < k)%Z ->
  (Err (i-1) k = eps (i-1) k +
                     sum_f_z
                       (fun l => sum_f_z
                             (fun j => (alpha a (j+1) (Zabs_nat l)
                                        * (eps (i+j)%Z (k-l)%Z))%R)
                            (-l-1) (l-1))
                      1 k)%R.

Lemma alpha_err_ok_aux4: forall i:Z, forall k:nat,
  (0 < k)%Z ->
  (Err (i+1) k = eps (i+1) k +
                     sum_f_z
                       (fun l => sum_f_z
                             (fun j => (alpha a (j-1) (Zabs_nat l)
                                        * (eps (i+j)%Z (k-l)%Z))%R)
                            (-l+1) (l+1))
                      1 k)%R.

Lemma alpha_err_ok_aux5: forall i:Z, forall k:nat,
   (0 < k)%Z ->
  (Err i k = eps i k +
                     sum_f_z
                       (fun l => sum_f_z
                             (fun j => (alpha a j (Zabs_nat l)
                                        * (eps (i+j)%Z (k-l)%Z))%R)
                            (-l) l)
                      1 k)%R.

Lemma alpha_err_ok_nat:
  forall i:Z, forall k:nat,
   (0 < k)%Z ->
  (Err i (k+1)= eps i (k+1)+a*(Err (i+1) k + Err (i-1) k)
         + 2*(1-a)*Err i k - Err i (k-1))%R.

Lemma alpha_err_ok:
  forall i k:Z,
   (0 < k)%Z ->
  (Err i (k+1)= eps i (k+1)+a*(Err (i+1) k + Err (i-1) k)
         + 2*(1-a)*Err i k - Err i (k-1))%R.

End alphaOK.