text:
(*
(C) Copyright 2010, COQTAIL team

Project Info: http://sourceforge.net/projects/coqtail/

This library is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301,
USA.
*)

Require Import Max.
Require Import Reals.
Require Import Rintegral.
Require Import Rintegral_usual.
Require Import Rintegral_tactic.
Require Import Fourier.
Require Import Rsequence_facts.
Require Import Rsequence_subsequence.
Require Import Rpser_facts.
Require Import RTaylor.
Require Import MyRfunctions.
Require Import Rsequence_tactics.

Open Local Scope R_scope.

Section Wallis.

Local Coercion INR : nat >-> R.

Definition sin_n n := fun x => pow (sin x) n.

Lemma integrable_sin_n n : Riemann_integrable (comp (fun x => pow x n) sin) 0 (PI/2).
Proof.
intro n.
apply RiemannInt_P6.
pose proof PI_RGT_0; fourier.
intros x Hx; apply continuity_pt_comp.
  apply continuity_sin.
  apply derivable_continuous; apply derivable_pow.
Qed.

Definition W_even n := ((PI/2)*(fact (2 * n))/(2 ^ (2 * n) * (fact n) ^2 )).
Definition W_odd n := ((2^(2*n) * (fact n) ^2)/(fact (S (2 * n)))).

Lemma Wallis_0 : Rint (sin_n 0)  0 (PI/2) (PI/2).
Proof.
apply Rint_eq_compat with (fct_cte 1).
  trivial.
assert(Hrew : (PI/2) = (1*(PI/2 -0))) by ring; rewrite Hrew at 2; clear Hrew.
auto with Rint.
Qed.

Lemma Wallis_1 : Rint (sin_n 1)  0 (PI/2) 1.
Proof.
unfold sin_n.
apply Rint_eq_compat with sin.
  auto with *.
assert(Heq : 1 = (cos 0 - cos (PI/2))).
  rewrite cos_PI2, cos_0; ring.
rewrite Heq at 3.
auto with Rint.
Qed.

(** Recurrence formula *)
Lemma Wallis_formula Wn n:
  Rint (sin_n n) 0 (PI/2) Wn->
    Rint (sin_n (S (S n))) 0 (PI/2) (Wn * (S n)/(S (S n))).
Proof.
intros Wn n H.
pose (RiemannInt (integrable_sin_n (S (S n)))) as W2n.
pose proof (Rint_RiemannInt_link _ _ _ (integrable_sin_n (S (S n)))) as HW2.
assert (W2n = (S n) * (Wn - W2n)) as Heq.
replace (S n * (Wn - W2n)) with
  (comp (fun x => x ^ (S n)) sin (PI/2) *((- cos)%F (PI/2)) - comp (fun x => x ^ (S n)) sin 0 *((- cos)%F 0) - (-(S n) * (Wn - W2n))).
eapply Rint_parts with (f' := fun x => ((S n) * sin x ^ n)* cos x) (g' := sin).
  left; apply PI2_RGT_0.

  intros.
  apply derivable_pt_lim_comp.
    auto with Rcont.
    apply derivable_pt_lim_pow.

  intros.
  replace (sin x) with (- - sin x) by ring.
  auto with Rcont.

  intros.
  apply continuity_pt_mult.
    apply continuity_pt_mult; auto with Rcont.
      apply continuity_pt_const; intros u v; reflexivity.
      reg.
    auto with Rcont.
    auto with Rcont.

  apply Rint_eq_compat with (comp (fun x : R => x ^ S (S n)) sin).
  intros; unfold comp, mult_fct; simpl; ring.
  apply Rint_RiemannInt_link.

  apply Rint_eq_compat with (fun x => (- S n) * (sin x ^ n - sin x ^ (S (S n)))).
  intros.
  unfold mult_fct, opp_fct.
  rewrite Rmult_assoc.
  replace (cos x * - cos x) with (-(1 - (sin x)^2)).
  unfold Rsqr; simpl; ring.
  replace ((sin x) ^ 2) with ((sin x)²) by (unfold Rsqr; ring).
  rewrite <- (cos2 x); unfold Rsqr; ring.
  apply Rint_scalar_mult_compat_l.
  apply Rint_minus; trivial.

  unfold opp_fct, comp;
  rewrite cos_0, cos_PI2, sin_PI2, sin_0;
  simpl; ring.

replace (Wn * S n / S (S n)) with (W2n).
  trivial.
rewrite Rmult_comm.
replace W2n with (/(S (S n)) * (S (S n))* W2n).
replace Wn with ((Wn - W2n) + W2n) by ring.
rewrite Rmult_plus_distr_l.
rewrite <- Heq.
unfold Rdiv.
rewrite Rmult_comm with (r1 := (W2n + S n * W2n)).
rewrite Rmult_assoc.
apply Rmult_eq_compat_l.
simpl; ring.
field.
auto with *.
Qed.

Lemma Wallis_even n : Rint (sin_n (2 * n)) 0 (PI/2) (W_even n).
Proof.
unfold W_even.
induction n.
  simpl; unfold Rsqr; field_simplify; apply Wallis_0.

apply Rint_eq_compat with (sin_n (S (S (2 *  n)))).
  intros; simpl;
  replace (S (n + S (n + 0))) with (S (S (n + (n + 0)))) by auto;
  reflexivity.

replace (PI / 2 * fact (2 * S n) / (2 ^ (2 * S n) * (fact (S n)) ^2)) with
  ((PI / 2 * fact (2 * n) / (2 ^ (2 * n) * (fact n)^2)) * S (2 * n) / (S (S (2 * n)))).
apply Wallis_formula; apply IHn.

replace (2 * S n)%nat with (S (S (2 * n))) by auto with *.
repeat rewrite fact_simpl.
repeat rewrite <- tech_pow_Rmult.
unfold Rsqr.
replace (S (S (2*n))) with (2 * (S n))%nat by auto with *.
repeat rewrite mult_INR.
replace (INR 2) with 2 by trivial; field.

repeat split.
apply not_0_INR; apply fact_neq_0.
auto with *.
apply pow_nonzero; pose proof Rlt_R0_R2; intro; fourier.
Qed.

Lemma Wallis_odd n : Rint (sin_n (S (2 * n))) 0 (PI/2) ((2^(2*n) * (fact n) ^2)/(fact (S (2 * n)))).
Proof.
unfold W_odd.
induction n.
  simpl; replace (1 * (1 * (1*1)) / 1) with 1 by field; apply Wallis_1.
apply Rint_eq_compat with (sin_n (S (S (S (2 * n))))).
  intros; simpl.
replace (n + S (n + 0))%nat with (S (n + (n + 0))) by auto;
reflexivity.
replace (2 ^ (2 * S n) * (fact (S n))^2 / fact (S (2 * S n))) with
  ((2 ^ (2 * n) * (fact n)^2 / fact (S (2 * n))) * (S (S (2 * n)))/(S (S (S (2 * n))))).
apply Wallis_formula; apply IHn.

replace (2 * S n)%nat with (S (S (2 * n))) by auto with *.
repeat rewrite fact_simpl.
repeat rewrite <- tech_pow_Rmult.
unfold Rsqr.
replace (S (S (2*n))) with (2 * (S n))%nat by auto with *.
repeat rewrite mult_INR.
replace (INR 2) with 2 by trivial; field.

repeat split;apply not_0_INR; try (auto with *).
apply fact_neq_0.
Qed.

Lemma Wallis_odd_le_even n : W_odd n <= W_even n.
Proof.
intro n.
apply (Rint_le_compat (sin_n (S (2 * n))) (sin_n (2 * n)) 0 (PI/2)).
  left; apply PI2_RGT_0.

  intros; unfold sin_n.
  replace (sin u ^ (2 * n)) with (1 * sin u ^ (2 * n)) by ring.
  rewrite <- tech_pow_Rmult.
  apply Rmult_le_compat_r.
    apply pow_le.
    apply sin_ge_0; intuition; fourier.
    pose proof (SIN_bound u); intuition.

apply Wallis_odd.
apply Wallis_even.
Qed.

Lemma Wallis_even_le_odd n : W_even (S n) <= W_odd n.
Proof.
intro n.
apply (Rint_le_compat (sin_n (2 * (S n))) (sin_n (S (2 * n))) 0 (PI/2)).
  left; apply PI2_RGT_0.

  intros; unfold sin_n.
  replace (2 * (S n))%nat with (S (S (2 * n))) by auto with *.
  replace (sin u ^ (S (2 * n))) with (1 * sin u ^ (S (2 * n))) by ring.
  rewrite <- tech_pow_Rmult.
  apply Rmult_le_compat_r.
    apply pow_le.
    apply sin_ge_0; intuition; fourier.
    pose proof (SIN_bound u); intuition.

apply Wallis_even.
apply Wallis_odd.
Qed.

Lemma Wallis_maj (n : nat) : (2*n/ (S (2 * n))) * (W_even n) <= W_odd n.
Proof.
intro n.
destruct n.

  unfold W_even, W_odd, Rsqr.
  simpl; field_simplify; fourier.

   replace 2 with (INR 2) by trivial; rewrite <- mult_INR.
  replace (2 * S n)%nat with (S (S (2 * n))) by auto with *.

  replace (W_odd (S n)) with (INR (S (S (2 * n))) / INR (S (S (S (2 * n)))) * W_odd n).

  apply Rmult_le_compat_l.
  unfold Rdiv.
  repeat apply Rmult_le_pos; auto with *.

  apply Wallis_even_le_odd.

  rewrite Rmult_comm.
  unfold Rdiv.
  rewrite <- Rmult_assoc.
  eapply (Rint_uniqueness _ 0 (PI/2)).
  apply Wallis_formula.
  apply Wallis_odd.
  replace (S (S (2 * n))) with (2 * (S n))%nat by auto with *.
  apply Wallis_odd.
Qed.

Lemma W_even_pos : forall n, 0 < W_even n.
Proof.
intro n.
 unfold W_even, Rsqr, Rdiv.
 repeat apply Rmult_lt_0_compat.
 fourier.
 fourier.
 apply Rmult_lt_reg_l with 4.
 fourier.
 rewrite Rmult_0_r ; apply PI_RGT_0.
 fourier.
 apply INR_fact_lt_0.
 repeat rewrite Rinv_mult_distr ; repeat apply Rmult_lt_0_compat.
 rewrite Rinv_pow.
 apply pow_lt ; fourier.
 apply Rgt_not_eq ; fourier.
 apply Rinv_0_lt_compat, pow_lt, INR_fact_lt_0.
 apply Rgt_not_eq, pow_lt ; fourier.
apply Rgt_not_eq, pow_lt, INR_fact_lt_0.
Qed.

Lemma Wallis_bound : forall (n : nat),  2*n / (S(2 * n)) <= W_odd n / W_even n <= 1.
Proof.
intro n ; split.
 apply Rle_trans with (2 * n / S (2 * n) * (W_even n / W_even n)).
 right ; field ; split.
 apply Rgt_not_eq ; apply W_even_pos.
 apply not_0_INR ; discriminate.
 unfold Rdiv ; rewrite <- Rmult_assoc ; apply Rmult_le_compat_r.
 left ; apply Rinv_0_lt_compat ; apply W_even_pos.
 apply Wallis_maj.
 rewrite <- Rinv_r with (W_even n).
 unfold Rdiv ; apply Rmult_le_compat_r.
 left ; apply Rinv_0_lt_compat ; apply W_even_pos.
 apply Wallis_odd_le_even.
 apply Rgt_not_eq ; apply W_even_pos.
Qed.

Lemma Wallis_quotient_lim1 : Rseq_cv (fun n => W_odd n / W_even n) 1.
Proof.
apply Rseq_sandwich_theorem
  with (Un := fun (n : nat) => 2 * n / (S (2 * n))) (Wn := fun n => 1).
  apply Rseq_cv_eq_compat with (Un := fun n => 2 * n / (2 * n + 1)).
    intros n.
    rewrite S_INR; rewrite mult_INR; reflexivity.
  apply Rseq_equiv_cv_div.
    apply Rseq_equiv_sym.
    apply Rseq_equiv_plus_little_O_compat_l.
    apply Rseq_equiv_refl.
    apply Rseq_little_O_Rmult_compat_r; [|intro H; fourier].
    eapply Rseq_little_O_eq_compat
      with (Un := Rseq_poly 0) (Vn := Rseq_poly 1).
      intros n; unfold Rseq_poly; apply pow_O.
      intros n; unfold Rseq_poly; apply pow_1.
      apply Rseq_poly_little_O; constructor.
    intros n; assert (H := pos_INR n); intros Hc; fourier.
  apply Rseq_constant_cv.
  apply Wallis_bound.
Qed.

Lemma Wallis_quotient :
  forall n, W_odd n / W_even n = (2*2 ^ (4*n) * (fact n) ^ 4)/(PI*fact (2 * n) * fact (S (2 * n))).
Proof.
intro n; unfold W_odd, W_even, Rsqr, Rdiv.
replace (2 ^ (4 * n)) with ((2 ^ (2 * n)) * (2 ^ (2 * n))).
field; repeat split.
  apply not_0_INR; apply fact_neq_0.
  apply not_0_INR; apply fact_neq_0.
  apply PI_neq0.
  apply not_0_INR; apply fact_neq_0.
  apply pow_nonzero; intros H; fourier.
  replace (4 * n)%nat with (2 * n + 2 * n)%nat by omega.
  rewrite pow_add.
  reflexivity.
Qed.

Lemma Rseq_cv_eq_compat1 : forall Un Vn l,
  {m | forall n, (n >= m)%nat -> Un n = Vn n} ->
    Rseq_cv Un l -> Rseq_cv Vn l.
Proof.
intros Un Vn l Heq Hseq.
intros eps Heps.
destruct Heq as (m, Heq).
destruct (Hseq eps Heps) as (N, Hseq2).
exists (max N m).
intros n Hn.
assert (Hm : (n >= m)%nat). apply le_trans with (max N m) ; intuition.
unfold R_dist in *.
rewrite <- (Heq n) ; intuition.
apply Hseq2; intuition. apply le_trans with (max N m) ; intuition.
Qed.

Lemma sqrt_id : forall n : nat, (INR n <> 0)%R -> sqrt (2 * n) / (2 * n) = /sqrt (2 * n).
Proof.
intros n H1.
rewrite <- (sqrt_sqrt (2 * n)) at 2.
field.
intro H. apply H1. apply Rmult_eq_reg_l with 2. rewrite Rmult_0_r. apply sqrt_eq_0 ; intuition.
apply Rmult_le_pos ; intuition. intro ; fourier.
apply Rmult_le_pos ; intuition.
Qed.

Lemma Rseq_equiv_eq : forall Un Vn,
  {m | forall n, (n >= m)%nat -> Un n = Vn n} -> Un ~ Vn.
Proof.
intros Un Vn Heq eps Heps.
destruct Heq as (N, Heq).
exists N. intros n HN.
unfold Rseq_constant, Rseq_minus, Rseq_plus.
rewrite (Heq n).
ring_simplify (Vn n - Vn n).
rewrite Rabs_R0. apply Rmult_le_pos.
intuition.
apply Rabs_pos.
assumption.
Qed.

Lemma DL_sqrt_1 : forall Un, Rseq_cv Un 0 -> (fun n => sqrt (1 + Un n) - 1) = o(1).
Proof.
intros Un Un0.
intros eps Heps.
destruct (Un0 eps Heps) as (N, HUn).
assert (H01 : 0 < 1) by fourier.
destruct (Un0 1 H01) as (N1, HUn1).
exists (max N N1).
intros n HNmax.
unfold Rseq_constant. rewrite Rabs_R1. rewrite Rmult_1_r.
apply Rle_trans with (Rabs (Un n)).
apply sqrt_var_maj.
unfold R_dist in *. left.
assert (HN : (n >= N1)%nat). apply le_trans with (max N N1) ; intuition.
generalize (HUn1 n HN) ; intros HU1.
rewrite Rminus_0_r in HU1. assumption.
left.
assert (HN : (n >= N)%nat). apply le_trans with (max N N1) ; intuition.
generalize (HUn n HN) ; intros HU1.
unfold R_dist in HU1.
rewrite Rminus_0_r in HU1. assumption.
Qed.

Lemma Rseq_cv_inv_INR : Rseq_cv (fun n => /INR (n + 1)) 0.
Proof.
generalize RinvN_cv. intros useful.
intros eps Heps.
destruct (useful eps Heps) as (N, H).
exists N.
intros n Hn.
generalize (H n Hn). intros H1.
unfold pos in  H1. simpl in H1.
rewrite plus_INR. apply H1.
Qed.

Lemma Rinv_plus : forall a b c : R, c <> 0 -> (a + b) / c = a / c + (b / c).
Proof.
intros a b c d.
field ; apply d.
Qed.

Lemma Rinv_eq_1 : forall a, a <> 0 -> a / a = 1.
Proof.
intros.
field ; assumption.
Qed.

Lemma pow_exp_ln : forall x n, 0 < x -> x ^ n = exp (n * ln x).
Proof.
intros x n H.
induction n.
 simpl. rewrite Rmult_0_l. rewrite exp_0. reflexivity.

 rewrite S_INR. rewrite Rmult_plus_distr_r.
 rewrite <- tech_pow_Rmult. rewrite exp_plus.
 rewrite Rmult_1_l. rewrite exp_ln ; [ | apply H].
 rewrite IHn. ring.
Qed.

Lemma Rseq_equiv_ln : forall Un, Rseq_cv Un 0 -> (fun n => ln (1 + Un n)) ~ Un.
Proof.
intros Un Hu.
destruct (Hu 1) as [M HM]; [fourier|].
apply Rseq_equiv_sym.
intros eps Heps.
assert (H1 : 0 < 1) by fourier.
destruct (Rpser_little_O_partial_sum _ Un 1 1 H1 Hu ln_plus_cv_radius eps Heps) as [N HN].
exists (Max.max M N); intros n Hn.
unfold Rseq_minus; simpl.
pattern Un at 2; replace (Un n) with (Un n ^ 1) by field.
eapply Rle_trans; [right|apply HN].
  simpl; rewrite ln_plus_taylor_sum.
  rewrite Rabs_minus_sym.
  apply f_equal; field.
  replace (Un n) with (Un n - 0) by field; apply HM.
  eapply le_trans; [apply Max.le_max_l|eassumption].
eapply le_trans; [apply Max.le_max_r|eassumption].
Qed.

Lemma Rseq_cv_inv_INR_0_1 : Rseq_cv (fun n => - / (2 * INR n + 1))%R 0%R.
Proof.
replace 0 with (-0)%R by intuition. apply Rseq_cv_opp_compat.
generalize RinvN_cv. intros useful.
intros eps Heps;
destruct (Rseq_cv_inv_INR eps Heps) as (N, Hun).
exists N. intros n Hn.
generalize (Hun n Hn) ; intros Hun1.
apply Rle_lt_trans with (/INR (n + 1))%R.
 unfold R_dist. rewrite Rminus_0_r. rewrite Rabs_pos_eq.
  apply Rle_Rinv.
   generalize (pos_INR n) ; intros ; rewrite plus_INR ; intuition.
   generalize (pos_INR n) ; intros ; fourier.
   rewrite plus_INR. simpl. apply Rplus_le_compat_r. replace (INR n)%R with ((INR n ) * 1)%R by ring. rewrite Rmult_comm. apply Rmult_le_compat.
    intuition.
    apply pos_INR.
    intuition.
    intuition.
 left. apply Rinv_0_lt_compat. generalize (pos_INR n) ; intros ; fourier.
unfold R_dist in Hun1. rewrite Rminus_0_r in Hun1. rewrite Rabs_right in Hun1. apply Hun1.
left. apply Rgt_lt. apply Rinv_0_lt_compat. intuition.
Qed.

Lemma Rseq_equiv_continuity : forall Un Vn l f, continuity_pt f l ->
  f l <> 0 -> Rseq_cv Un l -> Rseq_cv Vn l ->
    (fun n : nat => f (Un n)) ~ (fun n : nat => f (Vn n)).
Proof.
intros Un Vn l f Hcont H0 Hun Hvn.
apply Rseq_equiv_trans with (f l).
 apply Rseq_cv_equiv_constant.
  assumption.
  apply Rseq_cv_continuity_compat ; [assumption | reg].

 apply Rseq_equiv_sym.
 apply Rseq_cv_equiv_constant.
  assumption.
  apply Rseq_cv_continuity_compat ; [assumption | reg].
Qed.

Lemma Wallis_quotient_lim2 l :
l <> 0 ->
  (fun n => fact n) ~ (fun n => (n /(exp 1)) ^ n * sqrt n * l) ->
    Rseq_cv (fun n => W_odd n / W_even n) (l^2/(2*PI)).
Proof.
intros l Hneq Hl.
apply Rseq_cv_eq_compat with (fun n => (2*2 ^ (4*n) * (fact n) ^ 4)/(PI*fact (2 * n) * fact (S (2 * n)))).
  intro; rewrite Wallis_quotient; reflexivity.
assert (H2n : (fun n => fact (2 * n)) ~ (fun n => (2 * n / exp 1) ^ (2 * n) * sqrt (2 * n) * l)).
  assert (Hex : is_extractor (mult 2)).
    intros n; omega.
  pose (exist _ _ Hex) as db.
  assert (Hrw1 : extracted db fact == (fun n => fact (2 * n))).
    intros n; reflexivity.
  assert (Hrw2 : extracted db (fun n => (n / exp 1) ^ n * sqrt n * l) == (fun n => (2 * n / exp 1) ^ (2 * n) * sqrt (2 * n) * l)).
    intros n; unfold extracted. simpl.
    repeat rewrite plus_INR; simpl.
    replace (n + (n + 0)) with (2 * n) by field.
    reflexivity.
  eapply Rseq_equiv_eq_compat; [eassumption|eassumption|].
  apply Rseq_equiv_subseq_compat with (Un := fact) (Vn := fun k : nat => (k / exp 1) ^ k * sqrt k * l) (phi := db).
  assumption.
assert (H2n1 : (fun n : nat => fact (S(2*n))) ~ (fun n : nat => (S(2*n) / exp 1) ^ (S(2*n)) * sqrt (S(2*n)) * l)).
  assert (Hex : is_extractor (fun n => S (mult 2 n))).
    intros n; omega.
  pose (exist _ _ Hex) as db.
  apply Rseq_equiv_eq_compat with (extracted db fact) (extracted db (fun n : nat => (n / exp 1) ^ n * sqrt n * l)).
    intro n; unfold db; reflexivity.
    intro n; unfold db; reflexivity.
  apply Rseq_equiv_subseq_compat with (Un := fact) (Vn := fun k : nat => (k / exp 1) ^ k * sqrt k * l)(phi := db).
  assumption.

apply Rseq_equiv_cv_constant.
Open Local Scope Rseq_scope.
apply Rseq_equiv_eq_compat with
  (Un := 2 * (fun n => 2 ^ (4 * n)) * fact * fact * fact * fact *
  / (PI * (fun n => fact (2 * n)) * (fun n => fact (S (2 * n))))) (Vn := (l^2 / (2 * PI))%R).
intro n; unfold Rseq_mult, Rseq_inv, Rdiv, Rseq_constant; ring.
intro; reflexivity.
apply Rseq_equiv_trans with
(2 * (fun n => 2 ^ (4 * n)) * ((fun n : nat => ((n / exp 1) ^ n * sqrt n * l)%R))*((fun n : nat => ((n / exp 1) ^ n * sqrt n * l)%R))*((fun n : nat => ((n / exp 1) ^ n * sqrt n * l)%R))*((fun n : nat => ((n / exp 1) ^ n * sqrt n * l)%R))*
  / (PI * (fun n : nat => ((2 * n / exp 1) ^ (2 * n) * sqrt (2 * n) * l)%R) *  (fun n => (S (2 * n) / exp 1) ^ S (2 * n) * sqrt (S (2 * n)) * l)%R)).

repeat apply Rseq_equiv_mult_compat;
  try apply Rseq_equiv_refl; try assumption.
apply Rseq_equiv_inv_compat.
exists O. intros n Hn.
unfold Rseq_mult, Rseq_constant.
apply Rmult_integral_contrapositive ; split ;
[ apply Rmult_integral_contrapositive ; split ;
[ apply PI_neq0 | apply INR_fact_neq_0] | apply INR_fact_neq_0 ].

exists (S O). intros n Hn; unfold Rseq_mult, Rseq_constant.

assert(H : {m | n = S m}). exists (pred n). intuition.
destruct H as (m, Subst).
rewrite Subst.

apply Rmult_integral_contrapositive ; split ;
[ apply Rmult_integral_contrapositive ; split ;
[ apply PI_neq0 | apply Rmult_integral_contrapositive ; split ;
[ apply Rmult_integral_contrapositive ; split ;
[ (apply pow_nonzero ; unfold Rdiv ; apply Rmult_integral_contrapositive ; split ;
[ apply Rmult_integral_contrapositive ; split ;
[ (intro ; fourier) | (apply not_0_INR ; intuition) ] | (apply Rinv_neq_0_compat ; generalize (exp_pos 1) ; intros ; intro ; fourier)])
| (intro H ; apply sqrt_eq_0 in H ; [ (apply Rmult_integral in H ; destruct H as [H|H] ; [ fourier | (generalize H ; apply not_0_INR ; intuition) ])
| apply Rmult_le_pos ; intuition ]) ]
| apply Hneq ] ] | apply Rmult_integral_contrapositive ; split ;
[ apply Rmult_integral_contrapositive ; split ;
[ (apply pow_nonzero ; unfold Rdiv ; apply Rmult_integral_contrapositive ; split ;
[ apply not_0_INR ; intuition | (apply Rinv_neq_0_compat ; generalize (exp_pos 1) ; intros ; intro ; fourier) ])
| (intro H ; apply sqrt_eq_0 in H ; [ (apply not_0_INR in H ; intuition )
| apply pos_INR ])
] | apply Hneq ] ].

apply Rseq_equiv_mult_compat.
apply Rseq_equiv_mult_compat.
apply Rseq_equiv_refl.
assumption.
assumption.
clear Hl H2n H2n1.
apply Rseq_cv_equiv_constant.
unfold Rdiv.
apply Rmult_integral_contrapositive ; split .
 apply pow_nonzero; assumption.
 apply Rinv_neq_0_compat.
 apply Rmult_integral_contrapositive; split.
  intro ; fourier.
  apply PI_neq0.
unfold Rseq_constant, Rseq_mult, Rseq_div, Rseq_plus, Rseq_minus, Rseq_inv.
(* Simplification of the expression ! *)
apply Rseq_cv_eq_compat1 with
(fun n:nat => (sqrt (2 * n) / sqrt (S (2 * n))) * (Rsqr l / PI * exp 1) * ((2 * n) / S (2 * n)) ^ (S (2 * n)) * /2)%R.
exists (S O). intros n Hn.
(* need to have n <> 0 *)
assert(H : {m | n = S m}). exists (pred n). intuition.
destruct H as (m, Subst).
(* Solving the simplification equation *)
unfold Rseq_constant.
unfold Rdiv.
ring_simplify.
replace (sqrt n ^ 4) with (n ^ 2) by (rewrite <- (sqrt_sqrt n) at 1 ; [ring | intuition]).
repeat rewrite Rpow_mult_distr.
repeat rewrite <- pow_mult.
repeat rewrite Rinv_mult_distr.
replace (2 ^ (4 * n))%R with (2 ^ (2 * n) * 2 ^ (2 * n))%R by (repeat rewrite pow_mult ;
rewrite <- Rpow_mult_distr ; replace (2 ^ 2 * 2 ^ 2)%R with (2 ^ 4)%R by (unfold pow ; ring) ; ring).
rewrite (mult_comm n 4).
replace (n ^ (4 * n))%R with (n ^ (2 * n) * n ^ (2 * n))%R by
(rewrite <- pow_add ; ring_simplify (2 * n + 2 * n)%nat ; ring).
replace ((/exp 1) ^ (4 * n))%R with ((/exp 1) ^ (2 * n) * (/exp 1) ^ (2 * n))%R by
(rewrite <- pow_add ; ring_simplify (2 * n + 2 * n)%nat ; ring).
repeat rewrite Rinv_pow. rewrite Rinv_involutive.
rewrite <- tech_pow_Rmult with (exp 1) (2 * n)%nat.
rewrite <- (Rinv_pow (exp 1) _).
do 2 rewrite <- tech_pow_Rmult.
replace (/sqrt (2 * n))%R with (sqrt (2 * n) / (2 * n))%R by (apply sqrt_id ; inversion Hn ; [intuition | apply not_0_INR ; intuition]).
unfold Rsqr. unfold Rdiv.
repeat rewrite <- Rinv_pow. field.

(* Solving the "<> 0" equations *)

split. assumption.
split. intro H. apply sqrt_eq_0 in H. apply not_0_INR in H. assumption.
intuition.
intuition.
split. apply pow_nonzero. apply not_0_INR. intuition.
split. apply not_0_INR. intuition.
split. apply pow_nonzero. apply not_0_INR. intuition.
split. apply pow_nonzero. intro. fourier.
split. apply PI_neq0.
apply pow_nonzero. intro. generalize (exp_pos 1) ; intros ; fourier.
apply not_0_INR. intuition.
intro ; fourier.
apply not_0_INR. intuition.
intro. generalize (exp_pos 1) ; intros ; fourier.
intro. generalize (exp_pos 1) ; intros ; fourier.
apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
apply not_0_INR ; intuition.
apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
apply not_0_INR ; intuition.
intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply not_0_INR ; intuition.
apply pos_INR.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply not_0_INR ; intuition.
apply pos_INR.
apply Hneq.
apply pow_nonzero. intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply Rmult_integral_contrapositive ; split.
intro ; fourier.
apply not_0_INR. intuition.
apply Rmult_le_pos. intuition.
intuition.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply Rmult_integral_contrapositive ; split.
intro ; fourier.
apply not_0_INR. intuition.
apply Rmult_le_pos. intuition. intuition.
assumption.
apply PI_neq0.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply Rmult_integral_contrapositive ; split.
intro ; fourier.
apply not_0_INR. intuition.
apply Rmult_le_pos. intuition. intuition.
assumption.
apply Rmult_integral_contrapositive ; split.
apply PI_neq0.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. intro ; fourier.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply Rmult_integral_contrapositive ; split.
intro ; fourier.
apply not_0_INR. intuition.
apply Rmult_le_pos. intuition. intuition.
assumption.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply Rmult_integral_contrapositive ; split.
apply pow_nonzero. apply not_0_INR ; intuition.
apply pow_nonzero. apply Rinv_neq_0_compat. intro. generalize (exp_pos 1) ; intros ; fourier.
intro H. apply sqrt_eq_0 in H. generalize H. apply not_0_INR. intuition.
apply pos_INR. assumption.
(* End of "<> 0" *)

(* end of simplification *)

eapply Rseq_equiv_cv_compat.
  2: reflexivity.
 symmetry; instantiate (1 := (1 * (Rsqr l / PI * exp 1) * /exp 1 * /2)).
 apply Rseq_equiv_mult_compat ; [ | apply Rseq_equiv_eq ; exists O ; intro ; intuition ].
 apply Rseq_equiv_mult_compat.
  apply Rseq_equiv_mult_compat ; [ | apply Rseq_equiv_eq ; exists O ; intro ; intuition ].
  pose (Un := (fun n => - / (2 * INR n + 1))%R).
  assert (Hcv0 : Rseq_cv Un 0). unfold Un.
   apply Rseq_cv_inv_INR_0_1.
  intros eps Heps.
  destruct (DL_sqrt_1 Un Hcv0 eps Heps) as (N, DL).
  exists N.
  intros n HN. generalize (DL n HN). intros DL1.
  unfold Rseq_constant, Rseq_plus, Rseq_minus, Rseq_mult, Rseq_inv, Un in *.
  rewrite <-sqrt_div ; [ | (apply Rmult_le_pos ; intuition) | intuition ].
  rewrite S_INR. rewrite mult_INR. do 2 rewrite S_INR.
  replace (2 * n)%R with ((2 * n + 1) - 1)%R by ring. rewrite Rplus_0_l.
  unfold Rminus. rewrite Rinv_plus.
   rewrite Rinv_eq_1. unfold Rdiv in *. rewrite Ropp_mult_distr_l_reverse. rewrite Rmult_1_l. rewrite Rabs_minus_sym in DL1. apply DL1.
    generalize (pos_INR n) ; intuition ; fourier.
    generalize (pos_INR n) ; intuition ; fourier.
 apply Rseq_equiv_trans with (fun n : nat => exp ((2 * n + 1) * ln (1 - /(2 * INR n + 1)))).
 pose (Un := (fun n:nat => - /(2 * (INR n) + 1))%R).
 assert (Hcv0 : Rseq_cv Un 0).
  apply Rseq_cv_inv_INR_0_1.
 apply Rseq_equiv_trans with (fun n : nat => exp ((2 * n + 1) * (- / (2 * n + 1)))).
  apply Rseq_equiv_eq. exists O. intros n Hn.
  field_simplify ((2 * n + 1) * - / (2 * n + 1))%R. unfold Rseq_constant, Rseq_inv.
   rewrite <- exp_Ropp. unfold Rdiv. rewrite Rinv_1. rewrite Rmult_1_r. reflexivity.
   generalize (pos_INR n) ; intuition ; fourier.
  apply Rseq_equiv_continuity with ((-1)).
   reg.
   generalize (exp_pos (-1)) ; intros l1 H1 ; fourier.
   apply Rseq_cv_eq_compat with (-R1).
   intros n. unfold Rseq_constant, Rseq_minus, Rseq_plus, Rseq_opp.
   field. intros H1 ; generalize (pos_INR n) ; intros ; fourier.
   intuition.
   eapply Rseq_equiv_cv_compat.
     2: reflexivity.
    symmetry; instantiate (1 := (fun n => (2 * INR n + 1) * - / (2 * INR n + 1))%R).
    apply Rseq_equiv_mult_compat.
     apply Rseq_equiv_refl.
     apply Rseq_equiv_sym. apply Rseq_equiv_ln.
     apply Hcv0.
    apply Rseq_cv_eq_compat with (-R1).
    unfold Rseq_opp, Rseq_constant, Rseq_minus.
    intros n. field. generalize (pos_INR n) ; intuition ; fourier.
    intuition.
 apply Rseq_equiv_eq.
 exists 1%nat.
 intros n Hn.
 rewrite pow_exp_ln.
  replace (2 * n)%R with ((2 * n + 1) - 1)%R by ring.
  rewrite S_INR. rewrite mult_INR. repeat rewrite S_INR.
  rewrite Rplus_0_l.
  unfold Rminus. rewrite Rinv_plus. rewrite Rinv_eq_1.
   ring_simplify (2 * n + 1 + -1 + 1)%R. unfold Rdiv. ring_simplify (1 + -1 * / (2 * n + 1))%R.
   rewrite (Rplus_comm (- / (2 * n + 1)) _). reflexivity.
   generalize (pos_INR n) ; intuition ; fourier.
   generalize (pos_INR n) ; intuition ; fourier.
   unfold Rdiv. apply Rmult_lt_0_compat. apply Rmult_lt_0_compat ; intuition.
   apply Rinv_0_lt_compat.
   rewrite S_INR. generalize (pos_INR n) ; intuition ; fourier.
apply Rseq_cv_eq_compat with (Rsqr l / (2 * PI)).
intro. unfold Rseq_mult, Rseq_plus, Rseq_constant, Rseq_div, Rseq_inv.
field. split.
intro. generalize (exp_pos 1) ; intros ; fourier.
apply PI_neq0.
intros eps Heps. exists O. intros n Hn.
unfold R_dist. unfold Rseq_constant, Rseq_div, Rseq_inv, Rseq_mult.
unfold Rdiv, Rsqr. rewrite Rminus_diag_eq.
rewrite Rabs_R0. apply Heps.
ring.
Qed.

End Wallis.


text:
Inductive list (A : Set) : Set :=
  | Nil : list A
  | Cons : A -> list A -> list A
.

Check Nil.
Check Cons.

Implicit Arguments Nil [A].
Implicit Arguments Cons [A].

Infix "::" := Cons (at level 60, right associativity).

Check (@Nil bool).

Check 3 :: 4 :: Nil.

Fixpoint length  {A : Set} (xs : list A) :=
  match xs with
  | Nil => 0
  | x :: xs' => 1 + length xs'
end.

Definition xs := 3 :: 4 ::1 :: Nil.

Eval compute in length xs.

Fixpoint app {A : Set} (xs ys : list A) :=
  match xs with
  | Nil => ys
  | x :: xs' => x :: app xs' ys
end.

Infix "++" := app (at level 60, right associativity).

Fixpoint rev {A : Set} (xs : list A) :=
  match xs with
  | Nil => Nil
  | x :: xs' => (rev  xs') ++ (x :: Nil)
end.

Eval compute in rev xs.

Fixpoint rev' (A : Set) (xs : list A) :=
  match xs with
  | Nil => Nil
  | x :: xs' => (rev' A xs') ++ (x :: Nil)
  end.

Eval compute in rev' nat xs.

Fixpoint sum (ns : list nat) :=
  match ns with
  | Nil => 0
  | n :: ns' => n + sum ns'
end.

Eval compute in sum xs.

Print Nil.

Definition propo := forall xs, length xs = 0 -> xs = (@Nil nat).
Eval compute in propo.

Theorem app_length : forall {A : Set} (xs ys : list A),
    length (xs ++ ys) = length xs + length ys.
Proof.
  intro A.
  intro xs.
  intro ys.
  induction xs.
    simpl.
    reflexivity.

    simpl.
    rewrite <- IHxs.
    reflexivity.
Qed.

Require Import Arith.

Theorem rev_length : forall {A : Set} (xs : list A),
    length (rev xs) = length xs.
Proof.
  intro A.
  intro xs.
  induction xs.
    simpl.
    reflexivity.

    simpl.
    rewrite app_length.
    simpl.
    rewrite plus_comm.
    rewrite IHxs.
    simpl.
    reflexivity.
Qed.

Theorem app_nil : forall {A : Set} (xs : list A), xs ++ Nil = xs.
Proof.
  intro A.
  intro xs.
  induction xs.
    simpl.
    reflexivity.

    simpl.
    rewrite IHxs.
    reflexivity.
Qed.

Theorem  app_assoc : forall {A : Set} (xs ys zs: list A),
  xs ++ (ys ++ zs) = (xs ++ ys) ++ zs.
Proof.
  intro A.
  intro xs.
  intro ys.
  intro zs.
  induction xs.
    simpl.
    reflexivity.

    simpl.
    rewrite IHxs.
    reflexivity.
Qed.

Theorem app_rev : forall {A : Set} (xs ys: list A),
    rev (xs ++ ys) = rev ys ++ rev xs.
Proof.
  intro A.
  intro xs.
  intro ys.
  induction xs.
    simpl.
    rewrite app_nil.
    reflexivity.

    simpl.
    rewrite IHxs.
    rewrite app_assoc.
    reflexivity.
Qed.

Theorem rev_rev : forall {A : Set} (xs : list A),
    rev (rev xs) = xs.

  intros A xs.
  induction xs.
    simpl.
    reflexivity.

    simpl.
    rewrite app_rev.
    simpl.
    rewrite IHxs.
    reflexivity.
Qed.

(* 2nd day *)

Require Import Recdef.

Function foo (xs : list nat) {measure length xs} : list nat :=
  match xs with
  | Nil => Nil
  | x :: xs' => x :: foo (rev xs')
  end.

  intro xs.
  intro x.
  intro xs'.
  intro H.
  rewrite rev_length.
  simpl.
  SearchPattern(_ < S _).
  apply lt_n_Sn.
Qed.

author: seizans

repo: coqtest

save_path: github-repos/coq/seizans-coqtest


text:
(************************************************************************)
(* Copyright 2006 Milad Niqui                                           *)
(* This file is distributed under the terms of the                      *)
(* GNU Lesser General Public License Version 2.1                        *)
(* A copy of the license can be found at                                *)
(*                  <http://www.gnu.org/licenses>                       *)
(************************************************************************)

Require Import digits.
 Require Import R_addenda.
Require Import Fourier_solvable_ineqs.
Require Import Fourier.
Require Import Bounded_M.

(** Properties of the predicate [Refining_M]. *)

Open Scope Q_scope.

Lemma Bounded_M_auxiliary_5: forall c d r, Zero <= d + c -> Zero <= d - c ->(-1<=r<=1)%R->(0<=c*r+d)%R.
Proof.
 intros c d r H1 H2 Hr.
 (* TP : Zero<=d *)
 assert (Hd: Zero<=d);
 [apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto; stepl Zero; auto; stepr ((d+c)+(d-c)); auto; ring|].

 destruct (Qtrichotomy_inf Zero c) as [[Hc|Hc]|Hc]; realify_Q.
   intros; apply Bounded_M_pos_auxiliary_3; trivial.
   intros Hc _ _ Hd; rewrite <- Hc; stepr d; trivial; ring.
   intros; apply Bounded_M_pos_auxiliary_4; trivial.
Qed.

Lemma Bounded_M_auxiliary_6: forall c d r, d + c<=Zero -> d - c<=Zero ->(-1<=r<=1)%R->(c*r+d<=0)%R.
Proof.
 intros c d r H1 H2 Hr.
 (* TP : d<=Zero *)
 assert (Hd: d<=Zero);
 [apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto; stepr Zero; auto; stepl ((d+c)+(d-c)); auto; ring|].

 destruct (Qtrichotomy_inf Zero c) as [[Hc|Hc]|Hc]; realify_Q.
   intros; apply Bounded_M_neg_auxiliary_3; trivial.
   intros Hc _ _ Hd; rewrite <- Hc; stepl d; trivial; ring.
   intros; apply Bounded_M_neg_auxiliary_4; trivial.
Qed.

Lemma Is_refining_M_auxiliary_1: forall (a b c d:Q), Zero <= a + b + c + d -> Zero <= a - b - c + d ->
                                                     Zero <= - a - b + c + d -> Zero <= - a + b - c + d ->
         Zero<=d/\Zero<=d-c/\Zero<=d+c/\Zero<=d-b/\Zero<=d+b/\Zero<=d-a/\Zero<=d+a.
Proof.
 let
 t_local:= apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto; stepl Zero; auto
 in
 ( intros a b c d H1 H2 H3 H4; repeat split;
 [ apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone+Qone); auto; stepl Zero; auto;
   stepr ((a + b + c + d)+(- a - b + c + d)+(a - b - c + d)+(- a + b - c + d))
 | t_local; stepr ((a - b - c + d)+(- a + b - c + d))
 | t_local; stepr ((a + b + c + d)+(- a - b + c + d))
 | t_local; stepr ((a - b - c + d)+(- a - b + c + d))
 | t_local; stepr ((a + b + c + d)+(- a + b - c + d))
 | t_local; stepr ((- a - b + c + d)+(- a + b - c + d))
 | t_local; stepr ((a + b + c + d)+(a - b - c + d))
 ];
 auto; ring).
Qed.

Lemma Is_refining_M_auxiliary_2: forall (a b c d:Q), a + b + c + d <= Zero -> a - b - c + d <= Zero ->
                                                     - a - b + c + d <= Zero -> - a + b - c + d <= Zero ->
                 d<=Zero/\d-c<=Zero/\d+c<=Zero/\d-b<=Zero/\d+b<=Zero/\d-a<=Zero/\d+a<=Zero.
Proof.
 let
 t_local:= apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto; stepr Zero; auto
 in
 ( intros a b c d H1 H2 H3 H4; repeat split;
 [ apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone+Qone); auto; stepr Zero; auto;
   stepl ((a + b + c + d)+(- a - b + c + d)+(a - b - c + d)+(- a + b - c + d))
 | t_local; stepl ((a - b - c + d)+(- a + b - c + d))
 | t_local; stepl ((a + b + c + d)+(- a - b + c + d))
 | t_local; stepl ((a - b - c + d)+(- a - b + c + d))
 | t_local; stepl ((a + b + c + d)+(- a + b - c + d))
 | t_local; stepl ((- a - b + c + d)+(- a + b - c + d))
 | t_local; stepl ((a + b + c + d)+(a - b - c + d))
 ];
 auto; ring).
Qed.

Lemma Is_refining_M_d_minusplus_c_sign_dec: forall a b c d, Is_refining_M ((a,b),(c,d))->
              (Qle Zero (Qplus (Qplus (Qplus a b) c) d) /\ Qle Zero (Qplus (Qminus (Qminus a b) c) d) /\
               Qle Zero (Qplus (Qplus (Qminus (Qopp a) b) c) d) /\ Qle Zero (Qplus (Qminus (Qplus (Qopp a) b) c) d) /\
               Qlt Zero (Qminus d c)/\Qlt Zero (Qplus d c)/\Qlt Zero d) \/
              (Qle (Qplus (Qplus (Qplus a b) c) d) Zero /\ Qle (Qplus (Qminus (Qminus a b) c) d) Zero /\
	       Qle (Qplus (Qplus (Qminus (Qopp a) b) c) d) Zero /\ Qle (Qplus (Qminus (Qplus (Qopp a) b) c) d) Zero /\
               Qlt (Qminus d c) Zero/\Qlt (Qplus d c) Zero/\Qlt d Zero).
Proof.
 intros a b c d [[[H0 H0']|[H0 H0']] [[H1 [H2 [H3 H4]]]|[H1 [H2 [H3 H4]]]]]; unfold fst, snd in H0, H0', H1, H2, H3, H4;
 first [generalize (Bounded_M_auxiliary_1 _ _ H0 H0')|generalize (Bounded_M_auxiliary_2 _ _ H0 H0')] ;intros Hd;
 tauto ||
 first [decompose record (Is_refining_M_auxiliary_1 _ _ _ _ H1 H2 H3 H4)
       |decompose record (Is_refining_M_auxiliary_2 _ _ _ _ H1 H2 H3 H4)]; contradiction.
Qed.

Close Scope Q_scope.

Lemma Is_refining_M_denom_nonvanishing_M:forall mu r, Is_refining_M mu -> ((-1)<=r<=1)%R -> denom_nonvanishing_M mu r.
Proof.
 intros ((a,b),(c,d)) r H_refining Hr;
 generalize (Is_refining_M_d_minusplus_c_sign_dec _ _ _ _ H_refining); clear H_refining.
 intros [[_[_[_[_[H1 [H2 H3]]]]]]|[_[_[_[_[H1 [H2 H3]]]]]]]; red; simpl; realify_Q; intros H1 H2 H3;
 destruct (Qtrichotomy_inf Zero c) as [[Hc|Hc]|Hc]; realify_Q; intros Hc.
  apply Rlt_not_eq'; apply Bounded_M_pos_auxiliary_1; trivial...
  rewrite <- Hc; stepl d; auto; ring...
  apply Rlt_not_eq'; apply Bounded_M_pos_auxiliary_2; trivial...
  apply Rlt_not_eq; apply Bounded_M_neg_auxiliary_1; trivial...
  rewrite <- Hc; stepl d; auto; ring...
  apply Rlt_not_eq; apply Bounded_M_neg_auxiliary_2; trivial...
Qed.

Lemma Is_refining_M_property:forall mu r, ((-1)<=r<=1)%R->Is_refining_M mu -> ((-1)<=as_Moebius mu r<=1)%R.
Proof.
 intros ((a,b),(c,d)) r Hr H_refining;
 generalize (Is_refining_M_denom_nonvanishing_M _ _ H_refining Hr).
 unfold denom_nonvanishing_M, as_Moebius, fst, snd.
 destruct (Is_refining_M_d_minusplus_c_sign_dec _ _ _ _ H_refining)
     as [[H1 [H2 [H3 [H4 [H7 [H5 H6]]]]]]|[H1 [H2 [H3 [H4 [H7 [H5 H6]]]]]]];
 [ generalize (Is_refining_M_auxiliary_1 _ _ _ _ H1 H2 H3 H4)
 | generalize (Is_refining_M_auxiliary_2 _ _ _ _ H1 H2 H3 H4)
 ];
 intros [H12 [H13 [H14 [H8 [H9 [H10 H11]]]]]];
 realify_Q; intros H1 H2 H3 H4 H7 H5 H6 H12 H13 H14 H8 H9 H10 H11 H_denom.
  (* everything is positive *)
   (* TP: 0 < c*r+d *)
   assert (H_crd:0<c*r+d);
   [ destruct (Qtrichotomy_inf Zero c) as [[Hc|Hc]|Hc]; realify_Q; intros Hc;
     [ (* 0<c *) apply Bounded_M_pos_auxiliary_1; trivial
     | (* c=0 *) rewrite <- Hc; stepr d; auto; ring
     | (* c<0 *) apply Bounded_M_pos_auxiliary_2; trivial]
   |].

  split.

   (* lower bound *)
   destruct (Qtrichotomy_inf Zero (Qplus c a)) as [[Hca|Hca]|Hca]; realify_Q; intros Hca;
   (stepl (-1/1); [|field; auto]); apply Rmult_Rdiv_pos_Rle; auto; apply Rle_zero_Rminus; (stepr ((c+a)*r+(d+b)); [|ring]).
    (* 0<c+a *)
    apply Bounded_M_pos_auxiliary_3; trivial; ring_exact_R H4.
    (* c+a=0 *)
    rewrite <- Hca; stepr (d+b); auto; ring.
    (* 0<c+a *)
    apply Bounded_M_pos_auxiliary_4; trivial; ring_exact_R H1...
   (* upper bound *)
   destruct (Qtrichotomy_inf Zero (Qminus c a)) as [[Hca|Hca]|Hca]; realify_Q; intros Hca;
   (stepr (1/1); [|field; auto]); apply Rmult_Rdiv_pos_Rle; auto; apply Rle_zero_Rminus; (stepr ((c-a)*r+(d-b)); [|ring]).
    (* 0<c+a *)
    apply Bounded_M_pos_auxiliary_3; trivial; ring_exact_R H2.
    (* c+a=0 *)
    rewrite <- Hca; stepr (d-b); auto; ring.
    (* 0<c+a *)
    apply Bounded_M_pos_auxiliary_4; trivial; ring_exact_R H3...
  (* everything is negative *)
   (* TP: c*r+d < 0 *)
   assert (H_crd:c*r+d<0);
   [ destruct (Qtrichotomy_inf Zero c) as [[Hc|Hc]|Hc]; realify_Q; intros Hc;
     [ (* 0<c *) apply Bounded_M_neg_auxiliary_1; trivial
     | (* c=0 *) rewrite <- Hc; stepl d; auto; ring
     | (* c<0 *) apply Bounded_M_neg_auxiliary_2; trivial]
   |].

  split.

   (* lower bound *)
   destruct (Qtrichotomy_inf Zero (Qplus c a)) as [[Hca|Hca]|Hca]; realify_Q; intros Hca;
   (stepl (1/(-1)); [|field; auto]); apply Rmult_Rdiv_neg_Rle; auto; apply Rminus_le; (stepl ((c+a)*r+(d+b)); [|ring]).
    (* 0<c+a *)
    apply Bounded_M_neg_auxiliary_3; trivial; ring_exact_R H1.
    (* c+a=0 *)
    rewrite <- Hca; stepl (d+b); auto; ring.
    (* 0<c+a *)
    apply Bounded_M_neg_auxiliary_4; trivial; ring_exact_R H4...
   (* upper bound *)
   destruct (Qtrichotomy_inf Zero (Qminus c a)) as [[Hca|Hca]|Hca]; realify_Q; intros Hca;
   (stepr (-1/-1); [|field; auto]); apply Rmult_Rdiv_neg_Rle; auto; apply Rminus_le; (stepl ((c-a)*r+(d-b)); [|ring]).
    (* 0<c+a *)
    apply Bounded_M_neg_auxiliary_3; trivial; ring_exact_R H3.
    (* c+a=0 *)
    rewrite <- Hca; stepl (d-b); auto; ring.
    (* 0<c+a *)
    apply Bounded_M_neg_auxiliary_4; trivial; ring_exact_R H2...
Qed.

Lemma denom_nonvanishing_M_product:forall mu mu' r,Is_refining_M mu ->((-1)<=as_Moebius mu' r<=1)%R->
                denom_nonvanishing_M mu' r-> denom_nonvanishing_M (product mu mu') r.
Proof.
 intros ((a,b),(c,d)) ((a',b'),(c',d')) r.
 unfold as_Moebius, denom_nonvanishing_M, product, fst, snd.
 intros  H_refining Hr H_denom.
 repeat rewrite Q_to_Rplus; repeat rewrite Q_to_Rmult.
 stepl (c*(a'*r+b')+d*(c'*r+d')); [|ring].
 apply Rlinear_non_zero_2; auto.
 apply (Is_refining_M_denom_nonvanishing_M _ _ H_refining Hr).
Qed.

Lemma Is_refining_M_Bounded_M_product:forall mu mu',Is_refining_M mu' -> Bounded_M mu ->
                   Bounded_M (product mu mu').
Proof.
 intros ((a,b),(c,d)) ((a',b'),(c',d')) H_refining' H_bounded.
 apply denom_nonvanishing_M_Bounded_M; intros r Hr.
 assert (H_mu':=Is_refining_M_property _ r Hr H_refining').
 generalize (Bounded_M_denom_nonvanishing_M _ H_bounded _ H_mu').
 unfold as_Moebius, denom_nonvanishing_M, product, fst, snd.
 intros H_denom.
 repeat rewrite Q_to_Rplus; repeat rewrite Q_to_Rmult.
 stepl (c*(a'*r+b')+d*(c'*r+d')); [|ring].
 apply Rlinear_non_zero_2; auto.
 apply (Is_refining_M_denom_nonvanishing_M _ _ H_refining' Hr).
Qed.

Lemma as_Moebius_product : forall mu mu' r, denom_nonvanishing_M mu' r -> ((-1)<=as_Moebius mu' r<=1)%R ->
                     Is_refining_M mu ->  as_Moebius (product mu mu') r = as_Moebius mu (as_Moebius mu' r).
Proof.
 intros ((a,b),(c,d)) ((a',b'),(c',d')) r H_denom H_mu'_r H_refining.
 generalize (denom_nonvanishing_M_product _ _ _ H_refining H_mu'_r H_denom).
 generalize H_denom; clear H_denom.
 unfold as_Moebius, denom_nonvanishing_M, product, fst, snd.
 realify_Q.
 intros H_denom H_denom_prod.
 repeat rewrite Rdiv_Rmult_numerator; trivial.
 repeat rewrite Rdiv_Rplus_Rmult; trivial.
 rewrite Rdiv_Rdiv_simplify; trivial.
  apply (f_equal2 Rdiv); ring.
  stepl ((c * a' + d * c') * r + (c * b' + d * d')); trivial; ring.
Qed.

Open Scope Z_scope.
Open Scope Q_scope.

Lemma Incl_M_absorbs_Is_refining_M_L:forall mu, Incl_M mu LL -> Is_refining_M (product inv_LL mu).
Proof.
 intros ((a,b),(c,d));
 unfold Incl_M, Is_refining_M, Bounded_M, product, inv_digit, map_digits, inv_LL, fst, snd.
 replace ((-1)/2 - 1/2) with (Qopp Qone); [|qZ_numerals; field; auto];
 replace (1/2 + 3/2) with (Qone +Qone); [|qZ_numerals; field; auto];
 replace (3/2 - 1/2) with Qone; [|qZ_numerals; field; auto];
 replace (1/2 + -1/2) with Zero; [|qZ_numerals; field; auto];
 repeat rewrite Qmult_one_right;
 intros [[[H1 H2]|[H1 H2]] [H3 [H4 [H5 H6]]]].

  (* Zero < d-c,d+c *)
  generalize (Qmult_resp_Qle_pos_l _ _ _ H1 H5); generalize (Qmult_resp_Qle_pos_l _ _ _ H1 H6);
  generalize (Qmult_resp_Qle_pos_l _ _ _ H2 H3); generalize (Qmult_resp_Qle_pos_l _ _ _ H2 H4); clear H3 H4 H5 H6;
  repeat rewrite Qmult_zero_right;
  intros H3 H4 H5 H6;
  split; left; repeat split.
   (* TP: boundedness *)
   stepr ((1/2)*((d+c)-(a+b))); [|qZ_numerals; field; auto];
    apply Qlt_mult_pos_pos; auto; apply Qlt_Qminus_Zero; apply Qle_lt_trans with Zero; trivial;
    apply Qmult_resp_Qle_pos_r with (Qone+Qone);  auto.
   stepr ((1/2)*((d-c)-(b-a))); [|qZ_numerals; field; auto];
    apply Qlt_mult_pos_pos; auto; apply Qlt_Qminus_Zero; apply Qle_lt_trans with Zero; trivial;
    apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto.
   (* TP: refining *)
   stepr ((c+d)+(a+b)); [|qZ_numerals; field; auto];
    stepl ((c+d)+(c + d)*(Qopp Qone)); [|ring]; apply Qle_plus_plus; auto.
   stepr (Zero-((b-a)*(Qone+Qone))); [|qZ_numerals; field; auto];
    apply Qle_Qminus_Zero; trivial.
   stepr (Zero-((a+b)*(Qone+Qone))); [|qZ_numerals; field; auto];
    apply Qle_Qminus_Zero; trivial.
   stepr ((d-c)+(b-a)); [|qZ_numerals; field; auto];
    stepl ((d-c)+(d-c)*(Qopp Qone)); [|ring]; apply Qle_plus_plus; auto.
  (* d-c,d+c < Zero *)
  generalize (Qmult_resp_Qle_neg_l _ _ _ H1 H5); generalize (Qmult_resp_Qle_neg_l _ _ _ H1 H6);
  generalize (Qmult_resp_Qle_neg_l _ _ _ H2 H3); generalize (Qmult_resp_Qle_neg_l _ _ _ H2 H4); clear H3 H4 H5 H6;
  repeat rewrite Qmult_zero_right;
  intros H3 H4 H5 H6;
  split; right; repeat split.
   (* TP: boundedness *)
   stepl (((d+c)-(a+b))*(1/2)); [|qZ_numerals; field; auto];
    apply Qlt_mult_neg_pos; auto; apply Qlt_Qminus_Zero_neg; apply Qlt_le_trans with Zero; trivial;
    apply Qmult_resp_Qle_pos_r with (Qone+Qone);  auto.
   stepl (((d-c)-(b-a))*(1/2)); [|qZ_numerals; field; auto];
    apply Qlt_mult_neg_pos; auto; apply Qlt_Qminus_Zero_neg; apply Qlt_le_trans with Zero; trivial;
    apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto.
   (* TP: refining *)
   stepl ((c+d)+(a+b)); [|qZ_numerals; field; auto];
    stepr ((c+d)+(c + d)*(Qopp Qone)); [|ring]; apply Qle_plus_plus; auto.
   stepl (Zero-((b-a)*(Qone+Qone))); [|qZ_numerals; field; auto];
    apply Qle_Qminus_Zero_neg; trivial.
   stepl (Zero-((a+b)*(Qone+Qone))); [|qZ_numerals; field; auto];
    apply Qle_Qminus_Zero_neg; trivial.
   stepl ((d-c)+(b-a)); [|qZ_numerals; field; auto];
    stepr ((d-c)+(d-c)*(Qopp Qone)); [|ring]; apply Qle_plus_plus; auto.
Qed.

Lemma Incl_M_absorbs_Is_refining_M_R:forall mu, Incl_M mu RR -> Is_refining_M (product inv_RR mu).
Proof.
 intros ((a,b),(c,d));
 unfold Incl_M, Is_refining_M, Bounded_M, product, inv_digit, map_digits, inv_RR, fst, snd.
 replace (1/2 - 1/2) with Zero; [|qZ_numerals; field; auto];
 replace (-1/2 + 3/2) with Qone; [|qZ_numerals; field; auto];
 replace (3/2 - -1/2) with (Qone +Qone); [|qZ_numerals; field; auto];
 replace (1/2 + 1/2) with Qone; [|qZ_numerals; field; auto];
 repeat rewrite Qmult_one_right;
 intros [[[H1 H2]|[H1 H2]] [H3 [H4 [H5 H6]]]].
  (* Zero < d-c,d+c *)
  generalize (Qmult_resp_Qle_pos_l _ _ _ H1 H5); generalize (Qmult_resp_Qle_pos_l _ _ _ H1 H6);
  generalize (Qmult_resp_Qle_pos_l _ _ _ H2 H3); generalize (Qmult_resp_Qle_pos_l _ _ _ H2 H4); clear H3 H4 H5 H6;
  repeat rewrite Qmult_zero_right;
  intros H3 H4 H5 H6;
  split; left; repeat split.
   (* TP: boundedness *)
   stepr ((1/2)*((d+c)+(a+b))); [|qZ_numerals; field; auto];
   apply Qlt_mult_pos_pos; auto; apply Qlt_le_reg_pos; trivial; apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto.
   stepr ((1/2)*((d-c)+(b-a))); [|qZ_numerals; field; auto];
   apply Qlt_mult_pos_pos; auto; apply Qlt_le_reg_pos; trivial; apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto.
   (* TP: refining *)
   stepr ((a+b)*(Qone+Qone)); trivial; qZ_numerals; field; auto.
   stepr ((d-c)-(b-a)); [|qZ_numerals; field; auto]; apply Qle_Qminus_Zero; trivial.
   stepr ((c+d)-(a+b)); [|qZ_numerals; field; auto]; apply Qle_Qminus_Zero; trivial.
   stepr ((b-a)*(Qone+Qone)); trivial; qZ_numerals; field; auto.
  (* d-c,d+c < Zero *)
  generalize (Qmult_resp_Qle_neg_l _ _ _ H1 H5); generalize (Qmult_resp_Qle_neg_l _ _ _ H1 H6);
  generalize (Qmult_resp_Qle_neg_l _ _ _ H2 H3); generalize (Qmult_resp_Qle_neg_l _ _ _ H2 H4); clear H3 H4 H5 H6;
  repeat rewrite Qmult_zero_right;
  intros H3 H4 H5 H6;
  split; right; repeat split.
   (* TP: boundedness *)
   stepl (((d+c)+(a+b))*(1/2)); [|qZ_numerals; field; auto];
   apply Qlt_mult_neg_pos; auto; apply Qlt_le_reg_neg; trivial; apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto.
   stepl (((d-c)+(b-a))*(1/2)); [|qZ_numerals; field; auto];
   apply Qlt_mult_neg_pos; auto; apply Qlt_le_reg_neg; trivial; apply Qmult_resp_Qle_pos_r with (Qone+Qone); auto.
   (* TP: refining *)
   stepl ((a+b)*(Qone+Qone)); trivial; qZ_numerals; field; auto.
   stepl ((d-c)-(b-a)); [|qZ_numerals; field; auto]; apply Qle_Qminus_Zero_neg; trivial.
   stepl ((c+d)-(a+b)); [|qZ_numerals; field; auto]; apply Qle_Qminus_Zero_neg; trivial.
   stepl ((b-a)*(Qone+Qone)); trivial; qZ_numerals; field; auto.
Qed.

Lemma Incl_M_absorbs_Is_refining_M_M:forall mu, Incl_M mu MM -> Is_refining_M (product inv_MM mu).
Proof.
 intros ((a,b),(c,d));
 unfold Incl_M, Is_refining_M, Bounded_M, product, inv_digit, map_digits, inv_MM, fst, snd.
 replace (0/1 - 1/1) with (Qopp Qone); [|qZ_numerals; field; auto];
 replace (1/1 + 0/1) with Qone; [|qZ_numerals; field; auto];
 replace (3/1 - 0/1) with (Qone+Qone+Qone); [|qZ_numerals; field; auto];
 replace (0/1 + 3/1) with (Qone+Qone+Qone); [|qZ_numerals; field; auto];
 repeat rewrite Qmult_one_right;
 intros [[[H1 H2]|[H1 H2]] [H3 [H4 [H5 H6]]]].
  (* Zero < d-c,d+c *)
  generalize (Qmult_resp_Qle_pos_l _ _ _ H1 H5); generalize (Qmult_resp_Qle_pos_l _ _ _ H1 H6);
  generalize (Qmult_resp_Qle_pos_l _ _ _ H2 H3); generalize (Qmult_resp_Qle_pos_l _ _ _ H2 H4); clear H3 H4 H5 H6;
  intros H3 H4 H5 H6;
  split; left; repeat split.
   (* TP: boundedness *)
   stepr ((1/3)*(d+c)); [|qZ_numerals; field; auto]; apply Qlt_mult_pos_pos; auto.
   stepr ((1/3)*(d-c)); [|qZ_numerals; field; auto]; apply Qlt_mult_pos_pos; auto.
   (* TP: refining *)
   stepr ((a+b)+(1/3)*(d+c)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepl Zero; trivial;
    stepr (((a+b)*(Qone+Qone+Qone))-((c + d)*-Qone)); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero; trivial.
   stepr ((a-b)+(1/3)*(d-c)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepl Zero; trivial;
    stepr ((d-c)-((b-a)*(Qone+Qone+Qone))); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero; trivial.
   stepr ((1/3)*(c+d)-(a+b)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepl Zero; trivial;
    stepr ((c+d)-((a+b)*(Qone+Qone+Qone))); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero; trivial.
   stepr ((1/3)*(d-c)+(b-a)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepl Zero; trivial;
    stepr (((b-a)*(Qone+Qone+Qone))-((d-c)*-Qone)); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero; trivial.
  (* d-c,d+c < Zero *)
  generalize (Qmult_resp_Qle_neg_l _ _ _ H1 H5); generalize (Qmult_resp_Qle_neg_l _ _ _ H1 H6);
  generalize (Qmult_resp_Qle_neg_l _ _ _ H2 H3); generalize (Qmult_resp_Qle_neg_l _ _ _ H2 H4); clear H3 H4 H5 H6;
  intros H3 H4 H5 H6;
  split; right; repeat split.
   (* TP: boundedness *)
   stepl ((d+c)*(1/3)); [|qZ_numerals; field; auto]; apply Qlt_mult_neg_pos; auto.
   stepl ((d-c)*(1/3)); [|qZ_numerals; field; auto]; apply Qlt_mult_neg_pos; auto.
   (* TP: refining *)
   stepl ((a+b)+(1/3)*(d+c)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepr Zero; trivial;
    stepl (((a+b)*(Qone+Qone+Qone))-((c + d)*-Qone)); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero_neg; trivial.
   stepl ((a-b)+(1/3)*(d-c)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepr Zero; trivial;
    stepl ((d-c)-((b-a)*(Qone+Qone+Qone))); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero_neg; trivial.
   stepl ((1/3)*(c+d)-(a+b)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepr Zero; trivial;
    stepl ((c+d)-((a+b)*(Qone+Qone+Qone))); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero_neg; trivial.
   stepl ((1/3)*(d-c)+(b-a)); [|qZ_numerals; field; auto];
    apply Qmult_resp_Qle_pos_r with (Qone+Qone+Qone); auto; stepr Zero; trivial;
    stepl (((b-a)*(Qone+Qone+Qone))-((d-c)*-Qone)); [|qZ_numerals; field; auto];apply Qle_Qminus_Zero_neg; trivial.
Qed.

Lemma Incl_M_absorbs_Is_refining_M:forall mu d, Incl_M mu d -> Is_refining_M (product (inv_digit d) mu).
Proof.
 intros mu [ | | ] H_Incl; unfold inv_digit;
 [ apply Incl_M_absorbs_Is_refining_M_L
 | apply Incl_M_absorbs_Is_refining_M_R
 | apply Incl_M_absorbs_Is_refining_M_M
 ]; assumption.
Qed.

Close Scope Q_scope.
Close Scope Z_scope.

Lemma denom_nonvanishing_M_product_inv_L: forall mu r, denom_nonvanishing_M mu r -> ((-1) <= as_Moebius mu r <= 0)%R ->
                                                       denom_nonvanishing_M (product inv_LL mu) r.
Proof.
 intros ((a,b),(c,d)) r.
 unfold product, as_Moebius, denom_nonvanishing_M, inv_LL, fst, snd.
 intros H_denom [Hr1 Hr2].
 destruct (Q_eq_dec c a) as [Hca|Hca].
  (* a=c *)
  destruct (Q_eq_dec d b) as [Hbd|Hbd].
   (* d=b *)
   apply False_ind; subst c d; apply (Rlt_irrefl 1); apply Rle_lt_trans with 0; trivial;
   stepl ((a * r + b)/(a * r + b)); trivial; field; trivial.
   (* d<>b *)
   subst c; stepl ((1%Z/2%Z)*(d-b)).
    apply Rmult_resp_nonzero; realify_Q; auto; apply Rminus_eq_contra; apply Q_to_R_not_eq; trivial.
    realify_Q_goal; auto; change (IZR (Zneg xH)) with (- IZR (Zpos xH))%R; field; auto.
  (* c<>a *)
  realify_Q; auto.
  intros Hca.
  generalize (Rminus_eq_contra _ _ Hca); clear Hca; intro Hca.
  change (IZR 1) with 1; change (IZR (-1)) with (-1); change (IZR 2) with 2.
  stepl (1/2*((c-a)*r+(d-b))); [|field; auto].
  apply Rmult_resp_nonzero; auto.
  apply Rlinear_non_zero_3; trivial.
  intro H_false.
  apply (Rlt_irrefl 1).
  apply Rle_lt_trans with 0; trivial.
  stepl ((a * r + b)/(c * r + d)); trivial.
  subst r.
  repeat rewrite Rdiv_Rmult_numerator; trivial.
  repeat rewrite Rdiv_Rplus_Rmult; trivial.
  assert (H_linear:c * - (d - b) + (c - a) * d<>0);[rewrite (Rmult_comm (c-a) d); apply Rlinear_non_zero_2; trivial|].
  rewrite Rdiv_Rdiv_simplify; trivial.
  apply Req_Rdiv_Rone; trivial; ring.
Qed.

Lemma denom_nonvanishing_M_product_inv_R: forall mu r, denom_nonvanishing_M mu r -> (0 <= as_Moebius mu r <= 1)%R ->
                                                       denom_nonvanishing_M (product inv_RR mu) r.
Proof.
 intros ((a,b),(c,d)) r.
 unfold product, as_Moebius, denom_nonvanishing_M, inv_RR, fst, snd.
 intros H_denom [Hr1 Hr2].
 destruct (Q_eq_dec c (Qopp a)) as [Hca|Hca].
  (* -a=c *)
  destruct (Q_eq_dec d (Qopp b)) as [Hbd|Hbd].
   (* -d=b *)
   apply False_ind; subst c d; apply (Rlt_irrefl (-1)); apply Rlt_le_trans with 0; trivial;
    stepr ((a * r + b)/((-a)%Q * r + (-b)%Q)); trivial; generalize H_denom; clear H_denom; realify_Q;
    intros H_denom; field; trivial.
   (* -d<>b *)
   subst c; stepl ((1%Z/2%Z)*(d+b)).
    apply Rmult_resp_nonzero; realify_Q; auto; intros Hbd; stepl (d-(-b)); [|ring]; apply Rminus_eq_contra; trivial.
    realify_Q; auto; intros Hdb; field; auto.
  (* -c<>a *)
  realify_Q; auto; intros Hca.
  assert (Hca': c + a <> 0); [stepl (c-(-a));[|ring]; apply Rminus_eq_contra; assumption|].
  change (IZR 1) with 1; change (IZR (-1)) with (-1); change (IZR 2) with 2.
  stepl (1/2*((c+a)*r+(d+b))); [|field; auto].
  apply Rmult_resp_nonzero; auto.
  apply Rlinear_non_zero_3; trivial.
  intro H_false.
  apply (Rlt_irrefl (-1)).
  apply Rlt_le_trans with 0; trivial.
  stepr ((a * r + b)/(c * r + d)); trivial.
  subst r.
  repeat rewrite Rdiv_Rmult_numerator; trivial.
  repeat rewrite Rdiv_Rplus_Rmult; trivial.
  assert (H_linear:c * - (d + b) + (c + a) * d<>0);[rewrite (Rmult_comm (c+a) d); apply Rlinear_non_zero_2; trivial|].
  rewrite Rdiv_Rdiv_simplify; trivial.
  apply Req_Ropp_Rdiv_minus_Rone; trivial; ring.
Qed.

Lemma denom_nonvanishing_M_product_inv_M: forall mu r, denom_nonvanishing_M mu r -> ((-1/3) <= as_Moebius mu r <= 1/3)%R ->
                                                       denom_nonvanishing_M (product inv_MM mu) r.
Proof.
 intros ((a,b),(c,d)) r;
 unfold product, as_Moebius, denom_nonvanishing_M, inv_MM, fst, snd;
 intros H_denom [Hr1 Hr2];
 realify_Q; auto; simpl;
 stepl (1/3*(c*r+d)); auto;
 field; repeat apply Rmult_resp_nonzero; auto; stepl 3; auto; ring.
Qed.

Lemma Is_refining_M_L:Is_refining_M (LL:Matrix).
Proof.
 unfold map_digits, Is_refining_M, Bounded_M, fst, snd; split; left; repeat split; auto.
Qed.

Lemma Is_refining_M_R:Is_refining_M (RR:Matrix).
Proof.
 unfold map_digits, Is_refining_M, Bounded_M, fst, snd; split; left; repeat split; auto.
Qed.

Lemma Is_refining_M_M:Is_refining_M (MM:Matrix).
Proof.
 unfold map_digits, Is_refining_M, Bounded_M, fst, snd; split; left; repeat split; auto.
Qed.

Lemma as_Moebius_L: forall mu r, denom_nonvanishing_M mu r -> ((-1)<=as_Moebius mu r<=1)%R->
                                 as_Moebius (product LL mu) r= (((as_Moebius mu r)-1)/((as_Moebius mu r)+3))%R.
Proof.
 intros ((a,b),(c,d)) r H_denom H_r.
 generalize (denom_nonvanishing_M_product _ _ _ Is_refining_M_L H_r H_denom); generalize H_denom H_r; clear H_denom H_r.
 unfold product, as_Moebius, denom_nonvanishing_M, map_digits, fst, snd.
 realify_Q ;auto; simpl.
 replace (2+1) with 3; [|ring].
 intros H_denom Hr H_denom_product.
 assert (H_denom': a*r+b+(c*r+d)*3<>0);
  [apply Rmult_reg_nonzero_l with (1/2); stepl ((1/2*a+3/2*c)*r+(1/2*b+3/2*d)); trivial; field; auto|].
 rewrite Rdiv_Rminus_Rmult; trivial.
 rewrite Rdiv_Rplus_Rmult; trivial.
 rewrite Rdiv_Rdiv_simplify; trivial.
 transitivity (((a * r + b - (c * r + d) * 1)*(1/2))/((a * r + b + (c * r + d) * 3)*(1/2))).
  apply (f_equal2 Rdiv); field; auto.
  apply Rdiv_Rmult_simplify; auto.
Qed.

Lemma as_Moebius_R: forall mu r, denom_nonvanishing_M mu r -> ((-1)<=as_Moebius mu r<=1)%R->
                                 as_Moebius (product RR mu) r= (((as_Moebius mu r)+1)/((-as_Moebius mu r)+3))%R.
Proof.
 intros ((a,b),(c,d)) r H_denom H_r.
 generalize (denom_nonvanishing_M_product _ _ _ Is_refining_M_R H_r H_denom); generalize H_denom H_r; clear H_denom H_r.
 unfold product, as_Moebius, denom_nonvanishing_M, map_digits, fst, snd.
 realify_Q ;auto; simpl.
 replace (2+1) with 3; [|ring].
 intros H_denom Hr H_denom_product.
 assert (H_denom': -(a*r+b)+(c*r+d)*3<>0);
  [apply Rmult_reg_nonzero_l with (1/2); stepl ((-1/2*a+3/2*c)*r+(-1/2*b+3/2*d)); trivial; field; auto|].
 rewrite <- Rdiv_Ropp_numerator; trivial.
 repeat rewrite Rdiv_Rplus_Rmult; trivial.
 rewrite Rdiv_Rdiv_simplify; trivial.
 transitivity (((a * r + b + (c * r + d) * 1)*(1/2))/((-(a * r + b) + (c * r + d) * 3)*(1/2))).
  apply (f_equal2 Rdiv); field; auto.
  apply Rdiv_Rmult_simplify; auto.
Qed.

Lemma as_Moebius_M: forall mu r, denom_nonvanishing_M mu r -> ((-1)<=as_Moebius mu r<=1)%R->
                                 as_Moebius (product MM mu) r= ((as_Moebius mu r)/3)%R.
Proof.
 intros ((a,b),(c,d)) r H_denom H_r.
 generalize (denom_nonvanishing_M_product _ _ _ Is_refining_M_M H_r H_denom); generalize H_denom H_r; clear H_denom H_r.
 unfold product, as_Moebius, denom_nonvanishing_M, map_digits, fst, snd.
 realify_Q ;auto; simpl.
 replace (2+1) with 3; [|ring].
 intros H_denom Hr H_denom_product.
 rewrite Rdiv_Rdiv_Rmult_numerator; auto.
 apply (f_equal2 Rdiv); field; auto.
Qed.

Lemma as_Moebius_inv_L: forall mu r, denom_nonvanishing_M mu r -> ((-1) <= as_Moebius mu r <= 0)%R -> as_Moebius (product inv_LL mu) r = ((3 * as_Moebius mu r + 1) / (- as_Moebius mu r + 1))%R.
Proof.
 intros ((a,b),(c,d)) r H_denom H_r;
 generalize (denom_nonvanishing_M_product_inv_L _ _ H_denom H_r); generalize H_denom H_r; clear H_denom H_r;
 unfold product, as_Moebius, denom_nonvanishing_M, inv_LL, fst, snd;
 realify_Q ;auto; simpl;
 replace (2+1) with 3; [|ring];
 intros H_denom Hr H_denom_product.
 assert (H_denom':    - (a * r + b) + (c * r + d) * 1 <> 0);
 [apply Rmult_reg_nonzero_l with (1/2); stepl ((-1/2*a+1/2*c)*r+(-1/2*b+1/2*d)); trivial; field; auto|].
 rewrite <- Rdiv_Ropp_numerator; trivial.
 rewrite Rdiv_Rmult_numerator; trivial.
 repeat rewrite Rdiv_Rplus_Rmult; trivial.
 rewrite Rdiv_Rdiv_simplify; trivial.
 transitivity (((3 * (a * r + b) + (c * r + d) * 1)*(1/2))/((- (a * r + b) + (c * r + d) * 1)*(1/2))).
  apply (f_equal2 Rdiv); field; auto.
  apply Rdiv_Rmult_simplify; auto.
Qed.

Lemma as_Moebius_inv_R: forall mu r, denom_nonvanishing_M mu r -> (0 <= as_Moebius mu r <= 1)%R -> as_Moebius (product inv_RR mu) r = ((3 * as_Moebius mu r - 1) / (as_Moebius mu r + 1))%R.
Proof.
 intros ((a,b),(c,d)) r H_denom H_r;
 generalize (denom_nonvanishing_M_product_inv_R _ _ H_denom H_r); generalize H_denom H_r; clear H_denom H_r;
 unfold product, as_Moebius, denom_nonvanishing_M, inv_RR, fst, snd;
 realify_Q ;auto; simpl;
 replace (2+1) with 3; [|ring];
 intros H_denom Hr H_denom_product.
 assert (H_denom':    (a * r + b) + (c * r + d) * 1 <> 0);
 [apply Rmult_reg_nonzero_l with (1/2); stepl ((1/2*a+1/2*c)*r+(1/2*b+1/2*d)); trivial; field; auto|].
 rewrite Rdiv_Rmult_numerator; trivial.
 repeat rewrite Rdiv_Rplus_Rmult; trivial.
 rewrite Rdiv_Rminus_Rmult; trivial.
 rewrite Rdiv_Rdiv_simplify; trivial.
 transitivity (((3 * (a * r + b) - (c * r + d) * 1)*(1/2))/(((a * r + b) + (c * r + d) * 1)*(1/2))).
  apply (f_equal2 Rdiv); field; auto.
  apply Rdiv_Rmult_simplify; auto.
Qed.

Lemma as_Moebius_inv_M: forall mu r, denom_nonvanishing_M mu r -> ((-1/3) <= as_Moebius mu r <= 1/3)%R -> as_Moebius (product inv_MM mu) r = (3 * as_Moebius mu r)%R.
Proof.
 intros ((a,b),(c,d)) r H_denom H_r;
 generalize (denom_nonvanishing_M_product_inv_M _ _ H_denom H_r); generalize H_denom H_r; clear H_denom H_r;
 unfold product, as_Moebius, denom_nonvanishing_M, inv_MM, fst, snd;
 realify_Q ;auto; simpl;
 replace (2+1) with 3; [|ring];
 intros H_denom Hr H_denom_product.
 rewrite Rdiv_Rmult_numerator; trivial.
 rewrite (Rmult_comm 3); rewrite <- Rdiv_Rdiv_Rmult_denominator; auto.
 apply (f_equal2 Rdiv); field; auto.
Qed.

Lemma as_Moebius_L_min_one: as_Moebius LL (-1)%R = (-1)%R.
Proof.
 unfold as_Moebius, map_digits, fst, snd;
 realify_Q; auto; simpl; replace (2+1) with 3; [|ring]; replace (1 / 2 * -1 + 3 / 2) with 1; field; auto.
Qed.

Lemma as_Moebius_R_min_one: as_Moebius RR (-1)%R = (0)%R.
Proof.
 unfold as_Moebius, map_digits, fst, snd;
 realify_Q; auto; simpl; replace (2+1) with 3; [|ring]; replace (-1 / 2 * -1 + 3 / 2) with 2; field; auto.
Qed.

Lemma as_Moebius_M_min_one: as_Moebius MM (-1)%R = (-1/3)%R.
Proof.
 unfold as_Moebius, map_digits, fst, snd;
 realify_Q; auto; simpl; replace (2+1) with 3; [|ring]; replace (0/1 * -1 + 3 / 1) with 3; field; auto.
Qed.

Lemma as_Moebius_L_one: as_Moebius LL (1)%R = (0)%R.
Proof.
 unfold as_Moebius, map_digits, fst, snd;
 realify_Q; auto; simpl; replace (2+1) with 3; [|ring]; replace (1 / 2 * 1 + 3 / 2) with 2; field; auto.
Qed.

Lemma as_Moebius_R_one: as_Moebius RR (1)%R = (1)%R.
Proof.
 unfold as_Moebius, map_digits, fst, snd;
 realify_Q; auto; simpl; replace (2+1) with 3; [|ring]; replace (-1 / 2 * 1 + 3 / 2) with 1; field; auto.
Qed.

Lemma as_Moebius_M_one: as_Moebius MM (1)%R = (1/3)%R.
Proof.
 unfold as_Moebius, map_digits, fst, snd;
 realify_Q; auto; simpl; replace (2+1) with 3; [|ring]; replace (0/1 * 1 + 3 / 1) with 3; field; auto.
Qed.

Lemma as_Moebius_endpoints_digits_base_l: forall d:Digit, ((-1)%R<=as_Moebius d (-1)%R)%R.
Proof.
 intros [ | | ];
 [ rewrite as_Moebius_L_min_one
 | rewrite as_Moebius_R_min_one
 | rewrite as_Moebius_M_min_one
 ]; fourier.
Qed.

Lemma as_Moebius_endpoints_digits_base_r: forall d:Digit, (as_Moebius d 1%R<=1)%R.
Proof.
 intros [ | | ];
 [ rewrite as_Moebius_L_one
 | rewrite as_Moebius_R_one
 | rewrite as_Moebius_M_one
 ]; fourier.
Qed.

Open Scope Z_scope.
Open Scope Q_scope.

Lemma Is_refining_M_property_fold:  forall mu, (forall r, (-1 <= r <= 1)%R -> denom_nonvanishing_M mu r)->
  (forall (r : R),(-1 <= r <= 1)%R -> (-1<=as_Moebius mu r<=1)%R) ->  Is_refining_M mu.
Proof.
 intros ((a,b),(c,d)) H_denom H_mu.
 unfold denom_nonvanishing_M in H_denom; simpl in H_denom.
 unfold as_Moebius in H_mu; simpl in H_mu.
 assert (H_dcmin:d-c <>Zero);
 [apply Q_to_R_Qneq; realify_Q; stepl (c*(-1)+d)%R; [|ring]; apply (H_denom (-1)%R min_one_is_in_base_interval)|].
 assert (H_dc:d+c <>Zero);
 [apply Q_to_R_Qneq; realify_Q_goal; stepl (c*1+d)%R; [|ring]; apply (H_denom (1)%R one_is_in_base_interval)|].
 assert (Hd: d<>Zero);
 [apply Q_to_R_Qneq; realify_Q_goal; stepl (c*0+d)%R; [|ring]; apply H_denom; split; fourier|].
 split; [ apply denom_nonvanishing_M_Bounded_M; assumption|].
 simpl.
 destruct (denom_nonvanishing_M_Bounded_M ((a,b),(c,d)) H_denom) as [[H1 H2]|[H1 H2]]; simpl in H1, H2.
  (* Zero < d + c /\ Zero < d - c *)
  left; repeat split.
   (* 1 *)
   apply Qle_Qdiv_denom_pos_nonneg with (d+c); trivial;
   assert (Hcd:c*Qone+d <> Zero);[stepl (d+c); trivial; ring|];
   stepr (((a*Qone + b)/(c*Qone + d)) - (-  Qone)); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj1 (H_mu (1)%R one_is_in_base_interval)).
   (* 2 *)
   apply Qle_Qdiv_denom_pos_nonneg with (d-c); trivial;
   assert (Hcd:c*(-Qone)+d <> Zero);[stepl (d-c); trivial; ring|];
   stepr (Qone -((a*(-Qone) + b)/(c*(-Qone) + d))); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj2 (H_mu (-1)%R min_one_is_in_base_interval)).
   (* 3 *)
   apply Qle_Qdiv_denom_pos_nonneg with (d+c); trivial;
   assert (Hcd:c*Qone+d <> Zero);[stepl (d+c); trivial; ring|];
   stepr (Qone-((a*Qone + b)/(c*Qone + d))); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj2 (H_mu (1)%R one_is_in_base_interval)).
   (* 4 *)
   apply Qle_Qdiv_denom_pos_nonneg with (d-c); trivial;
   assert (Hcd:c*(-Qone)+d <> Zero);[stepl (d-c); trivial; ring|];
   stepr (((a*(-Qone) + b)/(c*(-Qone) + d))-(-Qone)); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj1 (H_mu (-1)%R min_one_is_in_base_interval)).
  (* d + c < Zero /\ d - c , Zero *)
  right; repeat split.
   (* 1 *)
   apply Qle_Qdiv_denom_neg_nonneg with (d+c); trivial;
   assert (Hcd:c*Qone+d <> Zero);[stepl (d+c); trivial; ring|];
   stepr (((a*Qone + b)/(c*Qone + d)) - (-  Qone)); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj1 (H_mu (1)%R one_is_in_base_interval)).
   (* 2 *)
   apply Qle_Qdiv_denom_neg_nonneg with (d-c); trivial;
   assert (Hcd:c*(-Qone)+d <> Zero);[stepl (d-c); trivial; ring|];
   stepr (Qone -((a*(-Qone) + b)/(c*(-Qone) + d))); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj2 (H_mu (-1)%R min_one_is_in_base_interval)).
   (* 3 *)
   apply Qle_Qdiv_denom_neg_nonneg with (d+c); trivial;
   assert (Hcd:c*Qone+d <> Zero);[stepl (d+c); trivial; ring|];
   stepr (Qone-((a*Qone + b)/(c*Qone + d))); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj2 (H_mu (1)%R one_is_in_base_interval)).
   (* 4 *)
   apply Qle_Qdiv_denom_neg_nonneg with (d-c); trivial;
   assert (Hcd:c*(-Qone)+d <> Zero);[stepl (d-c); trivial; ring|];
   stepr (((a*(-Qone) + b)/(c*(-Qone) + d))-(-Qone)); [|field; auto];
   apply Qle_Qminus_Zero;
   apply Q_to_R_Qle;
   realify_Q_goal; trivial;
   apply (proj1 (H_mu (-1)%R min_one_is_in_base_interval)).
Qed.

Lemma Is_refining_M_product: forall mu mu', Is_refining_M mu -> Is_refining_M mu' -> Is_refining_M (product mu mu').
Proof.
 intros ((a,b),(c,d)) ((a',b'),(c',d')) H_mu H_mu'.
 apply Is_refining_M_property_fold; intros r Ht.
  apply denom_nonvanishing_M_product; trivial;
  [ apply Is_refining_M_property
  | apply Is_refining_M_denom_nonvanishing_M
  ]; trivial.

  rewrite as_Moebius_product; trivial;
  repeat apply Is_refining_M_property; trivial;
  apply Is_refining_M_denom_nonvanishing_M; trivial.
Qed.

Lemma Is_refining_M_product_hd:forall mu alpha, Is_refining_M mu -> Is_refining_M (product mu (map_digits (hd alpha))).
Proof.
 intros mu [[ | | ] alpha] H_mu;
 unfold hd; apply Is_refining_M_product; trivial;
 [apply Is_refining_M_L|apply Is_refining_M_R|apply Is_refining_M_M].
Qed.

Lemma det_nonneg_refining_endpoints:forall a b c d, Is_refining_M ((a,b),(c,d)) -> Zero <= a*d-b*c  ->
       as_Moebius_Q ((a,b),(c,d)) (-Qone)<=as_Moebius_Q ((a,b),(c,d)) Qone.
Proof.
 intros a b c d [H_bounded _] H_det.
 unfold as_Moebius_Q; simpl.
 assert (H : (-1 <= 1)%R) by fourier.
 generalize (det_nonneg_nondecreasing _ _ _ _ (-1) 1 H_bounded H_det min_one_is_in_base_interval one_is_in_base_interval H).
  let t_local:=  (destruct H_bounded as [[H1 H2]|[H1 H2]]; simpl in H1, H2; auto) in rationalify_R_goal;
 [ intro H_Rle; apply Q_to_R_Qle; assumption
 | stepl (d-c); try ring; t_local
 | stepl (d+c); try ring;t_local
 ].
Qed.

Lemma det_nonpos_refining_endpoints:forall a b c d, Is_refining_M ((a,b),(c,d)) -> a*d-b*c<= Zero  ->
       as_Moebius_Q ((a,b),(c,d)) Qone<=as_Moebius_Q ((a,b),(c,d)) (-Qone).
Proof.
 intros a b c d [H_bounded _] H_det.
 unfold as_Moebius_Q; simpl.
 assert (H : (-1 <= 1)%R) by fourier.
 generalize (det_nonpos_nonincreasing _ _ _ _ (-1) 1 H_bounded H_det min_one_is_in_base_interval one_is_in_base_interval H).
 let t_local:=  (destruct H_bounded as [[H1 H2]|[H1 H2]]; simpl in H1, H2; auto) in rationalify_R_goal;
 [ intro H_Rle; apply Q_to_R_Qle; assumption
 | stepl (d-c); try ring; t_local
 | stepl (d+c); try ring;t_local
 ].
Qed.

Lemma as_Moebius_Q_product: forall mu mu' (q:Q),
  denom_nonvanishing_M mu' q -> -1 <= as_Moebius_Q mu' q -> as_Moebius_Q mu' q <= 1 -> Is_refining_M mu ->
    as_Moebius_Q (product mu mu') q = as_Moebius_Q mu (as_Moebius_Q mu' q).
Proof.
 intros mu mu' q H_denom.
 qZ_numerals.
 intros H_l H_r H_refining.
 apply Q_to_R_Qeq.
 realify_Q.
 intros H_l H_r.
 assert (H1: denom_nonvanishing_M mu (as_Moebius mu' q));
  [apply Is_refining_M_denom_nonvanishing_M; trivial; revert H_l H_r; rewrite Q_to_R_as_Moebius; auto|].
 assert (H2: denom_nonvanishing_M (product mu mu') q);
  [apply denom_nonvanishing_M_product; trivial;revert H_l H_r; rewrite Q_to_R_as_Moebius; auto|].
 revert H_l H_r; repeat rewrite Q_to_R_as_Moebius; trivial.
 intros H_l H_r.
 apply as_Moebius_product; auto.
Qed.

Lemma Is_refining_M_product_init_pure:forall alpha n, Is_refining_M (product_init_pure alpha n).
Proof.
 intros alpha n; revert alpha; induction n; intros [d alpha]; unfold product_init_pure.
  (* O *)
  simpl; apply Is_refining_M_idM.
  (* S n *)
  rewrite Streams_addenda.take_S_n;
  rewrite (Streams_addenda.fold_right_cons product).
  apply Is_refining_M_product.
   simpl; destruct d; [apply Is_refining_M_L|apply Is_refining_M_R|apply Is_refining_M_M].
   unfold map_reals; rewrite Streams_addenda.map_spelled; simpl; apply (IHn alpha).
Qed.

Close Scope Q_scope.
Close Scope Z_scope.

author: verimath

repo: real

save_path: github-repos/coq/verimath-real


text:
From Hammer Require Import Hammer Reconstr.

Require Import Arith.

(* disable the preliminary crush-like tactic *)
Set Hammer CrushLimit 0.

Set Hammer Eprover.
Set Hammer Vampire.
Set Hammer Z3.

Lemma good1: forall n, n + n = 2 * n.
Proof.
  hammer.
(*
Extracting features...
Running provers (using 8 threads)...
Z3 succeeded
 - dependencies: Coq.Arith.PeanoNat.Nat.double_twice
 - definitions: Coq.Init.Nat.double
Trying tactic Reconstr.reasy...
Tactic Reconstr.reasy succeeded.
Replace the hammer tactic with:
        Reconstr.reasy (@Coq.Arith.PeanoNat.Nat.double_twice) (@Coq.Init.Nat.double).
*)
  Restart.
  Reconstr.reasy (@Coq.Arith.PeanoNat.Nat.double_twice) (@Coq.Init.Nat.double).
Qed.

Lemma bad1: forall n, n + n + n = 2 * n.
Proof.
  hammer.

(*
Error: Hammer failed: ATPs failed to find a proof
 *)
Abort.

author: samuelgruetter

repo: counterexamples


text:
Require Import Nat.
Require Import PPS.Env.

Set Implicit Arguments.

Section ImpModell.

  Inductive Val : Type :=
  | VInt      : nat -> Val
  | Undefined : Val.

  Definition Ref := nat.
  Definition Name := nat.

  Inductive Ty : Type :=
  | Int    : Ty.

  Definition Env := listMap Name Ref.

  Definition Memory := totalMap Ref Val.

  Reserved Notation "E '|-l' val" (at level 80).

  Notation "E ; x" := (union (cons x nil) E) (at level 40).

  Section EnvRelation.

    Inductive Lookup : Env -> (Name * Ref) -> Prop :=
    | LAxiom : forall E x beta,
        E; (x, beta) |-l (x, beta)
    | LRule : forall E x y beta gamma,
        E |-l (x,beta) ->
        E; (y, gamma) |-l (x, beta)
    where "E '|-l' val" := (Lookup E val).

  End EnvRelation.
  Notation "E '|-l' val" := (Lookup E val).

  Section Exp.

    Inductive ArithExp : Type :=
    | Num : nat -> ArithExp
    | Var : Name -> ArithExp
    | Plus : ArithExp -> ArithExp -> ArithExp
    | Mult : ArithExp -> ArithExp -> ArithExp.

  End Exp.

  Reserved Notation "EM '|-R' e ⇓ v" (at level 80).
  Section EvalR.

    Inductive evalR : (Env * Memory) -> ArithExp -> nat -> Prop :=
    | EvVarR  : forall E M x l val,
        E |-l (x, l) ->
        M l = VInt val ->
        (E,M) |-R (Var x) ⇓ val
    | EvNumR  : forall E M n, (E,M) |-R Num n ⇓ n

    | EvPlusR : forall E M e1 e2 v1 v2 v,
        (E,M) |-R e1 ⇓ v1 ->
        (E,M) |-R e2 ⇓ v2 ->
        v = v1 + v2 ->
        (E,M) |-R Plus e1 e2 ⇓ v

    | EvMultR : forall E M e1 e2 v1 v2 v,
        (E,M) |-R e1 ⇓ v1 ->
        (E,M) |-R e2 ⇓ v2 ->
        v = v1 * v2 ->
        (E,M) |-R Mult e1 e2 ⇓ v

    where "EM '|-R' e ⇓ v" := (evalR EM e v).

  End EvalR.
  Notation "EM '|-R' e ⇓ v" := (evalR EM e v).

  Section Stm.

    Inductive Stm : Type :=
    | Ass    : Name -> ArithExp -> Stm
    | Decl   : Name -> Stm
    | Seq    : Stm -> Stm -> Stm.

    Reserved Notation "'<' E '|' M '>' alpha '<' E' '|' M' '>'"
             (at level 40, E at level 39, M at level 39, M' at level 40,
              E' at level 39, alpha at level 39).
    Notation "M '[' var ↦ val ']'" := (updateTMap eqb M var val)
                                    (at level 40, right associativity).

    Inductive IsFree : Ref -> Memory -> Prop :=
    | isFree : forall x M,
        M x = Undefined ->
        IsFree x M.

    Inductive eval : Env * Memory -> Stm -> Env * Memory -> Prop :=
    | EvDecl : forall E M l n,
        IsFree l M ->
        < E | M > Decl n < E; (n, l) | M [l ↦ VInt 0] >

    | EvAss   : forall E M n e l v,
        IsFree l M ->
        (E,M) |-R e ⇓ v ->
        < E | M > Ass n e < E | M [n ↦ VInt v] >

    | EvSeq   : forall E M E' M' E'' M'' S1 S2,
        < E | M > S1 < E' | M' > ->
        < E' | M' > S2 < E'' | M'' > ->
        < E | M > Seq S1 S2 < E'' | M'' >
    where "'<' E '|' M '>' alpha '<' E' '|' M' '>'" := (eval (E,M) alpha (E',M')).

  End Stm.

End ImpModell.

author: nbun

repo: pps-coq

save_path: github-repos/coq/nbun-pps-coq


text:
From Coq Require Import String List Program.Wf Arith.
From minilog Require Import utils set.
Import ListNotations.

(** * An simple algebra of first order terms to describe databases *)

(** ** Syntax *)

(** A [datum] is function symbol applied to a set of [datum] *)
Inductive datum :=
  | Dapp : string -> list datum -> datum.

Definition Dcst (c : string) := Dapp c [].

(** Patterns over [datum] *)
Inductive pattern :=
  | Pvar : nat -> pattern
  | Papp : string -> list pattern -> pattern.

Definition Pcst (c : string) := Papp c [].

(** ** Induction Principles *)

(** We need a stronger induction principle that
    the one infered by Coq.
*)
Fixpoint datum_induction
  (P : datum -> Prop)
  (Hind : forall f dats, Forall P dats -> P (Dapp f dats))
  (dat : datum) : P dat :=
  match dat with
  | Dapp f dats =>
    let helper := fix F l :=
      match l as l' return Forall P l' with
      | nil => Forall_nil P
      | cons x xs =>
        @Forall_cons _ P x xs (datum_induction P Hind x) (F xs)
      end
    in
    Hind f dats (helper dats)
  end.

(** We need a stronger induction principle that
  the one infered by Coq.
*)
Fixpoint pattern_induction
  (P : pattern -> Prop)
  (H1 : forall x, P (Pvar x))
  (Hind : forall f dats, Forall P dats -> P (Papp f dats))
  (pat : pattern) : P pat :=
  match pat with
  | Pvar x => H1 x
  | Papp f pats =>
    let helper := fix F l :=
      match l as l' return Forall P l' with
      | nil => Forall_nil P
      | cons x xs =>
        @Forall_cons _ P x xs (pattern_induction P H1 Hind x) (F xs)
      end
    in
    Hind f pats (helper pats)
  end.

(** Free variables of a pattern *)
Inductive free_var : pattern -> pset nat :=
| free_var_var x : free_var (Pvar x) x
| free_var_app x f l :
  Exists (fun p => free_var p x) l -> free_var (Papp f l) x.

Fixpoint is_free (p : pattern) (x : nat) : bool :=
  match p with
  | Pvar y => (y =? x)%nat
  | Papp _ pats => existsb (fun p => is_free p x) pats
end.

Lemma is_free_iff:
  forall p x, is_free p x = true <-> free_var p x.
Proof.
  induction p using pattern_induction.
  - intros y. split.
    + intros Hy. simpl in Hy.
      apply Nat.eqb_eq in Hy as ->.
      econstructor.
    + inversion 1; subst. simpl.
      apply Nat.eqb_refl.
  - intros x. split.
    + simpl.
      intros [pat [Hpat1 Hpat2]]%existsb_exists.
      rewrite Forall_forall in H.
      pose proof (proj1 (H _ Hpat1 _) Hpat2).
      econstructor. apply Exists_exists.
      now exists pat.
    + inversion_clear 1 as [|? ? ? [pat [Hpat1 Hpat2]]%Exists_exists]; subst.
      rewrite Forall_forall in H.
      pose proof (proj2 (H _ Hpat1 _) Hpat2).
      apply existsb_exists.
      now exists pat.
Qed.

Theorem is_free_false:
  forall p x, is_free p x = false -> ~free_var p x.
Proof.
  intros p x H Hcontr.
  apply is_free_iff in Hcontr.
  now rewrite Hcontr in H.
Qed.

Declare Scope datum.
Delimit Scope datum with datum.
Declare Scope pattern.
Delimit Scope pattern with pattern.

Fixpoint datum_size (d : datum) : nat :=
  match d with
  | Dapp _ dats =>
    1 + sum datum_size dats
  end.

Fixpoint pattern_size (p : pattern) : nat :=
  match p with
  | Pvar _ => 1
  | Papp _ pats =>
    1 + sum pattern_size pats
  end.

(** Boolean equality for [datum] *)
Fixpoint datum_eqb (d1 d2 : datum) : bool :=
  match d1, d2 with
  | Dapp f1 l1, Dapp f2 l2 =>
    let fix F l1 l2 :=
      match l1, l2 with
      | [], [] => true
      | x::xs, y::ys => andb (datum_eqb x y) (F xs ys)
      | _, _ => false
      end
    in andb (f1 =? f2)%string (F l1 l2)
  end.

Lemma datum_eqb_refl:
  forall x, datum_eqb x x = true.
Proof.
  refine (datum_induction _ _); intros; simpl.
  apply Bool.andb_true_iff; split.
  apply String.eqb_refl.
  induction dats; try easy.
  inversion H as [| ? ? H1 H2]; subst.
  rewrite H1. simpl.
  now apply IHdats.
Qed.

Lemma datum_eqb_eq:
  forall x y, datum_eqb x y = true <-> x = y.
Proof.
  refine (datum_induction _ _).
  intros f dats IH y. split; intros.
  + destruct y; simpl in H; try easy.
    apply Bool.andb_true_iff in H as [->%String.eqb_eq H].
    induction l in dats, H, IH, s |-*.
    destruct dats; try easy.
    destruct dats; try easy.
    apply Bool.andb_true_iff in H as [H1 H2].
    inversion_clear IH as [|? ? H3 H4].
    apply H3 in H1 as ->.
    assert (Dapp s dats = Dapp s l).
    apply IHl; auto.
    now inversion H.
  + destruct y; simpl in H; try easy.
    inversion_clear H.
    apply datum_eqb_refl.
Qed.

Lemma datum_eqb_neq:
  forall x y, datum_eqb x y = false -> x <> y.
Proof.
  intros * Hfalse Hcontr.
  apply datum_eqb_eq in Hcontr.
  now rewrite Hfalse in Hcontr.
Qed.

Infix "=?" := (datum_eqb)%datum.

author: acorrenson

repo: minilog

save_path: github-repos/coq/acorrenson-minilog


text:
(*
 * Vericert: Verified high-level synthesis.
 * Copyright (C) 2020 Yann Herklotz <yann@yannherklotz.com>
 *               2020 James Pollard <j@mes.dev>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 *)

(*Set Implicit Arguments.

Require Export compcert.lib.Maps.

Require Import compcert.common.Errors.

Require Import vericert.common.Vericertlib.

Import PTree.

Local Open Scope error_monad_scope.

(*|
Instance of traverse for ``PTree`` and ``Errors``. This should maybe be generalised in the future.
|*)

Module PTree.

Fixpoint xtraverse (A B : Type) (f : positive -> A -> res B) (m : PTree.t A) (i : positive)
         {struct m} : res (PTree.t B) :=
  match m with
  | Leaf => OK Leaf
  | Node l o r =>
    let newo :=
        match o with
        | None => OK None
        | Some x =>
          match f (prev i) x with
          | Error err => Error err
          | OK val => OK (Some val)
          end
        end in
    match newo with
    | OK no =>
      do nl <- xtraverse f l (xO i);
      do nr <- xtraverse f r (xI i);
      OK (Node nl no nr)
    | Error msg => Error msg
    end
  end.

Definition traverse (A B : Type) (f : positive -> A -> res B) m := xtraverse f m xH.

Definition traverse1 (A B : Type) (f : A -> res B) := traverse (fun _ => f).

End PTree.
*)

author: ymherklotz

repo: vericert

save_path: github-repos/coq/ymherklotz-vericert


text:
Require Import Setoid Program.Basics.

 Parameter P : nat -> Prop.
 Parameter R : nat -> nat -> Prop.

 Add Parametric Morphism : P
   with signature R ++> impl as PM1.
 Admitted.

 Add Parametric Morphism : P
   with signature R --> impl as PM2.
 Admitted.

 Goal forall x y, R x y -> P y -> P x.
 Proof.
   intros x y H1 H2.
   rewrite H1.
   auto.
 Qed.

repo: Coq-misc


text:
Require Import Coq.ZArith.ZArith.
Require Import Coq.micromega.Lia.
Require Import Crypto.Util.ZUtil.Definitions.
Require Import Crypto.Util.ZUtil.Notations.
Require Import Crypto.Util.ZUtil.Hints.Core.
Require Import Crypto.Util.ZUtil.Hints.Ztestbit.
Require Import Crypto.Util.ZUtil.Tactics.ZeroBounds.
Require Import Crypto.Util.ZUtil.Testbit.
Require Import Crypto.Util.Tactics.BreakMatch.
Local Open Scope Z_scope.

Module Z.
  Lemma pow2_mod_spec : forall a b, (0 <= b) -> Z.pow2_mod a b = a mod (2 ^ b).
  Proof.
    intros.
    unfold Z.pow2_mod.
    rewrite Z.land_ones; auto.
  Qed.
  Hint Rewrite <- Z.pow2_mod_spec using zutil_arith : convert_to_Ztestbit.

  Lemma pow2_mod_0_r : forall a, Z.pow2_mod a 0 = 0.
  Proof.
    intros; rewrite Z.pow2_mod_spec, Z.mod_1_r; reflexivity.
  Qed.

  Lemma pow2_mod_0_l : forall n, 0 <= n -> Z.pow2_mod 0 n = 0.
  Proof.
    intros; rewrite Z.pow2_mod_spec, Z.mod_0_l; try reflexivity; try apply Z.pow_nonzero; lia.
  Qed.

  Lemma pow2_mod_split : forall a n m, 0 <= n -> 0 <= m ->
                                       Z.pow2_mod a (n + m) = Z.lor (Z.pow2_mod a n) ((Z.pow2_mod (a >> n) m) << n).
  Proof.
    intros; cbv [Z.pow2_mod].
    apply Z.bits_inj'; intros.
    repeat progress (try break_match; autorewrite with Ztestbit zsimplify; try reflexivity).
    try match goal with H : ?a < ?b |- context[Z.testbit _ (?a - ?b)] =>
      rewrite !Z.testbit_neg_r with (n := a - b) by lia end.
    autorewrite with Ztestbit; reflexivity.
  Qed.

  Lemma pow2_mod_pow2_mod : forall a n m, 0 <= n -> 0 <= m ->
                                          Z.pow2_mod (Z.pow2_mod a n) m = Z.pow2_mod a (Z.min n m).
  Proof.
    intros; cbv [Z.pow2_mod].
    apply Z.bits_inj'; intros.
    apply Z.min_case_strong; intros; repeat progress (try break_match; autorewrite with Ztestbit zsimplify; try reflexivity).
  Qed.

  Lemma pow2_mod_pos_bound a b : 0 < b -> 0 <= Z.pow2_mod a b < 2^b.
  Proof.
    intros; rewrite Z.pow2_mod_spec by lia.
    auto with zarith.
  Qed.
  Hint Resolve pow2_mod_pos_bound : zarith.

  Lemma pow2_mod_id_iff : forall a n, 0 <= n ->
                                      (Z.pow2_mod a n = a <-> 0 <= a < 2 ^ n).
  Proof.
    intros a n H.
    rewrite Z.pow2_mod_spec by assumption.
    assert (0 < 2 ^ n) by Z.zero_bounds.
    rewrite Z.mod_small_iff by lia.
    split; intros; intuition lia.
  Qed.
End Z.

author: dip-proto

repo: fiat-crypto

save_path: github-repos/coq/dip-proto-fiat-crypto