Require Import Reals Ranalysis5 Psatz ClassicalEpsilon.
Require Import Rcomplements RInt.
Require Import Coquelicot.
Require Import Morphisms.

Lemma Rlt_1_2 : 1 < 2.  Proof. psatzl R. Qed.

Ltac lt0 :=
  match goal with
  | |- _ => assumption
  | |- 0 < exp _ => apply exp_pos
  | |- 0 < pos _ => apply cond_pos
  | |- _ > 0 => unfold Rgt; lt0
  | |- 0 < 1 =>  apply Rlt_0_1
  | |- 0 < 2 => apply Rlt_0_2
  | |- 0 < PI => apply PI_RGT_0
  | |- _ <> 0 => apply Rgt_not_eq; lt0
  | |- 0 < Rabs _ + _ => apply (Rplus_le_lt_0_compat _ _ (Rabs_pos _)); lt0
  | |- 0 < _ * _ => apply Rmult_lt_0_compat; lt0
  | |- 0 < _ + _ => apply Rplus_lt_0_compat; lt0
  | |- 0 < Rmin _ _ => apply Rmin_glb_lt; lt0
  | |- 0 < / _ => apply Rinv_0_lt_compat; lt0
  | |- 0 < sqrt _ => apply sqrt_lt_R0; lt0
  | |- 0 < _ / _ => unfold Rdiv; apply Rmult_lt_0_compat; lt0
  | |- 0 < _ ^ _ => apply pow_lt; lt0
  | |- 0 < _ ^ 2 + _ => apply Rplus_le_lt_0_compat;[apply pow2_ge_0 | lt0]
  | |- _ <= _ => apply Rlt_le; lt0
  | |- _ => psatzl R
  end.

Lemma derivable_pt_lim_ext : (* This one should belong in a common file. *)
  forall f g x d, (forall y, f y = g y) ->
    derivable_pt_lim f x d -> derivable_pt_lim g x d.
intros f g x d ext Pd eps ep; destruct (Pd _ ep) as [delta P'].
exists delta; intros h hn0 hdelta; rewrite <- !ext; auto.
Qed.

Lemma ln_ub : forall k x, 1 <= x ->
  ln x <= sum_f_R0 (fun n => (-1)^n * (x - 1) ^ (n + 1)/INR(n+1)) (2 * k).
intros k x [cx | cx]; rewrite Rminus_le_0.
 destruct (MVT_cor2
             (fun y => sum_f_R0 (fun n => (-1) ^ n * (y - 1) ^ (n + 1) /INR(n+1))
                          (2 * k) - ln y)
             (fun y => sum_f_R0 (fun n => (-1) ^ n * (y - 1) ^ n) (2 * k) - / y)
             1 x) as [c [qc cc]];[assumption | | ].
  intros c cc; apply derivable_pt_lim_minus.
   induction (2 * k)%nat as [ | k' IHk'].
    simpl; apply (derivable_pt_lim_ext (fun y => y - 1)).
     intros; field.
    replace (1 * 1) with (1 - 0) by ring.
    now apply derivable_pt_lim_minus;
     [apply derivable_pt_lim_id | apply derivable_pt_lim_const].
   apply (derivable_pt_lim_ext
      (fun y : R =>
      sum_f_R0 (fun n : nat => (-1) ^ n * (y - 1) ^ (n + 1) / INR (n + 1))
        k' + (-1) ^ (S k') * (/ INR (S k' + 1) * (y - 1) ^ (S k' + 1)))).
    intros y; rewrite (Rmult_comm (/ _)), <- Rmult_assoc; reflexivity.
   replace (sum_f_R0 (fun n => (-1) ^ n * (c - 1) ^ n) (S k'))  with
      (sum_f_R0 (fun n => (-1) ^ n * (c - 1) ^ n) k' +
          (-1) ^ (S k') * (/INR (S k' + 1) * 
            ((INR (S k' + 1) * (c - 1) ^ (pred (S k' + 1))) * (1 - 0)))).
   apply derivable_pt_lim_plus;[ assumption | ].
    apply derivable_pt_lim_scal, derivable_pt_lim_scal. 
    apply (derivable_pt_lim_comp (fun y => y - 1) (fun y => y ^ (S k' + 1)) c
         (1 - 0) (INR (S k' + 1) * (c - 1) ^ pred (S k' + 1))).
     apply derivable_pt_lim_minus;[apply derivable_pt_lim_id |].
     apply derivable_pt_lim_const.
    apply derivable_pt_lim_pow.   
    rewrite Rminus_0_r, Rmult_1_r, <- (Rmult_assoc (/ _)).
    rewrite Rinv_l, Rmult_1_l, Plus.plus_comm;[reflexivity| apply not_0_INR; lia].
   apply derivable_pt_lim_ln; psatzl R.
   replace (sum_f_R0 (fun n => (-1) ^ n * (1 - 1) ^ (n + 1) / (INR (n + 1)))
              (2 * k)) with 0 in qc.
   rewrite ln_1, !Rminus_0_r in qc; rewrite qc.
  rewrite <- (sum_eq (fun n => ((-1) * (c - 1))^n)), tech3; try psatzl R.
   replace (1 - -1 * (c - 1)) with c by ring.
   rewrite Rpow_mult_distr, pow_1_odd, Ropp_mult_distr_l_reverse, Rmult_1_l.
   replace (/c) with (1 * /c) by ring.
   unfold Rdiv, Rminus at 2; rewrite Ropp_involutive.
   rewrite <- Rmult_minus_distr_r.
   assert (0 < /c) by (apply Rinv_0_lt_compat; psatzl R).
   repeat (apply Rmult_le_pos;[ | psatzl R]).
   assert (0 <= (c - 1) ^ S (2 * k)) by (apply pow_le; psatzl R); psatzl R.
  intros; rewrite Rpow_mult_distr; reflexivity.
 symmetry; apply sum_eq_R0.
 intros; replace (n + 1)%nat with (S n) by ring; simpl. 
 rewrite Rminus_eq_0; unfold Rdiv; rewrite !Rmult_0_l, Rmult_0_r, Rmult_0_l.
 reflexivity.
rewrite <- cx, sum_eq_R0, ln_1;[psatzl R | ].
intros; rewrite Plus.plus_comm, Rminus_eq_0; simpl; unfold Rdiv.
rewrite Rmult_0_l, Rmult_0_r, Rmult_0_l; reflexivity.
Qed.

Lemma ln_lb : forall k x, 1 <= x ->
  sum_f_R0 (fun n => (-1)^n * (x - 1) ^ (n + 1)/INR(n+1)) (2 * k + 1) <= ln x.
intros k x [cx | cx]; rewrite Rminus_le_0.
 destruct (MVT_cor2
             (fun y => ln y - 
                sum_f_R0 (fun n => (-1) ^ n * (y - 1) ^ (n + 1) /INR(n+1))
                          (2 * k + 1))
             (fun y => / y - sum_f_R0 (fun n => (-1) ^ n * (y - 1) ^ n) (2 * k + 1))
             1 x) as [c [qc cc]];[assumption | | ].
  intros c cc; apply derivable_pt_lim_minus.
   apply derivable_pt_lim_ln; psatzl R.
  induction (2 * k + 1)%nat as [ | k' IHk'].
   simpl; apply (derivable_pt_lim_ext (fun y => y - 1)).
    intros; field.
   replace (1 * 1) with (1 - 0) by ring.
   now apply derivable_pt_lim_minus;
     [apply derivable_pt_lim_id | apply derivable_pt_lim_const].
  apply (derivable_pt_lim_ext
      (fun y : R =>
      sum_f_R0 (fun n : nat => (-1) ^ n * (y - 1) ^ (n + 1) / INR (n + 1))
        k' + (-1) ^ (S k') * (/ INR (S k' + 1) * (y - 1) ^ (S k' + 1)))).
   intros y; rewrite (Rmult_comm (/ _)), <- Rmult_assoc; reflexivity.
  replace (sum_f_R0 (fun n => (-1) ^ n * (c - 1) ^ n) (S k'))  with
      (sum_f_R0 (fun n => (-1) ^ n * (c - 1) ^ n) k' +
          (-1) ^ (S k') * (/INR (S k' + 1) * 
            ((INR (S k' + 1) * (c - 1) ^ (pred (S k' + 1))) * (1 - 0)))).
  apply derivable_pt_lim_plus;[ assumption | ].
   apply derivable_pt_lim_scal, derivable_pt_lim_scal. 
   apply (derivable_pt_lim_comp (fun y => y - 1) (fun y => y ^ (S k' + 1)) c
         (1 - 0) (INR (S k' + 1) * (c - 1) ^ pred (S k' + 1))).
    apply derivable_pt_lim_minus;[apply derivable_pt_lim_id |].
    apply derivable_pt_lim_const.
   apply derivable_pt_lim_pow.   
   rewrite Rminus_0_r, Rmult_1_r, <- (Rmult_assoc (/ _)).
   rewrite Rinv_l, Rmult_1_l, Plus.plus_comm;[reflexivity| apply not_0_INR; lia].
  replace (sum_f_R0 (fun n => (-1) ^ n * (1 - 1) ^ (n + 1) / (INR (n + 1)))
              (2 * k + 1)) with 0 in qc.
  rewrite ln_1, !Rminus_0_r in qc; rewrite qc.
  rewrite <- (sum_eq (fun n => ((-1) * (c - 1))^n)), tech3; try psatzl R.
   replace (1 - -1 * (c - 1)) with c by ring.
   replace (S (2 * k + 1)) with (2 * (S k))%nat by ring.
   rewrite Rpow_mult_distr, pow_1_even, Rmult_1_l.
   replace (/c) with (1 * /c) by ring.
   unfold Rdiv; rewrite <- Rmult_minus_distr_r.
   assert (0 < /c) by (apply Rinv_0_lt_compat; psatzl R).
   repeat (apply Rmult_le_pos;[ | psatzl R]).
   assert (0 <= (c - 1) ^ (2 * S k)) by (apply pow_le; psatzl R); psatzl R.
  intros; rewrite Rpow_mult_distr; reflexivity.
 symmetry; apply sum_eq_R0.
 intros; replace (n + 1)%nat with (S n) by ring; simpl. 
 rewrite Rminus_eq_0; unfold Rdiv; rewrite !Rmult_0_l, Rmult_0_r, Rmult_0_l.
 reflexivity.
rewrite <- cx, sum_eq_R0, ln_1;[psatzl R | ].
intros; rewrite Plus.plus_comm, Rminus_eq_0; simpl; unfold Rdiv.
rewrite Rmult_0_l, Rmult_0_r, Rmult_0_l; reflexivity.
Qed.
  
Section bound.

Variables bup_div blo_div bup_mul blo_mul : R -> R -> R.

Hypothesis bup_div_spec : forall a b, 0 <= a -> 1 <= b -> a / b <= bup_div a b.

Hypothesis blo_div_spec : forall a b, 0 <= a -> 1 <= b -> blo_div a b <= a / b.

Hypothesis bup_mul_spec : forall a b, 0 <= a -> 0 <= b -> a * b <= bup_mul a b.

Hypothesis blo_mul_spec : forall a b, 0 <= a -> 0 <= b -> 
    0 <= blo_mul a b <= a * b.

Fixpoint bup_pow x n :=
  match n with 0 => 1 | 1 => x | S p => bup_mul x (bup_pow x p) end.

Fixpoint blo_pow x n :=
  match n with 0 => 1 | 1 => x |S p => blo_mul x (blo_pow x p) end.

Lemma bup_pow_spec x n : 0 <= x -> x ^ n <= bup_pow x n.
Proof.
induction n as [ | p].
 now simpl; intros; apply Rle_refl.
case (eq_nat_dec p 0).
 now intros p0 x0; rewrite p0; simpl; rewrite Rmult_1_r; apply Rle_refl.
intros pn0 x0; replace (bup_pow x (S p)) with (bup_mul x (bup_pow x p)).
 assert (x ^ p <= bup_pow x p) by (apply IHp; assumption).
simpl; apply Rle_trans with (x * bup_pow x p).
  apply Rmult_le_compat_l; psatzl R.
 apply bup_mul_spec;[assumption | ].
 apply Rle_trans with (x ^ p);[apply pow_le | ]; assumption.
destruct p; [case pn0|]; reflexivity.
Qed.

Lemma blo_pow_spec x n : 0 <= x -> 0 <= blo_pow x n <= x ^ n.
Proof.
induction n as [ | p].
 now simpl; intros; split;[apply Rle_0_1 | apply Rle_refl].
case (eq_nat_dec p 0).
 now intros p0 x0; rewrite p0; simpl; rewrite Rmult_1_r;
 split; [assumption | apply Rle_refl].
intros pn0 x0; replace (blo_pow x (S p)) with (blo_mul x (blo_pow x p)).
 assert (0 <= blo_pow x p) by (apply IHp; assumption).
 split.
  simpl; apply blo_mul_spec; assumption.
 apply Rle_trans with (x * blo_pow x p).
  apply blo_mul_spec; assumption.
 simpl; apply Rmult_le_compat_l;[assumption | ].
 apply IHp; assumption.
destruct p; [case pn0|]; reflexivity.
Qed.

Fixpoint bln_interval x rank :=
  match rank with
    0 => (x, x, true)
  | S p => 
    let '(b1, b2, even) := bln_interval x p in
    if even then (* S p is now odd *)
       (b1 - bup_div (bup_pow x (S (S p))) (INR (S (S p))), 
        b2 - blo_div (blo_pow x (S (S p))) (INR (S (S p))), negb even)
    else (* S p is now even *)
       (b1 + blo_div (blo_pow x (S (S p))) (INR (S (S p))),
        b2 + bup_div (bup_pow x (S (S p))) (INR (S (S p))), negb even)
  end.

Lemma bln_interval_step x rank :
  bln_interval x rank =
  match rank with
    0 => (x, x, true)
  | S p => 
    let '(b1, b2, even) := bln_interval x p in
    if even then (* S p is now odd *)
       (b1 - bup_div (bup_pow x (S (S p))) (INR (S (S p))), 
        b2 - blo_div (blo_pow x (S (S p))) (INR (S (S p))), negb even)
    else (* S p is now even *)
       (b1 + blo_div (blo_pow x (S (S p))) (INR (S (S p))),
        b2 + bup_div (bup_pow x (S (S p))) (INR (S (S p))), negb even)
  end.
Proof.
destruct rank as [| rank]; reflexivity.
Qed.

Lemma bln_interval_even x rank :
  forall v1 v2 b, bln_interval x rank = (v1, v2, b) ->
  (b = true -> exists k, rank = (2 * k)%nat) /\
  (b = false -> exists k, rank = (S (2 * k))).
Proof.
induction rank as [ | rank IHrank].
 simpl; intros v1 v2 b h; split; injection h; intros; subst;
  try discriminate; exists 0%nat; ring.
intros v1 v2 b; rewrite bln_interval_step; simpl. 
destruct (bln_interval x rank) as [[v'1 v'2] b'].
destruct (IHrank v'1 v'2 b' (refl_equal _)) as [IH1 IH2].
case_eq b'; [intros b't | intros b'f]; intros h; injection h;
 intros bq _ _; rewrite <- bq; split; try discriminate.
 apply IH1 in b't; destruct b't as [k qk]; rewrite qk;
 exists k; ring.
apply IH2 in b'f; destruct b'f as [k qk]; rewrite qk;
 exists (S k); ring.
Qed.

Lemma bln_correct x r : 
    0 <= x ->
    let '(v1, v2, b) := bln_interval x r in
    v1 <= sum_f_R0 (fun k => (-1) ^ k * x ^ (S k) / INR (S k)) r
       <= v2.
induction r; intros x0;[simpl; psatzl R | ].
assert (t := bln_interval_even x r).
rewrite bln_interval_step, tech5.
destruct (bln_interval x r) as [[v1 v2] b].
destruct (t v1 v2 b (refl_equal _)) as [t1 t2]; clear t.
assert (1 <= INR (S (S r))) by (apply (le_INR 1); lia).
assert (0 <= x ^ S (S r)) by (apply pow_le; assumption).
case_eq b; intros bq; try apply t1 in bq; try apply t2 in bq;
 destruct bq as [k qk]; 
 try replace ((-1) ^ S r) with ((-1) ^ (S (2 * k))) by
      (rewrite qk; reflexivity);
 try replace ((-1) ^ S r) with ((-1) ^ (2 * S k)) by
      (rewrite qk; apply f_equal; ring);
 rewrite ?pow_1_odd, ?pow_1_even, ?Ropp_mult_distr_l_reverse,
   ?Ropp_div, ?Rmult_1_l;
 split; (apply Rplus_le_compat; [tauto | ]); try apply Ropp_le_contravar.
   apply Rle_trans with ((bup_pow x (S (S r))) / (INR (S (S r)))).
    apply Rmult_le_compat_r;[apply Rlt_le, Rinv_0_lt_compat;psatzl R | ].
   apply bup_pow_spec; psatzl R.
   apply bup_div_spec.
    apply Rle_trans with (x ^ (S (S r)));
       [assumption | apply bup_pow_spec; psatzl R].
   assumption.
  apply Rle_trans with (blo_pow x (S (S r)) / (INR (S (S r)))).
   apply blo_div_spec;[apply blo_pow_spec; assumption | psatzl R].
  apply Rmult_le_compat_r;[apply Rlt_le, Rinv_0_lt_compat; psatzl R | ].
  apply blo_pow_spec; assumption.
 apply Rle_trans with (blo_pow x (S (S r)) / (INR (S (S r)))).
  apply blo_div_spec;[apply blo_pow_spec; assumption | psatzl R].
 apply Rmult_le_compat_r;[apply Rlt_le, Rinv_0_lt_compat; psatzl R | ].
 apply blo_pow_spec; assumption.
apply Rle_trans with (bup_pow x (S (S r)) / INR(S (S r))). 
 apply Rmult_le_compat_r;[apply Rlt_le, Rinv_0_lt_compat;psatzl R | ].
 apply bup_pow_spec; psatzl R.
apply bup_div_spec.
 apply Rle_trans with (x ^ (S (S r)));
       [assumption | apply bup_pow_spec; psatzl R].
assumption.
Qed.

End bound.

Section integers.

Variable precision : Z.

Hypothesis precision_big : (1 < precision)%Z.

Definition lo_div x y := Z.div (precision * x) y.

Definition lo_mul x y := Z.div (x * y) precision.

Definition up_div x y :=
  let (q, r) := Z.div_eucl (precision * x) y in
  if Z.eqb r 0 then q else (q + 1)%Z.

Definition up_mul x y :=
  let (q, r) := Z.div_eucl (x * y) precision in
  if Z.eqb r 0 then q else (q + 1)%Z.

Definition floor : R -> Z := Rcomplements.floor.

Lemma floorP (x : R) :  (IZR (floor x) <= x < IZR (floor x) + 1)%R.
Proof.
unfold floor, Rcomplements.floor.
destruct floor_ex as [v vp]; simpl; psatzl R.
Qed.

Definition ceiling (x : R) := (Rcomplements.floor1 x + 1)%Z.

Lemma ceilingP (x : R) :  (IZR (ceiling x) - 1 < x <= IZR (ceiling x))%R.
Proof.
unfold ceiling, Rcomplements.floor1.
destruct floor1_ex as [v vp]; rewrite plus_IZR; simpl; psatzl R.
Qed.

Definition hR (x : Z) := IZR x/IZR precision.

Definition Rh (v : R) : Z := floor ( v * IZR precision).

Definition Rh' (v : R) : Z := ceiling (v * IZR precision).

Lemma floor_pos : forall x, 0 <= x -> (0 <= floor x)%Z.
Proof.
intros x x0; destruct (floorP x) as [_ t].
assert (0 < floor x + 1)%Z;[|lia].
apply lt_IZR; rewrite plus_IZR; apply Rle_lt_trans with (2 := t).
assumption.
Qed.

Lemma floor_int : forall y, floor (IZR y) = y.
intros; apply eq_IZR.
destruct (floorP (IZR y)) as [H H1].
revert H1; change 1 with (IZR 1); rewrite <- plus_IZR.
intros H1; apply lt_IZR in H1. 
apply Rle_antisym; auto; apply IZR_le; lia.
Qed.

Lemma ceiling_int : forall y, ceiling (IZR y) = y.
intros; apply eq_IZR.
destruct (ceilingP (IZR y)) as [H H1].
revert H; change 1 with (IZR 1); rewrite <- minus_IZR.
intros H; apply lt_IZR in H. 
apply Rle_antisym; auto; apply IZR_le; lia.
Qed.

Lemma floor_mult :
   forall x y, (1 <= y)%Z -> (floor x * y <= floor (x * IZR y))%Z.
Proof.
intros x y cy.
assert (floor x * y < floor (x * IZR y) + 1)%Z;[ | lia].
 destruct (floorP (x * IZR y)) as [H H1].
apply lt_IZR; rewrite plus_IZR; simpl IZR; apply Rle_lt_trans with (2 := H1).
rewrite mult_IZR; apply Rmult_le_compat_r.
 apply (IZR_le 0); lia.
assert (t := floorP x); tauto.
Qed.

Lemma ceiling_mult :
   forall x y, (1 <= y)%Z -> (ceiling (x * IZR y) <= ceiling x * y)%Z.
Proof.
intros x y cy.
assert (ceiling (x * IZR y) - 1 < ceiling x * y)%Z;[ | lia].
 destruct (ceilingP (x * IZR y)) as [H H1].
apply lt_IZR; rewrite minus_IZR; simpl IZR; apply Rlt_le_trans with (1 := H).
rewrite mult_IZR; apply Rmult_le_compat_r.
 apply (IZR_le 0); lia.
assert (t := ceilingP x); tauto.
Qed.

Lemma up_div_spec : forall x y, 0 <= x -> 1 <= y ->
      x / y <= hR (up_div (Rh' x) (Rh y)).
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros x y x0 y1; unfold hR, up_div, Rh', Rh.
case_eq (Z.div_eucl (precision * ceiling (x * IZR precision))
            (floor (y * IZR precision))).
intros q r deq.
assert (H : (floor y * floor (IZR precision) <= floor (y * IZR precision))%Z).
 rewrite floor_int; apply floor_mult; lia.
assert (yp0 : (0 < floor (y * IZR precision))%Z).
 apply Z.lt_le_trans with (2 := H).
 apply Z.mul_pos_pos.
   elim (floorP y); change 1 with (IZR 1); rewrite <- plus_IZR; intros _ t.
  assert (tmp : IZR 1 < IZR (floor y + 1)) by (simpl IZR; psatzl R).
  apply lt_IZR in tmp; lia.
 rewrite floor_int; lia.
apply Rmult_le_reg_r with (1 := posprecision).
unfold Rdiv at 2; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [ | apply Rgt_not_eq, posprecision].
apply Rmult_le_reg_r with (IZR (floor (y * IZR precision))).
 apply (IZR_lt 0); assumption.
rewrite <- mult_IZR.
unfold Rdiv; rewrite (Rmult_assoc x), (Rmult_comm (/y)), <- (Rmult_assoc x).
rewrite (Rmult_assoc (x * _)), (Rmult_comm (x * _)).
rewrite mult_IZR.
destruct (floorP (y * IZR precision)).
apply Rle_trans with ((/y * y * IZR precision) * (x * IZR precision)).
 apply Rmult_le_compat_r.
  apply Rmult_le_pos; [ | apply Rlt_le]; assumption.
 rewrite Rmult_assoc; apply Rmult_le_compat_l;[lt0 | assumption].
rewrite Rinv_l, Rmult_1_l;[ | lt0].
rewrite <- mult_IZR.
assert (deq'' := Zdiv.Z_div_mod (precision * ceiling (x * IZR precision))
   (floor (y * IZR precision))); rewrite deq in deq''.
destruct deq'' as [deq2 deq3];[apply Z.lt_gt; assumption | ].
assert (xu : x * IZR precision <= IZR (ceiling (x * IZR precision))).
 destruct (ceilingP (x * IZR precision)); tauto.
apply Rle_trans with (IZR precision * IZR(ceiling (x * IZR precision))).
 apply Rmult_le_compat_l with (1 := Rlt_le _ _ posprecision); assumption.
rewrite <- mult_IZR, deq2, Z.mul_comm.
 apply IZR_le.
case_eq (r =? 0)%Z; intros remP.
 now rewrite Z.eqb_eq in remP; rewrite remP, Z.add_0_r; apply Z.le_refl.
rewrite Z.mul_add_distr_r, Z.mul_1_l; lia.
Qed.

Lemma lo_div_spec : forall x y, 0 <= x -> 1 <= y ->
      hR (lo_div (Rh x) (Rh' y)) <= x / y.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros x y x0 y1; unfold hR, lo_div, Rh', Rh.
unfold Z.div.
case_eq (Z.div_eucl (precision * floor (x * IZR precision))
            (ceiling (y * IZR precision))).
intros q r deq; simpl fst.
assert (H : (ceiling (y * IZR precision)
          <= ceiling y * ceiling (IZR precision))%Z).
 rewrite ceiling_int; apply ceiling_mult; lia.
assert (yp0 : (0 < ceiling (y * IZR precision))%Z).
 apply Z.lt_le_trans with (ceiling (IZR precision)).
  rewrite ceiling_int; apply lt_IZR; assumption.
 destruct (ceilingP (y * IZR precision)) as [_ t].
 apply le_IZR, Rle_trans with (2 := t); rewrite ceiling_int.
 pattern (IZR precision) at 1; rewrite <- Rmult_1_l.
 apply Rmult_le_compat_r; psatzl R.
apply Rmult_le_reg_r with (1 := posprecision).
unfold Rdiv at 1; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [ | apply Rgt_not_eq, posprecision].
apply Rmult_le_reg_r with (IZR (ceiling (y * IZR precision))).
 apply (IZR_lt 0); assumption.
rewrite <- mult_IZR.
unfold Rdiv; rewrite (Rmult_assoc x), (Rmult_comm (/y)), <- (Rmult_assoc x).
rewrite (Rmult_assoc (x * _)), (Rmult_comm (x * _)).
rewrite mult_IZR.
destruct (ceilingP (y * IZR precision)).
apply Rle_trans with ((/y * y * IZR precision) * (x * IZR precision)).
 rewrite Rinv_l, Rmult_1_l;[ | lt0].
 rewrite <- mult_IZR.
 assert (deq'' := Zdiv.Z_div_mod (precision * floor (x * IZR precision))
   (ceiling (y * IZR precision))); rewrite deq in deq''.
 destruct deq'' as [deq2 deq3];[apply Z.lt_gt; assumption | ].
 assert (xu : IZR (floor (x * IZR precision)) <= x * IZR precision).
  destruct (floorP (x * IZR precision)); tauto.
 apply Rle_trans with (IZR precision * IZR (floor (x * IZR precision))).
  rewrite <- mult_IZR, deq2.
  assert (0 <= IZR r) by (apply (IZR_le 0); tauto).
  rewrite Z.mul_comm, plus_IZR; psatzl R.
 apply Rmult_le_compat_l with (1 := Rlt_le _ _ posprecision); assumption.
apply Rmult_le_compat_r.
 apply Rmult_le_pos; [ | apply Rlt_le]; assumption.
rewrite Rmult_assoc; apply Rmult_le_compat_l;[lt0 | assumption ].
Qed.

Lemma up_mul_spec a b : 0 <= a -> 0 <= b ->
  a * b <= hR (up_mul (ceiling (a * IZR precision))
                   (ceiling (b * IZR precision))).
Proof.
intros a0 b0.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
unfold hR, up_mul.
case_eq (Z.div_eucl (ceiling (a * IZR precision) * 
                        ceiling (b * IZR precision)) precision).
intros q r deq.
apply (Rmult_le_reg_r (IZR precision));[assumption | ].
unfold Rdiv; rewrite (Rmult_assoc _ (/_)), Rinv_l, Rmult_1_r;[|psatzl R].
assert (deq'' := (Zdiv.Z_div_mod (ceiling (a * IZR precision) * 
                               ceiling (b * IZR precision)) precision) 
          (Z.lt_gt _ _ (lt_IZR 0 _ posprecision)));
  rewrite deq in deq''.
apply Rmult_le_reg_r with (IZR precision);[assumption | ].
apply Rle_trans with (IZR (ceiling (a * IZR precision) * 
                        ceiling (b * IZR precision))).
 assert (help : forall a b c : R, a * b * c * c = (a * c) * (b * c))
    by (intros; ring); rewrite help, mult_IZR.
 apply Rmult_le_compat; try (apply Rmult_le_pos; psatzl R).
  destruct (ceilingP (a * IZR precision)); assumption.
 destruct (ceilingP (b * IZR precision)); assumption.
destruct deq'' as [deq2 deq3].
rewrite deq2; clear deq deq2.
case_eq (r =? 0)%Z; intros testr.
 apply Z.eqb_eq in testr.
 rewrite testr, Zplus_0_r, mult_IZR, Rmult_comm; psatzl R.
rewrite !plus_IZR, mult_IZR, Rmult_comm, Rmult_plus_distr_r.
apply Rplus_le_compat_l.
simpl (IZR 1); rewrite Rmult_1_l; apply IZR_le; lia.
Qed.

Lemma lo_mul_spec a b : 0 <= a -> 0 <= b ->
  0 <= hR (lo_mul (floor (a * IZR precision))
                   (floor (b * IZR precision))) <= a * b.
Proof.
intros a0 b0.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
 unfold hR, lo_mul.
case_eq (Z.div_eucl (floor (a * IZR precision) * 
                        floor (b * IZR precision)) precision).
intros q r deq.
assert (deq'' := (Zdiv.Z_div_mod (floor (a * IZR precision) * 
                              floor (b * IZR precision)) precision) 
          (Z.lt_gt _ _ (lt_IZR 0 _ posprecision)));
  rewrite deq in deq''.
split.
 apply Rmult_le_pos;[ | apply Rlt_le, Rinv_0_lt_compat; assumption].
 apply (IZR_le 0).
 apply Zdiv.Z_div_pos; [apply Z.lt_gt; lia | ].
 apply Z.mul_nonneg_nonneg.
  apply floor_pos, Rmult_le_pos; psatzl R.
 apply floor_pos, Rmult_le_pos; psatzl R.
apply (Rmult_le_reg_r (IZR precision));[assumption | ].
unfold Rdiv; rewrite (Rmult_assoc _ (/_)), Rinv_l, Rmult_1_r;[|psatzl R].
apply Rmult_le_reg_r with (IZR precision);[assumption | ].
apply Rle_trans with (IZR (floor (a * IZR precision) * 
                        floor (b * IZR precision))).
 destruct deq'' as [deq2 deq3].
 unfold Zdiv.Zdiv; rewrite deq, deq2, <- mult_IZR.
 apply IZR_le; lia.
assert (help : forall a b c : R, a * b * c * c = (a * c) * (b * c))
    by (intros; ring); rewrite help, mult_IZR.
assert (help' : (forall a, 0 < a + 1 -> 0 <= a)%Z) by (clear; intros; lia).
destruct (floorP (a * IZR precision)) as [t1a t2a].
destruct (floorP (b * IZR precision)) as [t1b t2b].
apply Rmult_le_compat; try (apply Rmult_le_pos; psatzl R); auto.
 apply (IZR_le 0), floor_pos, Rmult_le_pos; psatzl R.
apply (IZR_le 0), floor_pos, Rmult_le_pos; psatzl R.
Qed.

Fixpoint lo_pow a n :=
  match n with 0%nat => precision | S p => lo_mul a (lo_pow a p) end.

Lemma RbZ_hR : forall n, (floor (hR n * IZR precision) = n).
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros n; unfold hR, Rdiv.
now rewrite Rmult_assoc, Rinv_l, Rmult_1_r, floor_int; [ | psatzl R].
Qed.

Lemma lo_mul_spec' :
 forall a b, 0 <= a -> 0 <= b ->
 exists c,
   (lo_mul (floor (a * IZR precision))
       (floor (b * IZR precision)) =
   floor (c * IZR precision)) /\ 0 <= c /\
   hR (floor (c * IZR precision)) = c.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a b a0 b0.
unfold lo_mul.
exists (IZR (floor (a * IZR precision) * floor (b * IZR precision) / precision)
           / IZR precision).
unfold Rdiv at 1.
rewrite Rmult_assoc, Rinv_l, Rmult_1_r;[|psatzl R].
rewrite floor_int; split;[reflexivity | ].
split.
 apply Rmult_le_pos;[|lt0].
  apply (IZR_le 0), Z.div_pos.
  apply Z.mul_nonneg_nonneg.
   apply floor_pos, Rmult_le_pos; psatzl R.
  apply floor_pos, Rmult_le_pos; psatzl R.
 apply Zlt_trans with (2 := precision_big); reflexivity.
unfold Rdiv at 1; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [|psatzl R].
unfold hR; rewrite floor_int; reflexivity.
Qed.

Lemma lo_pow_spec' :
  forall a n, 0 <= a ->
  exists c, (lo_pow (floor (a * IZR precision)) n =
     floor (c * IZR precision)) /\ 0 <= c /\
    hR (floor (c * IZR precision)) = c.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a n a0; induction n as [|n IH].
 simpl; exists 1; split;
  [rewrite Rmult_1_l, floor_int; reflexivity |split; [psatzl R|]].
 unfold hR, Rdiv; rewrite Rmult_1_l, floor_int.
 apply Rinv_r; psatzl R.
destruct IH as [c [cq [c0 ci]]].
simpl.
destruct (lo_mul_spec' a c a0 c0) as [d [dq d0]].
exists d; split;[|assumption].
rewrite <- dq, cq; reflexivity.
Qed.

Lemma lo_pow_spec :
  forall a n, 0 <= a ->
   0 <= hR (lo_pow (floor (a * IZR precision)) n) <= a ^ n.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a n; induction n as [ | n IH]; intros a0.
 simpl; unfold hR, Rdiv; rewrite Rinv_r; psatzl R.
split.
 destruct (lo_pow_spec' a (S n) a0) as [c [cq c0]].
 unfold hR; apply Rmult_le_pos;
  [ | apply Rlt_le, Rinv_0_lt_compat; assumption].
 rewrite cq.
 apply (IZR_le 0), floor_pos, Rmult_le_pos; psatzl R.
simpl.
destruct (lo_pow_spec' a n a0) as [c [cq [c0 ci]]].
(* Bug? why do I have to do this by hand? *)
rewrite cq.
apply Rle_trans with 
  (a * hR (floor (c * IZR precision))).
 now rewrite ci; apply lo_mul_spec.
apply Rmult_le_compat_l;[assumption | rewrite <- cq; tauto].
Qed.

Fixpoint up_pow a n :=
  match n with 0%nat => precision | S p => up_mul a (up_pow a p) end.

Lemma up_mul_spec' :
 forall a b, 0 <= a -> 0 <= b ->
 exists c,
   (up_mul (ceiling (a * IZR precision))
       (ceiling (b * IZR precision)) =
   ceiling (c * IZR precision)) /\ 
   hR (ceiling (c * IZR precision)) = c.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a b a0 b0.
unfold up_mul.
destruct (Z.div_eucl
            (ceiling (a * IZR precision) * ceiling (b * IZR precision))
            precision) as [q r].
exists (hR (if (r =? 0)%Z then q else q + 1)).
unfold hR, Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;[|psatzl R].
rewrite !ceiling_int; split; reflexivity.
Qed.

Lemma up_pow_spec' :
  forall a n, 0 <= a ->
  exists c, (up_pow (ceiling (a * IZR precision)) n =
     ceiling (c * IZR precision)) /\ a ^ n <= c /\
    hR (ceiling (c * IZR precision)) = c.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a n a0; induction n as [|n IH].
 simpl; exists 1.
 unfold hR, Rdiv; rewrite Rmult_1_l, ceiling_int, Rinv_r;[ | psatzl R].
 split; [reflexivity | split;[apply Rle_refl | reflexivity]].
destruct IH as [c [cq [cl ci]]].
simpl; destruct 
 (up_mul_spec' a (hR (ceiling (c * IZR precision))) a0) as [d [dq di]].
 rewrite ci; apply Rle_trans with (2 := cl), pow_le; assumption.
exists d; split.
 rewrite cq, <- dq, ci; reflexivity.
split.
 apply Rle_trans with (a * c);[apply Rmult_le_compat; try psatzl R | ].
  apply pow_le; assumption.
 rewrite <- di, <- dq, ci.
 apply up_mul_spec.
  assumption.
 apply Rle_trans with (2 := cl), pow_le; assumption.
assumption.
Qed.

Lemma up_pow_spec : 
  forall a n, 0 <= a ->
  a ^ n <= hR (up_pow (ceiling (a * IZR precision)) n).
Proof.
intros a n a0; destruct (up_pow_spec' a n a0) as [d [dq [dl di]]].
rewrite dq, di; apply dl.
Qed.

Fixpoint ln_i x rank :=
  match rank with
    0 => (x, x, true)
  | S p => 
    let '(b1, b2, even) := ln_i x p in
    if even then (* S p is now odd *)
       (b1 - up_div (up_pow x (S (S p))) 
                (Z.of_nat (S (S p)) * precision),
        b2 - lo_div (lo_pow x (S (S p)))
                (Z.of_nat (S (S p)) * precision),
          negb even)%Z
    else (* S p is now even *)
       (b1 + lo_div (lo_pow x (S (S p))) 
                (Z.of_nat (S (S p)) * precision),
        b2 + up_div (up_pow x (S (S p))) 
                (Z.of_nat (S (S p)) * precision),
          negb even)%Z
  end.

Lemma ln_i_step x rank :
  ln_i x rank =
  match rank with
    0 => (x, x, true)
  | S p => 
    let '(b1, b2, even) := ln_i x p in
    if even then (* S p is now odd *)
       (b1 - up_div (up_pow x (S (S p))) 
                (Z.of_nat (S (S p)) * precision),
        b2 - lo_div (lo_pow x (S (S p)))
                (Z.of_nat (S (S p)) * precision),
          negb even)%Z
    else (* S p is now even *)
       (b1 + lo_div (lo_pow x (S (S p))) 
                (Z.of_nat (S (S p)) * precision),
        b2 + up_div (up_pow x (S (S p))) 
                (Z.of_nat (S (S p)) * precision),
          negb even)%Z
  end.
Proof. case rank; reflexivity. Qed.

(* Compute the logarithm of a fixed precision number. *)

Definition ln_approx n p :=
 let v := Z.sqrt (Z.sqrt (Z.sqrt (Z.sqrt (n * precision ^ 15)))) in
 let '(_, b, _) := ln_i (v + 1 - precision) (2 * p) in
 let '(a, _, _) := ln_i (v - precision) (S (2 * p)) in
   (16 * a, 16 * b)%Z.

Lemma hR_sub : forall a b, hR (a - b) = hR a - hR b.
intros a b; unfold hR; unfold Rdiv.
rewrite minus_IZR, Rmult_minus_distr_r.
reflexivity.
Qed.

Lemma hR_add : forall a b, hR (a + b) = hR a + hR b.
intros a b; unfold hR; unfold Rdiv.
rewrite plus_IZR, Rmult_plus_distr_r.
reflexivity.
Qed.

Lemma bup_pow_i : forall a n, (0 <= a)%Z -> hR (up_pow a n) =
  bup_pow (fun x y => hR (up_mul (ceiling (x * IZR precision))
                            (ceiling (y * IZR precision)))) (hR a) n.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a n a0.
induction n.
 simpl; unfold hR; apply Rinv_r; psatzl R.
 destruct n as [ | n].
 simpl; unfold up_mul.
 case_eq (Z.div_eucl (a * precision) precision); intros q r deq.
 assert (nnzz: (precision <> 0)%Z) by lia.
 assert (t := Z.mod_mul a precision nnzz); unfold Z.modulo in t.
 rewrite deq in t.
 assert (t' := Z.div_mod (a * precision) precision nnzz).
 unfold Z.modulo, Z.div in t'; rewrite deq, t, Zplus_0_r in t'.
 rewrite Zmult_comm in t'; apply Z.mul_reg_l in t'; auto.
 unfold hR; rewrite t, Z.eqb_refl, t'; reflexivity.
change (bup_pow
     (fun x y : R =>
      hR
        (up_mul (ceiling (x * IZR precision))
           (ceiling (y * IZR precision)))) (hR a) 
     (S (S n))) with (hR (up_mul (ceiling (hR a * IZR precision))
                       (ceiling (bup_pow
     (fun x y : R =>
      hR
        (up_mul (ceiling (x * IZR precision))
           (ceiling (y * IZR precision)))) (hR a) 
     (S n) * IZR precision)))).
change (up_pow a (S (S n))) with
  (up_mul a (up_pow a (S n))).
unfold hR at 3, Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [ | psatzl R].
rewrite <- IHn; unfold hR at 3, Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [ | psatzl R].
symmetry.
rewrite !ceiling_int; reflexivity.
Qed.

Lemma blo_pow_i : forall a n, (0 <= a)%Z -> hR (lo_pow a n) =
  blo_pow (fun x y => hR (lo_mul (floor (x * IZR precision))
                            (floor (y * IZR precision)))) (hR a) n.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros a n a0.
induction n.
 simpl; unfold hR; apply Rinv_r; psatzl R.
 destruct n as [ | n].
 simpl; unfold lo_mul.
 rewrite Z.div_mul; auto; lia.
change (blo_pow
     (fun x y : R =>
      hR
        (lo_mul (floor (x * IZR precision))
           (floor (y * IZR precision)))) (hR a) 
     (S (S n))) with (hR (lo_mul (floor (hR a * IZR precision))
                       (floor (blo_pow
     (fun x y : R =>
      hR
        (lo_mul (floor (x * IZR precision))
           (floor (y * IZR precision)))) (hR a) 
     (S n) * IZR precision)))).
change (lo_pow a (S (S n))) with
  (lo_mul a (lo_pow a (S n))).
unfold hR at 3, Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [ | psatzl R].
rewrite <- IHn; unfold hR at 3, Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;
 [ | psatzl R].
now rewrite !floor_int.
Qed.

Lemma ln_i_correct :
 forall x r, (0 <= x)%Z ->
 let '(a, b, c) := ln_i x r in
   (hR a, hR b, c) =
   bln_interval 
     (fun x y => hR (up_div (ceiling (x * IZR precision)) 
                             (floor (y * IZR precision))))
     (fun x y => hR (lo_div (floor (x * IZR precision)) 
                             (ceiling (y * IZR precision))))
     (fun x y => hR (up_mul (ceiling (x * IZR precision))
                             (ceiling (y * IZR precision))))
     (fun x y => hR (lo_mul (floor (x * IZR precision))
                             (floor (y * IZR precision))))
     (hR x) r.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros x r x0; induction r.
 simpl; reflexivity.
rewrite ln_i_step, bln_interval_step.
destruct (ln_i x r) as [[A B] C].
destruct (bln_interval 
     (fun x y => hR (up_div (ceiling (x * IZR precision)) 
                             (floor (y * IZR precision))))
     (fun x y => hR (lo_div (floor (x * IZR precision)) 
                             (ceiling (y * IZR precision))))
     (fun x y => hR (up_mul (ceiling (x * IZR precision))
                             (ceiling (y * IZR precision))))
     (fun x y => hR (lo_mul (floor (x * IZR precision))
                             (floor (y * IZR precision))))
     (hR x) r) as [[D E] F].
injection IHr; intros CF EB AD.
assert (help : forall x y z t : R, x = z -> y = t -> (x, y) = (z, t)).
 clear; intros; subst; reflexivity.
rewrite <- CF, <- EB, <- AD; case C;
[apply f_equal with (f := fun x => (x, negb true))|
 apply f_equal with (f := fun x => (x, negb false))]; apply help.
   rewrite hR_sub; apply f_equal.
   rewrite INR_IZR_INZ, <- mult_IZR, floor_int.
   apply (f_equal (fun x => hR (up_div x (Z.of_nat (S (S r)) * precision)))).
   rewrite <- bup_pow_i;[ | assumption]; unfold hR, Rdiv. 
   rewrite Rmult_assoc, Rinv_l, Rmult_1_r, ceiling_int;
    [reflexivity | psatzl R].
  rewrite hR_sub; apply f_equal.
  rewrite INR_IZR_INZ, <- mult_IZR, ceiling_int.
  apply (f_equal (fun x => hR (lo_div x (Z.of_nat (S (S r)) * precision)))).
  rewrite <- blo_pow_i;[ | assumption]; unfold hR, Rdiv. 
  rewrite Rmult_assoc, Rinv_l, Rmult_1_r, floor_int;[reflexivity | psatzl R].
 rewrite hR_add; apply f_equal.
 rewrite INR_IZR_INZ, <- mult_IZR, ceiling_int.
 rewrite <- blo_pow_i;[ | assumption]; unfold hR, Rdiv. 
 rewrite Rmult_assoc, Rinv_l, Rmult_1_r, floor_int;[reflexivity | psatzl R].
rewrite hR_add; apply f_equal.
rewrite INR_IZR_INZ, <- mult_IZR, floor_int.
rewrite <- bup_pow_i;[ | assumption]; unfold hR, Rdiv. 
rewrite Rmult_assoc, Rinv_l, Rmult_1_r, ceiling_int;[reflexivity | psatzl R].
Qed.

Lemma sqrt_ap :
  forall x, (0 <= x)%Z ->
    hR (Z.sqrt (x * precision)) <= sqrt (hR x)
    < hR (Z.sqrt (x * precision)) + /IZR precision.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros x x0.
assert (pp : 0 <= IZR (x * precision)).
 apply (IZR_le 0).
 apply Z.mul_nonneg_nonneg;[assumption | ].
 apply le_IZR; change (IZR 0) with 0; psatzl R.
generalize (Z.sqrt_spec _ (le_IZR 0 _ pp)).
unfold hR.
set (s := Z.sqrt (x * precision)); intros [sl su].
split.
 apply Rmult_le_reg_r with (IZR precision);[assumption | ].
 unfold hR, Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r;[ | psatzl R].
 apply Rsqr_incr_0; unfold Rsqr. 
   rewrite <- mult_IZR.
   rewrite Rmult_assoc, !(Rmult_comm (sqrt (_))), <- Rmult_assoc,
      <- (mult_IZR precision), Rmult_assoc, sqrt_sqrt, Rmult_comm.
    rewrite Rmult_assoc, (mult_IZR precision), <- (Rmult_assoc (/ _)),
        Rinv_l, Rmult_1_l; [ | psatzl R].
    now rewrite <- mult_IZR; apply IZR_le, sl.
   apply Rmult_le_pos;[apply (IZR_le 0); assumption | ].
   apply Rlt_le, Rinv_0_lt_compat; assumption.
  apply (IZR_le 0), Z.sqrt_nonneg.
 apply Rmult_le_pos; [apply sqrt_pos |psatzl R].
apply Rmult_lt_reg_r with (IZR precision);[assumption | ].
 unfold Rdiv; rewrite Rmult_plus_distr_r,Rmult_assoc, Rinv_l, Rmult_1_r;
  [ | psatzl R].
assert (help : forall a b, 0 <= a -> 0 < b -> sqrt (a * / b) * b = sqrt (a * b)).
 intros a b a0 b0.
 pattern b at 2; rewrite <- (sqrt_sqrt b);[|psatzl R].
 rewrite <- !sqrt_mult; try lt0.
  apply f_equal; field; try psatzl R.
 apply Rmult_le_pos; lt0.
rewrite help;[ | apply (IZR_le 0); assumption| assumption].
change 1 with (IZR 1); rewrite <- plus_IZR.
assert (0 <= s)%Z by apply Z.sqrt_nonneg.
rewrite <- mult_IZR.
apply Rsqr_incrst_0;[ | apply sqrt_pos | apply (IZR_le 0); lia].
rewrite Rsqr_sqrt; auto.
unfold Rsqr; rewrite <- mult_IZR; apply IZR_lt; assumption.
Qed.

Lemma Zsqrt_pow :
  forall x z n t, (0 <= x)%Z -> 
  IZR z ^ (NPeano.pow 2 n) <= t <
      (IZR z + 1) ^ (NPeano.pow 2 n) ->
  IZR x * IZR x <= IZR z < (IZR x + 1) * (IZR x + 1) ->
  IZR x ^ (NPeano.pow 2 (n + 1)) <= t < (IZR x + 1) ^ (NPeano.pow 2 (n + 1)).
Proof.
intros x z n t x0 cz cx.
assert (0 <= IZR x) by (apply (IZR_le 0); assumption).
assert (tx : forall u, u ^ NPeano.pow 2 (n + 1) =
             (u * u) ^ NPeano.pow 2 n).
  intros u; rewrite NPeano.Nat.pow_add_r, mult_comm.
 simpl (NPeano.pow 2 1); rewrite pow_mult.
 apply (f_equal (fun x => x ^ (NPeano.pow 2 n))); field.
rewrite !tx; split.
 apply Rle_trans with (IZR z ^ NPeano.pow 2 n);[|tauto].
 apply pow_incr; split;[apply Rmult_le_pos |]; psatzl R.
apply Rlt_le_trans with ((IZR z + 1) ^ NPeano.pow 2 n);[tauto|].
apply pow_incr.
split.
 assert (0 <= IZR z);[|psatzl R].
 apply Rle_trans with (IZR x * IZR x);[|tauto].
 apply Rmult_le_pos;psatzl R.
assert (z < (x + 1) * (x + 1))%Z.
apply lt_IZR; rewrite mult_IZR, plus_IZR; tauto.
change 1 with (IZR 1); rewrite <- !plus_IZR, <- !mult_IZR;
apply IZR_le; lia.
Qed.

Lemma reduc_correct :
  forall x, (0 <= x)%Z -> 
   let s := Z.sqrt (Z.sqrt (Z.sqrt (Z.sqrt (x * precision ^ 15)))) in
   hR s ^ 16 <= hR x < (hR s + /IZR precision) ^ 16.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros x x0 s.
replace (hR s + / IZR precision) with (hR (s + 1)) by
 (unfold hR; rewrite plus_IZR; simpl (IZR 1); field; psatzl R).
unfold hR, Rdiv.
rewrite !Rpow_mult_distr, <- Rinv_pow;[ | psatzl R].
assert (0 < IZR precision ^ 16) by (apply pow_lt; assumption).
assert (0 < IZR precision ^ 8) by (apply pow_lt; assumption).
assert (0 < IZR precision ^ 4) by (apply pow_lt; assumption).
assert (0 < IZR precision ^ 2) by (apply pow_lt; assumption).
assert (IZR s ^ 16 <= IZR x * IZR precision ^ 15 <
        (IZR s + 1) ^ 16).
 apply (Zsqrt_pow _ (Z.sqrt (Z.sqrt (Z.sqrt (x * precision ^ 15)))) 3).
   apply Z.sqrt_nonneg.
  apply (Zsqrt_pow _ (Z.sqrt (Z.sqrt (x * precision ^ 15))) 2).
    apply Z.sqrt_nonneg.
   apply (Zsqrt_pow _ (Z.sqrt (x * precision ^ 15)) 1).
     apply Z.sqrt_nonneg.
    apply (Zsqrt_pow _ (x * precision ^ 15) 0).
      apply Z.sqrt_nonneg.
     rewrite NPeano.pow_0_r, !pow_1, pow_IZR, <- mult_IZR.
     simpl Z.of_nat; psatzl R.
    rewrite <- (plus_IZR _ 1), <- !mult_IZR.
    split;[apply IZR_le | apply IZR_lt].
     apply Z.sqrt_spec.
     apply Z.mul_nonneg_nonneg; [ | apply Z.pow_nonneg]; lia.
    apply Z.sqrt_spec.
    apply Z.mul_nonneg_nonneg; [ | apply Z.pow_nonneg]; lia.

   rewrite <- (plus_IZR _ 1), <- !mult_IZR.
   split;[apply IZR_le | apply IZR_lt].
    apply Z.sqrt_spec.
    apply Z.sqrt_nonneg.
   apply Z.sqrt_spec.
   apply Z.sqrt_nonneg.

  rewrite <- (plus_IZR _ 1), <- !mult_IZR.
  split;[apply IZR_le | apply IZR_lt].
   apply Z.sqrt_spec.
   apply Z.sqrt_nonneg.
  apply Z.sqrt_spec.
  apply Z.sqrt_nonneg.

 rewrite <- (plus_IZR _ 1), <- !mult_IZR.
 split;[apply IZR_le | apply IZR_lt].
  apply Z.sqrt_spec.
  apply Z.sqrt_nonneg.
 apply Z.sqrt_spec.
 apply Z.sqrt_nonneg.
split.
 apply Rmult_le_reg_r with (IZR precision ^ 16).
  apply pow_lt; assumption.
 rewrite Rmult_assoc, Rinv_l, Rmult_1_r;[|psatzl R].
 pattern 16%nat at 2; replace 16%nat with (1 + 15)%nat by ring.
 rewrite pow_add, pow_1, Rmult_assoc, <- (Rmult_assoc (/IZR precision)).
 rewrite Rinv_l, Rmult_1_l;[tauto | psatzl R].
apply Rmult_lt_reg_r with (IZR precision ^ 16).
 apply pow_lt; assumption.
rewrite !Rmult_assoc, Rinv_l, Rmult_1_r;[|psatzl R].
pattern 16%nat at 1; replace 16%nat with (1 + 15)%nat by ring.
rewrite pow_add, pow_1, <- (Rmult_assoc (/IZR precision)).
rewrite Rinv_l, Rmult_1_l;[ | psatzl R].
rewrite plus_IZR; tauto.
Qed.

Lemma Zsqrt_mul_sqr :
  forall x y, (0 <= x)%Z -> (0 < y)%Z ->
    IZR (Z.sqrt (x * y ^ 2)) <= sqrt (IZR x) * IZR y <
    IZR (Z.sqrt (x * y ^ 2) + 1).
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
intros x y x0 y0.
replace (y ^ 2)%Z with (y * y)%Z by ring.
assert (y20 : (0 <= y * y)%Z).
 apply Z.mul_nonneg_nonneg; lia.
assert (pp : 0 <= IZR (x * (y * y))).
 apply (IZR_le 0).
 apply Z.mul_nonneg_nonneg;[assumption | lia].
generalize (Z.sqrt_spec _ (le_IZR 0 _ pp)).
set (s := Z.sqrt (x * (y * y))); intros [sl su].
split.
 apply Rsqr_incr_0; unfold Rsqr. 
   rewrite <- mult_IZR.
   rewrite Rmult_assoc, !(Rmult_comm (sqrt (_))), <- Rmult_assoc,
      <- (mult_IZR y), Rmult_assoc, sqrt_sqrt, Rmult_comm.
    now rewrite <- mult_IZR; apply IZR_le, sl.
   apply (IZR_le 0); assumption.
  apply (IZR_le 0), Z.sqrt_nonneg.
 apply Rmult_le_pos; [apply sqrt_pos | apply (IZR_le 0); lia].
assert (0 <= s)%Z by (apply Z.sqrt_nonneg).
apply Rsqr_incrst_0;[ | repeat apply Rmult_le_pos; try (apply sqrt_pos); apply (IZR_le 0); lia
    | apply (IZR_le 0); lia].
rewrite Rsqr_mult, Rsqr_sqrt;[|apply (IZR_le 0); assumption].
unfold Rsqr; rewrite <- !mult_IZR; apply IZR_lt; assumption.
Qed.
(*
Lemma Zsqrt_mul_sqr :
  (forall x y, 0 <= x -> 0 < y -> 
    Z.sqrt x * y <= Z.sqrt (x * y ^ 2) < (Z.sqrt x + 1) * y)%Z.
Proof.
intros x y x0 y0.
assert (y0' : (0 <= y)%Z) by (apply Z.lt_le_incl; auto).
assert (xy0 : (0 <= x * y ^ 2)%Z).
 apply Z.mul_nonneg_nonneg;[|apply Z.pow_nonneg]; assumption.
destruct (Z.sqrt_spec _ xy0) as [t1 t2].
destruct (Z.sqrt_spec _ x0) as [t3 t4].
assert (ys0 : (0 <= y * y)%Z).
 apply Z.mul_nonneg_nonneg; assumption.
split.
 rewrite <- Z.lt_succ_r; unfold Z.succ.
 apply Z.square_lt_simpl_nonneg. 
  apply Z.add_nonneg_nonneg;[apply Z.sqrt_nonneg | compute; discriminate].
 replace (Z.sqrt x * y * (Z.sqrt x * y))%Z with
       ((Z.sqrt x * Z.sqrt x) * y ^ 2)%Z by ring.
 apply Z.le_lt_trans with (2 := t2).
apply Z.mul_le_mono_nonneg_r;[apply Z.pow_nonneg | ]; assumption.
apply Z.square_lt_simpl_nonneg.
 apply Z.mul_nonneg_nonneg;[ | assumption].
 apply Z.add_nonneg_nonneg;[apply Z.sqrt_nonneg | ].
 compute; discriminate.
apply Z.le_lt_trans with (1 := t1).
replace ((Z.sqrt x + 1) * y * ((Z.sqrt x + 1) * y))%Z with
        ((Z.sqrt x + 1) * (Z.sqrt x + 1) * y ^ 2)%Z by ring.
apply Z.mul_lt_mono_pos_r;[ | assumption].
apply Z.pow_pos_nonneg;[assumption | compute;discriminate].
Qed.
*)

Lemma Rle_pow_R1_cancel :
  forall x n, 0 <= x -> (0 < n)%nat -> 1 <= x ^ n -> 1 <= x.
Proof.
intros x n x0 n0 cx; destruct (Rle_dec 1 x) as [x1 | x1];[psatzl R | ].
assert (x1' : x < 1) by psatzl R.
destruct (pow_lt_1_compat x n (conj x0 x1') n0); psatzl R.
Qed.

Lemma Rlt_pow_R1_cancel :
  forall x n, 0 <= x -> 1 < x ^ n -> 1 < x.
Proof.
intros x [|n] x0 cx.
 simpl in cx; psatzl R.
destruct (Rlt_dec 1 x) as [x1 | x1];[assumption | ].
assert (x <= 1) by psatzl R.
assert (st : x ^ S n <= 1 ^ S n).
 apply pow_incr; psatzl R.
case (Rle_not_lt _ _ st); rewrite pow1; assumption.
Qed.

Lemma pow_incr_cancel :
  forall x y n, 0 < x -> 0 < y -> (0 < n)%nat
       -> x ^ n <= y ^ n -> x <= y.
Proof.
intros x y n x0 y0 n0 xy.
destruct (Rle_lt_dec x y);[assumption | ].
assert (y ^ n < x ^ n);[|psatzl R].
clear xy; induction n as [|k Ik];[lia | ].
destruct k.
 simpl; psatzl R.
replace (y ^ S (S k)) with (y * y ^ (S k)) by (simpl; ring).
replace (x ^ S (S k)) with (x * x ^ (S k)) by (simpl; ring).
apply Rmult_gt_0_lt_compat.
   apply pow_lt; assumption.
  assumption.
 assumption.
apply Ik; lia.
Qed.

Lemma ln_approx_correct :
  forall x, (precision < x)%Z ->
  forall p, let (a, b) := ln_approx x p in
    hR a <= ln (hR x) <= hR b.
Proof.
assert (posprecision : 0 < IZR precision).
 apply (IZR_lt 0); lia.
assert (posVpr : 0 < / IZR precision).
 apply Rinv_0_lt_compat; psatzl R.
intros x px64 p; unfold ln_approx.
assert (hrx1 : 1 < hR x).
 unfold hR, Rdiv; apply (Rmult_lt_reg_r (IZR precision)); auto.
 rewrite Rmult_assoc, Rinv_l, Rmult_1_r, Rmult_1_l; [|psatzl R].
 apply IZR_lt; assumption.
assert (x0 : (0 <= x)%Z) by lia.
generalize (reduc_correct _ x0).
set (s := Z.sqrt (Z.sqrt (Z.sqrt (Z.sqrt (x * precision ^ 15))))).
intros rc; lazy zeta in rc.
assert (smp0 : (0 <= s - precision)%Z).
 assert (precision < s + 1)%Z;[ | lia].
 apply lt_IZR; apply Rmult_lt_reg_r with (/ IZR precision).
  assumption.
 rewrite Rinv_r;[|psatzl R].
 rewrite plus_IZR, Rmult_plus_distr_r, Rmult_1_l.
 change (IZR s * /IZR precision) with (hR s).
 apply Rlt_pow_R1_cancel with 16%nat.
  apply Rplus_le_le_0_compat;[|psatzl R].
  apply Rmult_le_pos; [apply (IZR_le 0), Z.sqrt_nonneg| psatzl R].
 apply Rlt_trans with (hR x);[ | tauto].
 unfold hR; apply Rmult_lt_reg_r with (IZR precision);[assumption | ].
 unfold Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_l, Rmult_1_r;
   [apply IZR_lt;assumption | psatzl R].
assert (smp1 : (0 <= s + 1 - precision)%Z) by lia.
assert (smp0' : 0 <= hR (s - precision)).
 apply Rmult_le_pos; apply (IZR_le 0) || apply Rlt_le; assumption.
assert (smp1' : 0 <= hR (s + 1 - precision)).
 apply Rmult_le_pos; apply (IZR_le 0) || apply Rlt_le; assumption.
assert (ilnc1 := ln_i_correct 
                (s + 1 - precision) (2 * p) smp1).
case_eq (ln_i (s + 1 - precision) (2 * p));
 intros [A b d ilneq1].
rewrite ilneq1 in ilnc1.
assert (bln1 :=
  bln_correct _ _ _ _ up_div_spec lo_div_spec up_mul_spec lo_mul_spec
  _  (2 * p) smp1').
unfold Rh', Rh in bln1; rewrite <- ilnc1 in bln1; clear ilnc1 ilneq1.
assert (ilnc2 := ln_i_correct (s - precision) (S (2 * p)) smp0).
case_eq (ln_i (s - precision)  (S (2 * p)));
 intros [a B d' ilneq2].
rewrite ilneq2 in ilnc2.
assert (bln2 :=
  bln_correct _ _ _ _ up_div_spec lo_div_spec up_mul_spec lo_mul_spec
  _ (S (2 * p)) smp0').
unfold Rh', Rh in bln2; rewrite <- ilnc2 in bln2; clear ilnc2 ilneq2.
assert (sxp16 : hR x = sqrt (sqrt (sqrt (sqrt (hR x)))) ^ 16).
 assert (help : forall u, 0 <= u -> sqrt u ^ 2 = u).
  now clear; intros u u0; simpl; rewrite Rmult_1_r, sqrt_sqrt.
 replace 16%nat with (2 * 2 * 2 * 2)%nat by ring.
 rewrite !pow_mult, !help; try apply sqrt_pos.
  reflexivity.
 apply Rmult_le_pos;[apply (IZR_le 0); lia | psatzl R].
assert (smpa : hR (s - precision) <= sqrt (sqrt (sqrt (sqrt (hR x)))) - 1<=
        hR (s + 1 - precision)).
 assert (d1 : hR (s - precision) = hR s - 1).
  unfold hR, Rdiv; rewrite minus_IZR, Rmult_minus_distr_r, Rinv_r;psatzl R.
 assert (d2 : hR (s + 1 - precision) = hR (s + 1) - 1).
  unfold hR, Rdiv; rewrite minus_IZR, Rmult_minus_distr_r, Rinv_r;psatzl R.
 rewrite d1, d2.
 split; apply Rplus_le_compat_r.
  apply pow_incr_cancel with 16%nat.
     psatzl R.
    repeat apply sqrt_lt_R0; psatzl R.
   lia. 
  rewrite <- sxp16; tauto.
 apply pow_incr_cancel with 16%nat.
    repeat apply sqrt_lt_R0; psatzl R.
   psatzl R.
  lia. 
 rewrite <- sxp16.
 assert (app : hR (s + 1) = hR s + /IZR precision).
  unfold hR; rewrite plus_IZR; simpl (IZR 1); field; psatzl R.
 rewrite app; psatzl R.
assert (s1 : 1 <= hR s).
 unfold hR; apply Rmult_le_reg_r with (IZR precision); try psatzl R.
 unfold Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r, Rmult_1_l; try psatzl R.
 apply IZR_le; lia.
assert (s1' : 1 <= hR (s + 1)).
 unfold hR; apply Rmult_le_reg_r with (IZR precision); try psatzl R.
 unfold Rdiv; rewrite Rmult_assoc, Rinv_l, Rmult_1_r, Rmult_1_l; try psatzl R.
 apply IZR_le; lia.
assert (ln3 := ln_lb p (hR s) s1).
assert (ln4 := ln_ub p (hR (s + 1)) s1').
assert (s2p1 : (2 * p + 1 = S (2 * p))%nat) by ring.
rewrite s2p1 in ln3.
assert (h16 : forall u, hR (16 * u) = 16 * hR u).
 intros u; unfold hR, Rdiv; rewrite mult_IZR; simpl (IZR 16); ring.
rewrite !h16, sxp16, ln_pow; [ | repeat apply sqrt_lt_R0; psatzl R].
replace (INR 16) with 16 by (simpl; ring).
split; (apply Rmult_le_compat_l;[psatzl R|]).
 apply Rle_trans with (1 := proj1 bln2).
 apply Rle_trans with (ln (hR s)).
  apply Rle_trans with (2 := ln3).
  apply Req_le, sum_eq; intros i _; replace (i + 1)%nat with (S i) by ring.
  apply f_equal with (f := fun a => a / INR (S i)); apply f_equal.
  apply f_equal with (f := fun a => a ^ S i).
  unfold hR; rewrite minus_IZR; field; psatzl R.
 assert (sx : hR s <= sqrt (sqrt (sqrt (sqrt (hR x))))).
  apply (fun x y => pow_incr_cancel x y 16); try psatzl R.
  lia.
 destruct (Rle_lt_or_eq_dec _ _ sx) as [sx' | sx'].
  apply Rlt_le, ln_increasing; psatzl R.
 apply Req_le; rewrite sx'; reflexivity.
apply Rle_trans with (2 := proj2 bln1).
 apply Rle_trans with (ln (hR (s + 1))).
 assert (sx : sqrt (sqrt (sqrt (sqrt (hR x)))) <= hR (s + 1)).
  apply (fun x y => pow_incr_cancel x y 16); try psatzl R.
  lia.
  replace (hR (s + 1)) with (hR s + /IZR precision); try psatzl R.
  unfold hR; rewrite plus_IZR; simpl IZR; field; psatzl R.
 destruct (Rle_lt_or_eq_dec _ _ sx) as [sx' | sx'].
  apply Rlt_le, ln_increasing; psatzl R.
 apply Req_le; rewrite sx'; reflexivity.
apply Rle_trans with (1 := ln4).
 apply Req_le, sum_eq; intros i _; replace (i + 1)%nat with (S i) by ring.
 apply f_equal with (f := fun a => a / INR (S i)); apply f_equal.
 apply f_equal with (f := fun a => a ^ S i).
 unfold hR; rewrite minus_IZR; field; psatzl R.
Qed.

End integers.