Module FormalML.utils.RiemannAdd

Require Import Reals.Rbase Coq.Reals.RList.
Require Import Reals.Rfunctions.
Require Import Ranalysis_reg.
Require Import Reals.Integration.
Require Import Rtrigo_def.
Require Import List Permutation.
Require Import Sorted.

Require Import Lra Lia.
Require Import LibUtils Sums ListAdd RealAdd.

Local Open Scope R_scope.
Implicit Type f : R -> R.

Definition Partition (a b : R) (n : nat) : (list R) :=
  let inc := (b - a)/(INR n) in
  map (fun i => a + (INR i)*inc) (seq 0 (n+1)).


Lemma Partition_hd a b n d : hd d (Partition a b n) = a.

Lemma Partition_nnil a b n : Partition a b n <> nil.

Lemma Partition_nth a b n d nn :
  (nn <= n)%nat ->
  nth nn (Partition a b n) d = a + (INR nn)*(b - a)/(INR n).

Lemma Partition_length a b n : length (Partition a b n) = S n.

Lemma Partition_last a b n d :
  (0 < n)%nat ->
  last (Partition a b n) d = b.

Lemma list_hd_last {A} (l:list A) d1 d2 :
  (length l > 1)%nat ->
  l = hd d1 l :: tl (removelast l) ++ (last l d2::nil).
  
Lemma Partition_unfold_hd_tl (a b : R) (n : nat) :
  (n > 0)%nat ->
  Partition a b n =
  a::
  let inc := (b - a)/(INR n) in
  (map (fun i => a + (INR i)*inc) (seq 1 (n-1)))
    ++ (b::nil).

Lemma Partition_func_shift_nonneg a b n i:
  a <= b ->
  (0 < n)%nat ->
  0 <= INR i * (b - a) / INR n.
  
Lemma Partition_func_increasing a b n:
  a < b ->
  (0 < n)%nat ->
  Morphisms.Proper (Morphisms.respectful lt Rlt) (fun i : nat => a + INR i * ((b - a) / INR n)).

Lemma Partition_func_eq a n i :
  (0 < n)%nat ->
  a + INR i * ((a - a) / INR n) = a.

Lemma Partition_func_nondecreasing a b n:
  a <= b ->
  (0 < n)%nat ->
  Morphisms.Proper (Morphisms.respectful lt Rle) (fun i : nat => a + INR i * ((b - a) / INR n)).

Lemma Partition_StronglySorted_lt a b n :
  a < b ->
  (0 < n)%nat ->
  StronglySorted Rlt (Partition a b n).

Lemma Partition_StronglySorted_le a b n :
  a <= b ->
  (0 < n)%nat ->
  StronglySorted Rle (Partition a b n).
                                       
Lemma Partition_nth_le a b n idx1 idx2 d1 d2:
  a <= b ->
  (0 < n)%nat ->
  (idx2 <= n)%nat ->
  (idx1 <= idx2)%nat ->
  nth idx1 (Partition a b n) d1 <= nth idx2 (Partition a b n) d2.

Lemma Partition_lower_bound a b n idx :
  (a <= b) ->
  (0 < n)%nat ->
  (idx <= n)%nat ->
  a <= nth idx (Partition a b n) 0.

Lemma Partition_upper_bound a b n idx :
  (a <= b) ->
  (0 < n)%nat ->
  (idx <= n)%nat ->
  nth idx (Partition a b n) 0 <= b.

Lemma Partition_telescope f (a b : R) (n : nat) :
  (n > 0)%nat ->
  let pl := map f (Partition a b n) in
  (fold_right Rplus 0 (removelast pl)) - (fold_right Rplus 0 (tl pl)) =
  f(a) - f(b).

Lemma orderedPartition (a b : R) (n:nat) :
  (n>0)%nat ->
  a <= b ->
  let rpl := (Partition a b n) in
           pos_Rl rpl 0 = a /\
           pos_Rl rpl (pred (length rpl)) = b /\
           ordered_Rlist rpl.

Lemma find_bucket_Partition a b n idx d1 d2 needle:
  (n > 0)%nat ->
  (idx < n)%nat ->
  a <= b ->
  nth idx (Partition a b n) d1 < needle <= nth (S idx) (Partition a b n) d2 ->
  find_bucket Rle_dec needle (Partition a b n) = Some (nth idx (Partition a b n) d1, nth (S idx) (Partition a b n) d2).
    
Definition find_pt_le f a b n needle : R
  := match find_bucket Rle_dec needle (Partition a b n) with
     | Some (lower,upper) => f upper
     | None => 0
     end.

Lemma part2step (f:R -> R) (a b:R) (n : nat) :
  (n > 0)%nat ->
  a <= b ->
  IsStepFun (find_pt_le f a b n) a b.

Lemma Partition_f_increasing (f : R -> R) (a b x : R) (idx n : nat) :
  (0 < n)%nat ->
  (idx < n)%nat ->
  a <= b ->
  interval_increasing f a b ->
  nth idx (Partition a b n) 0 <= x <= nth (S idx) (Partition a b n) 0 ->
  R_dist (f x) (f (nth (S idx) (Partition a b n) 0)) <= R_dist (f (nth idx (Partition a b n) 0)) (f (nth (S idx) (Partition a b n) 0)).


Lemma Partition_f_decreasing (f : R -> R) (a b x : R) (idx n : nat) :
  (0 < n)%nat ->
  (idx < n)%nat ->
  a <= b ->
  interval_decreasing f a b ->
  nth idx (Partition a b n) 0 <= x <= nth (S idx) (Partition a b n) 0 ->
  R_dist (f x) (f (nth (S idx) (Partition a b n) 0)) <= R_dist (f (nth idx (Partition a b n) 0)) (f (nth (S idx) (Partition a b n) 0)).

Definition find_pt_le_psi f a b n needle : R
  := match find_bucket Rle_dec needle (Partition a b n) with
     | Some (lower,upper) => f lower - f upper
     | None => 0
     end.

Definition list_map_diffs f (l : list R) : (list R) :=
  map (fun '(x, y) => x-y) (adjacent_pairs (map f l)).

Lemma list_map_diffs_length f l : length (list_map_diffs f l) = pred (length l).

Lemma list_map_diff_nth_in n f (l:list R) d d1 d2 :
  (S n < length l)%nat ->
  nth n (list_map_diffs f l) d = (f (nth n l d1) - f (nth (S n) l d2)).

Lemma part2step_psi (f:R -> R) (a b:R) (n : nat) :
  (n > 0)%nat ->
  a <= b ->
  IsStepFun (find_pt_le_psi f a b n) a b.

Lemma adjacent_pairs_Partition a b n :
  adjacent_pairs (Partition a b n) =
  let inc := (b - a)/(INR n) in
  map (fun i => (a + (INR i)*inc, a + (INR (i+1))*inc)) (seq 0 n).

Lemma RiemannInt_SF_psi (f: R -> R) (a b:R) (n: nat) :
  forall (npos: (n > 0)%nat) (aleb: (a <= b)),
    RiemannInt_SF (mkStepFun (part2step_psi f a b n npos aleb)) = (f(a)-f(b))*(b-a)/(INR n).


Lemma find_bucket_Partition_exists a b n needle:
  (n > 0)%nat ->
  a <= needle <= b ->
  exists lower upper,
    find_bucket Rle_dec needle (Partition a b n) = Some (lower, upper).
            
Lemma StepBounded (f : R -> R) (a b : R) (n : nat) :
  forall (npos: (n > 0)%nat) (aleb: (a <= b)),
    interval_decreasing f a b ->
    let phi := mkStepFun (part2step f a b n npos aleb) in
    let psi := mkStepFun (part2step_psi f a b n npos aleb) in
    (forall t:R,
        a <= t <= b -> Rabs (f t - phi t) <= psi t).

Lemma natp1gz (n : nat) : (n+1 > 0)%nat.

Lemma INR_up_over_cancel r (epsilon:posreal) :
  r <> 0 ->
  Rabs (r / INR (Z.to_nat (up (Rabs r / epsilon)) + 1)) < epsilon.
                                             
Lemma RiemannInt_SF_psi_limit (f: R -> R) (a b:R) :
 forall (aleb: (a <= b)) (epsilon : posreal),
   {n:nat |
    Rabs (RiemannInt_SF (mkStepFun (part2step_psi f a b (n+1)%nat (natp1gz n) aleb))) < epsilon}.

Theorem RiemannInt_decreasing (f: R -> R) (a b:R) :
  a <= b ->
  interval_decreasing f a b ->
  Riemann_integrable f a b.

Theorem RiemannInt_increasing (f: R -> R) (a b:R) :
  a <= b ->
  interval_increasing f a b ->
  Riemann_integrable f a b.