Require Export Reals.
Require Export Field.
Require Export ZArith.
Require Export Omega.
Require Export ZArithRing.
Require Export Wellfounded.
Require Export Zwf.
Require Export Fourier.

Open Scope Z_scope.

CoInductive stream (A:Set) : Set :=
  Cons : A -> stream A -> stream A.

Implicit Arguments Cons.

Definition str_decompose (A:Set)(x:stream A) :=
  match x with Cons a tl => Cons a tl end.

Implicit Arguments str_decompose.

Theorem str_dec_thm : forall A:Set, forall x : stream A,
  x = str_decompose x.
intros; case x; auto.
Qed.

Implicit Arguments str_dec_thm.

Inductive idigit : Set :=  R | L | C.

CoInductive represents : stream idigit ->
                  Rdefinitions.R -> Prop :=
  reprL : forall s r, represents s r -> (0 <= r <= 1)%R ->
           represents (Cons L s) (r/2)
| reprR : forall s r, represents s r -> (0 <= r <= 1)%R ->
           represents (Cons R s) ((r+1)/2)
| reprC : forall s r, represents s r -> (0 <= r <= 1)%R ->
           represents (Cons C s) ((2*r+1)/4).


(* returns the interval represented by the n first digits of the stream,
   the value is given by (a, b, k) so that the interval actually is
   (a/(2**k), b/(2**k)) *)

Fixpoint bounds (n:nat) (s:stream idigit) {struct n} : Z*Z*Z :=
  match n with
    0%nat => (0,1,0)
  | S p =>
    let (d,tl) := s in
    let (p,e) := bounds p tl in
    let (a,b) := p in
    match d with
      L => (a,b,e+1)
    | R => (a+Zpower 2 e,b*Zpower 2 e,e+1)
    | C => (2*a+Zpower 2 e, 2*b+Zpower 2 e, e+2)
    end
  end.

Definition si_un (s:stream idigit) (n:nat) : Rdefinitions.R :=
   let (p,k) := bounds n s in
   let (a,b) := p in
    (IZR a/IZR(2^k))%R.

Definition lift (d:idigit) (x:Rdefinitions.R) : Rdefinitions.R :=
  match d with
    L => (x/2)%R
  | R => (x/2 + 1/2)%R
  | C => (x/2 + 1/4)%R
  end.

Lemma pow2_larger :
  forall n, (INR n <= 2^n)%R.
induction n.
simpl; fourier.
replace (S n) with (1+n)%nat.
rewrite pow_add; rewrite pow_1.
assert (Hlem: (forall a b, 1<= b -> a<=b -> 1+a<= 2*b)%R).
intros; fourier.
rewrite plus_INR; simpl (INR 1).
apply Hlem; auto.
apply pow_R1_Rle;fourier.
omega.
Qed.

Theorem INR_Zabs_nat_eq :
  forall z, 0 <= z -> (IZR z = INR(Zabs_nat z))%R.
unfold Zabs_nat.
intros z; case z.
simpl.
intros; auto.
unfold IZR.
intros; auto.
compute.
intros x H'; elim H'; auto.
Qed.

Lemma Rlt_neq_0 : forall r, (0 < r)%R -> (r <> 0)%R.
intros r Hpos; replace r with (0 + r)%R.
apply tech_Rplus.
fourier.
auto.
field.
Qed.

Lemma pow_half_smaller :
 forall r, (0 < r)%R -> (1/(2^(Zabs_nat(up (1/r)))) <= r)%R.
intros r Hpos.
elim (archimed (1/r)); intros Harc1 Harc2.
apply (Rmult_le_reg_l (2^(Zabs_nat (up (1/r))))).
apply pow_lt.
fourier.
replace (2 ^ Zabs_nat (up (1 / r)) * (1 / 2 ^ Zabs_nat (up (1 / r))))%R
 with 1%R.
apply Rle_trans with ((INR (Zabs_nat (up (1/r))))*r)%R.
rewrite <- INR_Zabs_nat_eq.
pattern 1%R at 1; replace 1%R with (1/r * r)%R.
apply Rmult_le_compat_r.
fourier.
fourier.
field.
apply Rlt_neq_0; assumption.
apply le_IZR.
apply Rle_trans with (1/r)%R.
replace (1/r)%R with (/r)%R.
simpl IZR.
assert (0 < / r)%R.
apply Rinv_0_lt_compat; auto.
fourier.
field.
apply Rlt_neq_0; auto.
fourier.
apply Rmult_le_compat_r.
fourier.
apply pow2_larger.
field.
apply pow_nonzero.
apply (not_O_IZR 2);omega.
Qed.

Lemma Zmult_non0 :
  forall x y, x <> 0-> y<> 0 -> x*y  <> 0.
  intros x y H H1 H2; elim (Zmult_integral _ _ H2); intuition.
Qed.

Lemma Zmult_non_0_r :
  forall x y, x*y  <> 0 -> y <> 0.
  intros x y H H1; apply H; rewrite H1; ring.
Qed.

Lemma Zpower_not_0 :
  forall x n, 0 <= n -> x <> 0 -> x ^ n <> 0.
  intros x n; elim n using (well_founded_ind (Zwf_well_founded 0)).
  clear n; intros n Hrec Hpos Hn0; elim (Z_le_gt_dec n 0).
  intros Hneg.
  replace n with 0; simpl; omega.
  intros Hgt; replace n with ((n-1)+1);[idtac|ring].
  rewrite Zpower_exp; try omega.
  replace (x ^ 1) with x.
  apply Zmult_non0.
  apply Hrec.
  unfold Zwf; omega.
  omega.
  assumption.
  assumption.
  unfold Zpower, Zpower_pos, iter_pos; ring.
Qed.

Theorem bounds_exp_non_0 :
  forall n s k p, bounds n s = (p,k) -> 0 <= k.
induction n.
  simpl.
  intros s k p Heq; injection Heq; intros; omega.
  simpl.
  intros s k p.
  destruct s as [d tl].
  generalize (IHn tl).
  destruct (bounds n tl) as [[a b] e].
  intros IHn'; generalize (IHn' e (a,b) (refl_equal (a,b,e))).
  intros IHn''; destruct d; intros Heq; injection Heq; intros; omega.
Qed.

Theorem si_un_hd_lift :
  forall d s n, si_un (Cons d s) (S n) = lift d (si_un s n).
  intros d s n; simpl.
  case d.
  unfold si_un; simpl.
  generalize (bounds_exp_non_0 n s).
  destruct (bounds n s) as [[a b] k].
  intros Hen0; generalize (Hen0 k (a,b) (refl_equal (a,b,k))).
  intros Hen0'.
  rewrite Zpower_exp.
  repeat (rewrite plus_IZR || rewrite mult_IZR).
  replace (IZR (2^1)) with 2%R.
  field.
  assert (2^k <> 0).
  apply Zpower_not_0; [assumption | omega].
  assert (2 <> 0)%R.
  apply (not_O_IZR 2); omega.
  repeat
   (apply prod_neq_R0 || apply not_O_IZR); auto.
  auto.
  omega.
  omega.
  unfold si_un; simpl.
  generalize (bounds_exp_non_0 n s).
  destruct (bounds n s) as [[a b] k].
  intros Hen0; generalize (Hen0 k (a,b) (refl_equal (a,b,k))).
  intros Hen0'.
  rewrite Zpower_exp.
  repeat (rewrite plus_IZR || rewrite mult_IZR).
  replace (IZR (2^1)) with 2%R.
  field.
  assert (2^k <> 0).
  apply Zpower_not_0; [assumption | omega].
  assert (2 <> 0)%R.
  apply (not_O_IZR 2); omega.
  repeat
   (apply prod_neq_R0 || apply not_O_IZR); auto.
  auto.
  omega.
  omega.

  unfold si_un.
    lazy beta iota zeta delta [bounds]; fold bounds.
  generalize (bounds_exp_non_0 n s).
  destruct (bounds n s) as [[a b] k].
  intros Hen0; generalize (Hen0 k (a,b) (refl_equal (a,b,k))).
  intros Hen0'.
  rewrite Zpower_exp.
  repeat (rewrite plus_IZR || rewrite mult_IZR).
  replace (IZR 2) with 2%R.
  replace (IZR (2^2)) with 4%R.
  replace (IZR (2^1)) with 2%R.
  unfold lift.
  field.
  assert (2^k <> 0).
  apply Zpower_not_0; [assumption | omega].
  assert (2 <> 0)%R.
  apply (not_O_IZR 2); omega.

  repeat
   (apply prod_neq_R0 || apply not_O_IZR); auto.
  auto.
  simpl; field.
  simpl; field.
  omega.
  omega.
Qed.

Theorem lift_monotonic :
  forall d x y, (x <= y)%R -> (lift d x <= lift d y)%R.
intros [ | | ] x y Hle; simpl; fourier.
Qed.

Theorem si_un_increase :
  forall n s, (si_un s n <= si_un s (n+1))%R.
  induction n.
  intros [[ | | ] s]; unfold si_un; simpl; fourier.
  intros [ d s]; rewrite si_un_hd_lift.
  simpl (S n + 1)%nat.
  rewrite si_un_hd_lift.
  apply lift_monotonic; auto.
Qed.

Theorem si_un_growing :  forall s, Un_growing (si_un s).
unfold Un_growing.
intros s n.
replace (S n) with (n+1)%nat.
exact (si_un_increase n s).
ring.
Qed.

Theorem si_un_lower_bound_0 : forall n s, (0 <= si_un s n)%R.
induction n.
intros s; unfold si_un, bounds; simpl; fourier.
intros s; apply Rle_trans with (1:= (IHn s)).
exact (si_un_growing s n).
Qed.

Theorem si_un_bounded_1 :  forall n s, (si_un s n <= 1)%R.
induction n.
intros s.
unfold si_un.
unfold bounds; simpl; fourier.
destruct s as [d s].
rewrite si_un_hd_lift.
case d; unfold lift; generalize (IHn s);intros; fourier.
Qed.

Theorem si_un_bounded : forall s, has_ub (si_un s).
intros s.
unfold has_ub, EUn, bound.
exists 1%R; unfold is_upper_bound.
intros x [n Heq]; rewrite Heq; clear Heq.
apply si_un_bounded_1.
Qed.

Definition real_value (s:stream idigit) :=
  let (v,_) := growing_cv (si_un s) (si_un_growing s) (si_un_bounded s) in
  v.

(* maps any x to 1-x *)

CoFixpoint complement (s:stream idigit): stream idigit :=
  match s with
    Cons d s' => 
    Cons (match d with L => R | R => L | C => C end) (complement s')
  end.

(* constructs the infinite stream representing a rational number.
   This function works even if the rational number is not between
   0 and 1, but the value is meaningless. *)

CoFixpoint rat_to_stream (a b:Z) : stream idigit :=
  if Z_le_gt_dec (2*a) b then
    Cons L (rat_to_stream (2*a) b)
  else
    Cons R (rat_to_stream (2*a-b) b).

Record affine_data : Set :=
 {m_a : Z; m_a' : Z; m_b : Z; m_b' : Z; m_c : Z; m_c' : Z;
  m_v1 : stream idigit; m_v2 : stream idigit}.

Definition positive_coefficients (x:affine_data) :=
  0 <= m_a x /\ 0 < m_a' x /\ 0 <= m_b x /\ 0 < m_b' x
  /\ 0 <= m_c x /\ 0 < m_c' x.

Definition af_real_value (x:affine_data) : Rdefinitions.R :=
      (IZR (m_a x)/IZR (m_a' x) * (real_value (m_v1 x)) +
       IZR (m_b x)/IZR (m_b' x) * (real_value (m_v2 x)) +
       IZR (m_c x)/IZR (m_c' x))%R.

Module Type axbyc_type.

Parameter axbyc : 
   forall x: affine_data, positive_coefficients x ->
   stream idigit.

Parameter
 axbyc_correct : forall x:affine_data, forall H :positive_coefficients x,
   (0 <= af_real_value x <= 1)%R -> real_value (axbyc x H) = af_real_value x.

End axbyc_type.

Definition axbyc_consume (x:affine_data) : affine_data :=
  let (a,a',b,b',c,c',v1,v2) := x in
  let (d1,v1') := v1 in let (d2,v2') := v2 in
  let (c1,c1') :=
      match d1,d2 with
      | L,L => (c, c')
      | L,R => (b*c'+2*c*b', 2*b'*c')
      | R,L => (a*c'+2*c*a', 2*a'*c')
      | L,C => (b*c'+4*c*b', 4*b'*c')
      | C,L => (a*c'+4*c*a', 4*a'*c')
      | R,C => (2*a*b'*c'+b*a'*c'+4*c*a'*b', 4*a'*b'*c')
      | C,R => (2*b*a'*c'+a*b'*c'+4*c*b'*a', 4*b'*a'*c')
      | R,R => (a*b'*c'+b*a'*c'+2*c*a'*b', 2*a'*b'*c')
      | C,C => (b*a'*c'+a*b'*c'+4*c*b'*a', 4*b'*a'*c')
      end in
   Build_affine_data a (2*a') b (2*b') c1 c1' v1' v2'.
      
Definition affine_measure (x:affine_data) :=
  Zmax 0 (Zmax (8 * m_a x - m_a' x) (8 * m_b x - m_b' x)).

Definition order : affine_data -> affine_data -> Prop :=
  fun x y => Zwf 0 (affine_measure x) (affine_measure y).

Definition prod_R x :=
  Build_affine_data (2*m_a x) (m_a' x) (2*m_b x) (m_b' x)
  (2*m_c x-m_c' x) (m_c' x) (m_v1 x) (m_v2 x).

Definition prod_L x :=
  Build_affine_data (2*m_a x) (m_a' x) (2*m_b x) (m_b' x)
  (2*m_c x) (m_c' x) (m_v1 x) (m_v2 x).

Definition prod_C x :=
  Build_affine_data (2*m_a x) (m_a' x) (2*m_b x) (m_b' x)
  (4*m_c x-m_c' x) (2*m_c' x) (m_v1 x) (m_v2 x).

Definition same_value (a b:affine_data) : Prop :=
  af_real_value a = af_real_value b.

Module Type axbyc_proofs.

  Parameter axbyc_test :
  forall x,
   positive_coefficients x ->
   {m_c' x <= 2*m_c x}+
   {2*(m_a x*m_b' x*m_c' x +
        m_b x*m_a' x*m_c' x + m_a' x*m_b' x*m_c x)<=
     m_a' x*m_b' x*m_c' x}+
   {4*(m_a x*m_b' x*m_c' x +
        m_b x*m_a' x*m_c' x + m_a' x*m_b' x*m_c x)<=
     3*m_a' x*m_b' x*m_c' x /\ m_c' x <= 4*m_c x}+
   {m_a' x < 8*m_a x \/ m_b' x < 8*m_b x}.

  Parameter order_wf : well_founded order.

  Parameter axbyc_consume_decrease :
     forall x, positive_coefficients x ->
         m_a' x < 8*m_a x \/ m_b' x < 8*m_b x -> 
         order (axbyc_consume x) x.
  
  Parameter axbyc_consume_pos :
     forall x, positive_coefficients x ->
       positive_coefficients (axbyc_consume x).

  Parameter prod_R_pos :
     forall x, positive_coefficients x ->
        m_c' x <= 2*m_c x ->
        positive_coefficients (prod_R x).

  Parameter prod_L_pos :
     forall x, positive_coefficients x ->
        positive_coefficients (prod_L x).

  Parameter prod_C_pos :
     forall x, positive_coefficients x ->
        m_c' x <= 4*m_c x ->
        positive_coefficients (prod_C x).

  Parameter  axbyc_consume_same_value :
    forall x, positive_coefficients x -> same_value x (axbyc_consume x).

  Parameter prod_R_in_bound :
    forall x, positive_coefficients x -> (0 <= af_real_value x <= 1)%R ->
    m_c' x <= 2 * m_c x -> (0 <= af_real_value (prod_R x) <= 1)%R.

  Parameter prod_R_expand :
    forall x, positive_coefficients x ->
     ((af_real_value (prod_R x) + 1)/2 = af_real_value x)%R.

  Parameter prod_L_in_bound :
    forall x, positive_coefficients x -> (0 <= af_real_value x <= 1)%R ->
    2 * (m_a x * m_b' x * m_c' x + m_b x * m_a' x * m_c' x +
        m_a' x * m_b' x * m_c x) <= m_a' x * m_b' x * m_c' x ->
    (0 <= af_real_value (prod_L x) <= 1)%R.

  Parameter prod_L_expand :
    forall x, positive_coefficients x ->
     (af_real_value (prod_L x)/2 = af_real_value x)%R.

  Parameter prod_C_in_bound :
    forall x, positive_coefficients x -> (0 <= af_real_value x <= 1)%R ->
    4 * (m_a x * m_b' x * m_c' x + m_b x * m_a' x * m_c' x +
         m_a' x * m_b' x * m_c x) <= 3 * m_a' x * m_b' x * m_c' x ->
    m_c' x <= 4 * m_c x ->
    (0 <= af_real_value (prod_C x) <= 1)%R.

  Parameter prod_C_expand :
    forall x, positive_coefficients x ->
     ((2*af_real_value (prod_C x)+1)/4 = af_real_value x)%R.

  Parameter represents_equal :
    forall s r, represents s r -> real_value s = r.

End axbyc_proofs.

Module axbyc_def(A:axbyc_proofs) : axbyc_type.

Inductive decision_data (ref:affine_data) : Set :=
  caseR : forall x:affine_data, 
          positive_coefficients x -> 
          m_c' x <= 2*m_c x -> 
          same_value ref x -> decision_data ref
| caseL : forall x:affine_data,
          positive_coefficients x ->
          2*(m_a x*m_b' x*m_c' x +
             m_b x*m_a' x*m_c' x + m_a' x*m_b' x*m_c x)<=
           m_a' x*m_b' x*m_c' x ->
           same_value ref x ->  decision_data ref
| caseC : forall x:affine_data,
          positive_coefficients x ->
          4*(m_a x*m_b' x*m_c' x +
             m_b x*m_a' x*m_c' x + m_a' x*m_b' x*m_c x)<=
          3*m_a' x*m_b' x*m_c' x -> m_c' x <= 4*m_c x ->
          same_value ref x -> decision_data ref.

Theorem axbyc_consume_same_value2 :
  forall x b, positive_coefficients x ->
      same_value (axbyc_consume x) b -> same_value x b.
unfold same_value; intros x b Hpos Hs1;
generalize (A.axbyc_consume_same_value x Hpos); unfold same_value;
intros Hs2; rewrite Hs2; exact Hs1.
Qed.




Definition axbyc_rec_aux (x:affine_data)
   : (forall y, order y x ->
         positive_coefficients y -> decision_data y)->
     positive_coefficients x -> decision_data x :=
  fun f Hp =>
  match A.axbyc_test x Hp with
    inleft (inleft (left H)) => caseR x x Hp H (refl_equal (af_real_value x))
  | inleft (inleft (right H)) => caseL x x Hp H (refl_equal (af_real_value x))
  | inleft (inright (conj H1 H2)) => 
         caseC x x Hp H1 H2 (refl_equal (af_real_value x))
  | inright H =>
    match f (axbyc_consume x)
            (A.axbyc_consume_decrease x Hp H)
            (A.axbyc_consume_pos x Hp) with
    | caseR b Hp' Hc Hsame =>
      caseR x b Hp' Hc (axbyc_consume_same_value2 x b Hp Hsame)
    | caseL b Hp' Hc Hsame =>
      caseL x b Hp' Hc (axbyc_consume_same_value2 x b Hp Hsame)
    | caseC b Hp' Hc1 Hc2 Hsame =>
      caseC x b Hp' Hc1 Hc2 (axbyc_consume_same_value2 x b Hp Hsame)
    end
  end.

Definition axbyc_rec :=
  well_founded_induction A.order_wf
  (fun x => positive_coefficients x -> decision_data x)
  axbyc_rec_aux.

CoFixpoint axbyc (x:affine_data) 
   (h:positive_coefficients x):stream idigit :=
  match axbyc_rec x h with
    caseR y Hpos Hc _ => Cons R (axbyc (prod_R y) (A.prod_R_pos y Hpos Hc))
  | caseL y Hpos _ _ => Cons L (axbyc (prod_L y) (A.prod_L_pos y Hpos))
  | caseC y Hpos H1 H2 _ =>
            Cons C (axbyc (prod_C y) (A.prod_C_pos y Hpos H2))
  end.

Theorem axbyc_correct_aux :
  forall x:affine_data, forall H :positive_coefficients x,
  (0 <= (af_real_value x) <= 1)%R ->
   represents (axbyc x H) (af_real_value x).
cofix.
intros x H Hint; rewrite (str_dec_thm (axbyc x H)); simpl.
destruct (axbyc_rec x H) as [v Hp Hc Hs | v Hp Hc Hs | v Hp Hc1 Hc2 Hs].

rewrite Hs; rewrite <- A.prod_R_expand.
assert (0 <= (af_real_value (prod_R v)) <= 1)%R.
apply A.prod_R_in_bound; auto.
rewrite <- Hs; exact Hint.
constructor.
apply axbyc_correct_aux;auto.
auto.
auto.

rewrite Hs; rewrite <- A.prod_L_expand.
assert (0 <= (af_real_value (prod_L v)) <= 1)%R.
apply A.prod_L_in_bound; auto.
rewrite <- Hs; exact Hint.
constructor.
apply axbyc_correct_aux;auto.
auto.
auto.

rewrite Hs; rewrite <- A.prod_C_expand.
assert (0 <= (af_real_value (prod_C v)) <= 1)%R.
apply A.prod_C_in_bound; auto.
rewrite <- Hs; exact Hint.
constructor.
apply axbyc_correct_aux;auto.
auto.
auto.
Qed.

Theorem axbyc_correct :
  forall x:affine_data, forall H :positive_coefficients x,
   (0 <= af_real_value x <= 1)%R -> real_value (axbyc x H) = af_real_value x.
intros x H Hint; apply A.represents_equal;
apply axbyc_correct_aux; assumption.
Qed.
End axbyc_def.


Module done_proofs : axbyc_proofs.

Theorem IZR4 : 4%R=IZR 4.
simpl; field.
Qed.

Theorem Rdiv_to_mult : forall a b, (b <> 0 -> a / b = a* /b)%R.
intros; field.
assumption.
Qed.

Theorem prod_R_in_bound_sup :
    forall x, (0 <= x <= 1)%R -> (2*x-1 <= 1)%R.
intros x [H1 H2]; fourier.
Qed.

Theorem prod_C_expand :
    forall x, positive_coefficients x ->
     ((2*af_real_value (prod_C x)+1)/4 = af_real_value x)%R.
intros [a a' b b' c c' v1 v2];
  unfold positive_coefficients, af_real_value, prod_C, m_a, m_a', m_b, m_b',
  m_c, m_c', m_v1, m_v2.
intros H; decompose [and] H; repeat (rewrite mult_IZR || rewrite <- Z_R_minus).
simpl (IZR 2).
simpl (IZR 4).
field.
repeat (apply prod_neq_R0|| apply prod_lt_neq_0|| apply Rlt_neq_0 || 
  apply (IZR_lt 0)); auto.
fourier.
fourier.
fourier.
Qed.

Theorem prod_L_expand :
    forall x, positive_coefficients x ->
     (af_real_value (prod_L x)/2 = af_real_value x)%R.
intros [a a' b b' c c' v1 v2];
  unfold positive_coefficients, af_real_value, prod_L, m_a, m_a', m_b, m_b',
  m_c, m_c', m_v1, m_v2.
intros H; decompose [and] H; repeat (rewrite mult_IZR || rewrite <- Z_R_minus).
simpl (IZR 2).
field.
repeat (apply prod_neq_R0|| apply prod_lt_neq_0|| apply Rlt_neq_0 || 
  apply (IZR_lt 0)); auto.
fourier.
Qed.


Theorem prod_C_in_bound_sup :
  forall a b, 0 < b -> 4*a <= 3*b -> (((4*IZR a)/IZR b -1)/IZR 2 <= 1)%R.
intros a b Hneq Hle; apply Rmult_le_reg_l with (2*IZR b)%R.
apply Rmult_lt_0_compat.
fourier.
apply (IZR_lt 0); assumption.
replace ((2*IZR b)*1)%R with (2*IZR b)%R.
replace ((2*IZR b)*((4*IZR a/IZR b -1)/IZR 2))%R with
  (4*IZR a - IZR b)%R.
replace 4%R with (IZR 4).
replace 2%R with (IZR 2).
repeat (rewrite Z_R_minus || rewrite <- mult_IZR).
apply IZR_le; omega.
simpl IZR;field.
simpl IZR;field.
simpl (IZR 2).
field.
apply prod_neq_R0.
apply Rlt_neq_0; apply (IZR_lt 0);assumption.
apply (not_O_IZR 2); omega.
field.
Qed.

Theorem prod_L_in_bound_sup :
    forall a b, 0 < b -> 2*a <= b -> ((2*IZR a)/IZR b <= 1)%R.
intros a b Hneq Hle; apply Rmult_le_reg_l with (IZR b).
apply (IZR_lt 0); assumption.
replace (IZR b * 1)%R with (IZR b).
replace (IZR b*((2 *IZR a)/IZR b))%R with (IZR (2*a)).
apply IZR_le;assumption.
rewrite mult_IZR;simpl (IZR 2).
field.
apply Rlt_neq_0; apply (IZR_lt 0);assumption.
field.
Qed.

Theorem prod_R_expand :
    forall x, positive_coefficients x ->
     ((af_real_value (prod_R x) + 1)/2 = af_real_value x)%R.
intros [a a' b b' c c' v1 v2];
  unfold positive_coefficients, af_real_value, prod_R, m_a, m_a', m_b, m_b',
  m_c, m_c', m_v1, m_v2.
intros H; decompose [and] H; repeat (rewrite mult_IZR || rewrite <- Z_R_minus).
simpl (IZR 2).
field.
repeat (apply prod_neq_R0|| apply prod_lt_neq_0|| apply Rlt_neq_0 || 
  apply (IZR_lt 0)); auto.
fourier.
Qed.

Theorem represent_in01 :
  forall x r, represents x r -> (0 <= r <= 1)%R.
intros x r Hr.
assert (H2notzero : (2<>0)%R).
apply (not_O_IZR 2); omega.
inversion Hr.
decompose [and] H0.
split; fourier.
decompose [and] H0.
replace ((r0+1)/2)%R with (r0/2+1/2)%R.
split; fourier.
field;assumption.
decompose [and] H0.
replace ((2*r0+1)/4)%R with (r0/2+1/4)%R.
split; fourier.
field.
replace 8%R with (IZR 8).
apply (not_O_IZR 8); omega.
simpl IZR.
field.
Qed.

Theorem Rdiv_le_0_compat :
  forall x y, (0 <= x -> 0 < y  -> 0 <= x / y)%R.
intros x y H H1; replace (x / y)%R with (x * /y)%R.
apply Rmult_le_pos.
assumption.
assert (0 < /y)%R.
apply Rinv_0_lt_compat;auto.
fourier.
field.
apply Rlt_neq_0;auto.
Qed.

Theorem represent_diff_2pow_n :
  forall n, forall x r1 r2, represents x r1 -> represents x r2 ->
        (-1/(2^n) <= r1 - r2 <= 1/(2^n))%R.
induction n.
simpl.
intros x r1 r2 Hr1 Hr2.
generalize (represent_in01 x r1 Hr1) (represent_in01 x r2 Hr2).
intros [H1 H2] [H3 H4].
split; fourier.
assert (H0lt2 : (0 < 2)%R).
fourier.
assert (H2n0: (2<> 0)%R).
apply (not_O_IZR 2); omega.
assert (Hpow_nz : (2^n <> 0)%R).
apply pow_nonzero;auto.
assert (Hrepl : forall a, (a<> 0)%R -> (2*(1/(2*a))= 1/a)%R).
intros a Hna; field.
apply prod_neq_R0;auto.
apply prod_neq_R0;auto.
assert (Hrepl' : forall a, (a<> 0)%R -> (2*(-1/(2*a))= -1/a)%R).
intros a Hna; field.
apply prod_neq_R0;auto.
apply prod_neq_R0;auto.
intros x r1 r2 Hr1; case Hr1.
intros tl r' Hr' _ Hr2; inversion Hr2.
simpl.
generalize (IHn tl r' r Hr' H0).
intros [IHn1 IHn2].
split.
apply Rmult_le_reg_l with 2%R; try assumption.
rewrite Hrepl'.
replace (2*(r'/2 - r/2))%R with  ((r'-r))%R.
assumption.
field; auto.
assumption.
apply Rmult_le_reg_l with 2%R; try assumption.
rewrite (Hrepl (2^n)%R).
replace (2*(r'/2 - r/2))%R with  ((r'-r))%R.
assumption.
field; auto.
assumption.
intros tl r' Hr' _ Hr2; inversion Hr2.
simpl.
generalize (IHn tl r' r Hr' H0).
intros [IHn1 IHn2].
split.
apply Rmult_le_reg_l with 2%R; try assumption.
rewrite Hrepl';auto.
replace (2*((r'+1)/2 - (r+1)/2))%R with (r' -r)%R.
assumption.
field; auto.
apply Rmult_le_reg_l with 2%R; try assumption.
rewrite Hrepl; auto.
replace (2*((r'+1)/2 - (r+1)/2))%R with (r' -r)%R.
assumption.
field;auto.
intros tl r' Hr' _ Hr2; inversion Hr2.
simpl.
generalize (IHn tl r' r Hr' H0).
intros [IHn1 IHn2].
split.
apply Rmult_le_reg_l with 2%R; try assumption.
rewrite Hrepl';auto.
replace (2 * ((2 * r' + 1) / 4 - (2 * r + 1) / 4))%R with (r' - r)%R.
assumption.
field.
replace 4%R with (IZR 4)%R.
apply not_O_IZR; omega.
simpl; field.
apply Rmult_le_reg_l with 2%R; try assumption.
rewrite Hrepl; auto.
replace (2 * ((2 * r' + 1) / 4 - (2 * r + 1) / 4))%R with (r' - r)%R.
assumption.
field.
replace 4%R with (IZR 4)%R.
apply not_O_IZR; omega.
simpl; field.
Qed.

Theorem Zmax_le1 :
  forall a b, a <= b -> Zmax a b = b.
intros a b Hle; unfold Zmax.
generalize (Zcompare_elim (a=b) True False a b); case (Zcompare a b);
 intuition.
Qed.

Theorem Zmax_le2 :
  forall a b, b <= a -> Zmax a b = a.
intros a b Hle; unfold Zmax.
generalize (Zcompare_elim (a=b) False True a b); case (Zcompare a b);
 intuition.
Qed.

Theorem cv_ext : forall f g : nat -> Rdefinitions.R, 
                     (forall n, f n = g n) ->
                 forall l, Un_cv f l -> Un_cv g l.
intros f g Hext l; unfold Un_cv.
intros Hcv1 eps Hpos; destruct (Hcv1 eps Hpos) as [n H];exists n.
intros n'; rewrite <- Hext; auto.
Qed.

Theorem cv_shift : forall f : nat -> Rdefinitions.R, forall l,
     Un_cv (fun n => f (S n)) l -> Un_cv f l.
intros f l Hcv eps Hpos; destruct (Hcv eps Hpos) as [n H]; exists (S n).
intros n'; case n'.
intro; assert (Habs :False).
omega.
elim Habs.
intros; apply H; omega.
Qed.

Theorem real_value_lift :
  forall d s, real_value (Cons d s) = lift d (real_value s).
intros d; case d; simpl; unfold real_value.
intros s; destruct (growing_cv (si_un (Cons R s)) (si_un_growing (Cons R s))
          (si_un_bounded (Cons R s))) as [x Hcv].
destruct (growing_cv (si_un s) (si_un_growing s) (si_un_bounded s)) 
   as [y Hcvy].
replace (y/2)%R with (y*(1/2))%R.
apply UL_sequence with (si_un (Cons R s)).
auto.
assert (Hcv_half : Un_cv (fun _ => (1/2)%R) (1/2)).
intros eps Hpos; exists 0%nat; intros n _; rewrite R_dist_eq.
fourier.
apply cv_shift.
apply cv_ext with 
  (fun n => (si_un s n *(1/2) + 1/2)%R).
intros n; rewrite si_un_hd_lift; unfold lift; field.
apply (not_O_IZR 2);omega.
apply (CV_plus (fun n => (si_un s n * (1/2))%R)(fun n => (1/2)%R)).
apply (CV_mult (si_un s) (fun n => 1/2)%R); auto.
auto.
field.
apply (not_O_IZR 2); omega.

intros s; destruct (growing_cv (si_un (Cons L s)) (si_un_growing (Cons L s))
          (si_un_bounded (Cons L s))) as [x Hcv].
destruct (growing_cv (si_un s) (si_un_growing s) (si_un_bounded s)) 
   as [y Hcvy].
replace (y/2)%R with (y*(1/2))%R.
apply UL_sequence with (si_un (Cons L s)).
auto.
assert (Hcv_half : Un_cv (fun _ => (1/2)%R) (1/2)).
intros eps Hpos; exists 0%nat; intros n _; rewrite R_dist_eq.
fourier.
apply cv_shift.
apply cv_ext with 
  (fun n => (si_un s n *(1/2))%R).
intros n; rewrite si_un_hd_lift; unfold lift; field.
apply (not_O_IZR 2);omega.
apply (CV_mult (si_un s) (fun n => 1/2)%R); auto.
auto.
field.
apply (not_O_IZR 2); omega.

intros s; destruct (growing_cv (si_un (Cons C s)) (si_un_growing (Cons C s))
          (si_un_bounded (Cons C s))) as [x Hcv].
destruct (growing_cv (si_un s) (si_un_growing s) (si_un_bounded s)) 
   as [y Hcvy].
replace (y/2)%R with (y*(1/2))%R.
apply UL_sequence with (si_un (Cons C s)).
auto.
assert (Hcv_half : Un_cv (fun _ => (1/2)%R) (1/2)).
intros eps Hpos; exists 0%nat; intros n _; rewrite R_dist_eq.
fourier.
apply cv_shift.
apply cv_ext with 
  (fun n => (si_un s n *(1/2) + 1/4)%R).
intros n; rewrite si_un_hd_lift; unfold lift; field.
replace 8%R with (IZR 8).
apply not_O_IZR; omega.
simpl;field.
apply (CV_plus (fun n => (si_un s n * (1/2))%R)(fun n => (1/4)%R)).
apply (CV_mult (si_un s) (fun n => 1/2)%R); auto.
intros eps Hpos; exists 0%nat; intros n _; rewrite R_dist_eq.
fourier.
field.
apply (not_O_IZR 2);omega.
Qed.

Theorem real_value_in01 :
  forall s, (0 <= real_value s <= 1)%R.
intros s; unfold real_value;
   destruct (growing_cv (si_un s)) as [v Hcv].
split.
apply Rle_cv_lim with (fun _:nat => 0%R) (si_un s).
intros n; apply si_un_lower_bound_0.
intros eps Hpos; exists 0%nat; intros n _.
rewrite R_dist_eq; fourier.
auto.
apply Rle_cv_lim with (si_un s) (fun _:nat => 1%R).
intros n; apply si_un_bounded_1.
assumption.
intros eps Hpos; exists 0%nat; intros n _.
rewrite R_dist_eq; fourier.
Qed.

Theorem Rmult_div_assoc :
  forall a b c, (c <> 0 -> a * (b/c) = a*b /c)%R.
intros; field; auto.
Qed.

Theorem prod_C_in_bound :
    forall x, positive_coefficients x -> (0 <= af_real_value x <= 1)%R ->
    4 * (m_a x * m_b' x * m_c' x + m_b x * m_a' x * m_c' x +
         m_a' x * m_b' x * m_c x) <= 3 * m_a' x * m_b' x * m_c' x ->
    m_c' x <= 4 * m_c x ->
    (0 <= af_real_value (prod_C x) <= 1)%R.
intros [a a' b b' c c' v1 v2]; unfold positive_coefficients, af_real_value,
   prod_C, m_a, m_a',  m_b, m_b', m_c, m_c', m_v1, m_v2.
intros H Hint Hc1 Hc2; decompose [and] H.
generalize (real_value_in01 v1);intros H'; destruct H'.
generalize (real_value_in01 v2);intros H'; destruct H'.
assert (IZR a' <> 0)%R.
apply Rlt_neq_0; apply (IZR_lt 0); assumption.
assert (IZR b' <> 0)%R.
apply Rlt_neq_0; apply (IZR_lt 0); assumption.
assert (IZR c' <> 0)%R.
apply Rlt_neq_0; apply (IZR_lt 0); assumption.
assert (IZR a'*IZR b'*IZR c' <> 0)%R.
repeat apply prod_neq_R0;auto.
repeat (rewrite (mult_IZR 2) || rewrite <- Z_R_minus || rewrite (mult_IZR 4)
        || rewrite <- Rmult_div_assoc || rewrite Rmult_assoc ||
        rewrite <- Rmult_plus_distr_l);auto.
split.
unfold Rdiv.
repeat (assumption || apply Rplus_le_le_0_compat || apply Rmult_le_pos ||
          apply (IZR_le 0) || (apply Rlt_le ; apply Rinv_0_lt_compat; 
                               apply (IZR_lt 0)) || omega).
rewrite <- mult_IZR; rewrite Z_R_minus; apply (IZR_le 0); omega.
apply Rlt_le; apply Rinv_0_lt_compat; apply Rmult_lt_0_compat;apply (IZR_lt 0).
omega.
assumption.
apply Rle_trans with 
    ((4*IZR(a * b' * c' + b * a' * c' + a' * b' * c)/IZR(a'*b'*c')-1)/IZR 2)%R.
replace (IZR 2 * (IZR a / IZR a' * real_value v1 +
         IZR b / IZR b' * real_value v2) +
    (IZR 4 * IZR c - IZR c') / (IZR 2 * IZR c'))%R with
    ((4* (IZR(a*b'*c')*real_value v1 +(IZR (b*a'*c')*real_value v2)+
         IZR(a'*b'*c))/IZR(a'*b'*c') -1)/IZR 2)%R.
unfold Rdiv; apply Rmult_le_compat_r.
apply Rlt_le; apply Rinv_0_lt_compat; apply (IZR_lt 0);omega.
unfold Rminus.
apply Rplus_le_compat_r.
apply Rmult_le_compat_r.
apply Rlt_le; apply Rinv_0_lt_compat.
apply (IZR_lt 0); repeat apply Zmult_lt_O_compat; auto.
apply Rmult_le_compat_l.
fourier.
repeat (rewrite plus_IZR).
repeat apply Rplus_le_compat.
match goal with |- (_ <= ?x)%R => pattern x at 2; rewrite <- Rmult_1_r end.
apply Rmult_le_compat_l.
repeat (rewrite plus_IZR || rewrite mult_IZR); simpl (IZR 2); simpl (IZR 4).
repeat (apply Rmult_le_pos || 
        (apply Rlt_le;apply (IZR_lt 0);assumption) ||
        (apply (IZR_le 0); assumption)).
assumption.
match goal with |- (_ <= ?x)%R => pattern x at 2; rewrite <- Rmult_1_r end.
apply Rmult_le_compat_l.
repeat (rewrite plus_IZR || rewrite mult_IZR); simpl (IZR 2); simpl (IZR 4).
repeat (apply Rmult_le_pos || 
        (apply Rlt_le;apply (IZR_lt 0);assumption) ||
        (apply (IZR_le 0); assumption)).
assumption.
fourier.
repeat (rewrite plus_IZR || rewrite mult_IZR); simpl (IZR 2); simpl (IZR 4).
field.
repeat (assumption || apply prod_neq_R0 || (apply (not_O_IZR 2);omega)).
apply prod_C_in_bound_sup.
repeat apply Zmult_lt_O_compat; auto.
repeat rewrite Zmult_assoc; auto.
Qed.

Theorem prod_L_in_bound :
    forall x, positive_coefficients x -> (0 <= af_real_value x <= 1)%R ->
    2 * (m_a x * m_b' x * m_c' x + m_b x * m_a' x * m_c' x +
        m_a' x * m_b' x * m_c x) <= m_a' x * m_b' x * m_c' x ->
    (0 <= af_real_value (prod_L x) <= 1)%R.
intros [a a' b b' c c' v1 v2]; unfold positive_coefficients, af_real_value,
   prod_L, m_a, m_a',  m_b, m_b', m_c, m_c', m_v1, m_v2.
intros H Hint Hc; decompose [and] H.
generalize (real_value_in01 v1);intros H'; destruct H'.
generalize (real_value_in01 v2);intros H'; destruct H'.
assert (IZR a' <> 0)%R.
apply Rlt_neq_0; apply (IZR_lt 0); assumption.
assert (IZR b' <> 0)%R.
apply Rlt_neq_0; apply (IZR_lt 0); assumption.
assert (IZR c' <> 0)%R.
apply Rlt_neq_0; apply (IZR_lt 0); assumption.
assert (IZR a'*IZR b'*IZR c' <> 0)%R.
repeat apply prod_neq_R0;auto.
repeat rewrite (mult_IZR 2).
split.
repeat (apply Rplus_le_le_0_compat||apply Rmult_le_pos||
        apply Rdiv_le_0_compat|| apply (IZR_le 0)|| apply (IZR_lt 0));
  fourier || (try omega).
simpl (IZR 2).
repeat (rewrite <- Rmult_div_assoc|| rewrite Rmult_assoc ||
        rewrite <- Rmult_plus_distr_l);auto.
assert (Hprodpos: 0 < a'*b'*c').
repeat (apply Zmult_lt_O_compat);assumption.
generalize (prod_L_in_bound_sup _ _ Hprodpos Hc).
intros Hn.
apply Rle_trans with (2 := Hn).
replace (IZR a / IZR a' * real_value v1 + IZR b / IZR b' * real_value v2 +
     IZR c / IZR c')%R with ((IZR a *IZR b'*IZR c' * real_value v1 +
      IZR b *IZR a'*IZR c' * real_value v2 +     IZR a' *IZR b' *IZR c)/ 
(IZR a' *IZR b' * IZR c'))%R.
repeat (rewrite mult_IZR || rewrite plus_IZR).
repeat rewrite Rdiv_to_mult; auto.
repeat rewrite  Rmult_assoc.
apply Rmult_le_compat_l.
fourier.
apply Rmult_le_compat_r.
repeat rewrite <- mult_IZR.
assert (0 < /IZR (a' * (b' * c')))%R.
apply Rinv_0_lt_compat;apply (IZR_lt 0).
rewrite Zmult_assoc; auto.
fourier.
apply Rplus_le_compat_r.
apply Rplus_le_compat.
repeat apply Rmult_le_compat_l; 
  match goal with
     |- (0 <= IZR ?x)%R => (apply (IZR_le 0);assumption) || 
               (assert (0 < IZR x)%R;[ apply (IZR_lt 0);assumption | fourier])
  |  |- _ => idtac
  end.
pattern (IZR c') at 2; replace (IZR c') with ((IZR c')*1)%R.
apply Rmult_le_compat_l.
assert (0 < IZR c')%R;[apply (IZR_lt 0); assumption | fourier].
auto.
field.
repeat apply Rmult_le_compat_l; 
  match goal with
     |- (0 <= IZR ?x)%R => (apply (IZR_le 0);assumption) || 
               (assert (0 < IZR x)%R;[ apply (IZR_lt 0);assumption | fourier])
  |  |- _ => idtac
  end.
pattern (IZR c') at 2; replace (IZR c') with ((IZR c')*1)%R.
apply Rmult_le_compat_l;auto.
assert (0 < IZR c')%R;[apply (IZR_lt 0); assumption | fourier].
field.
field.
repeat apply prod_neq_R0;apply Rlt_neq_0;apply (IZR_lt 0);assumption.
Qed.

Theorem prod_R_in_bound :
    forall x, positive_coefficients x -> (0 <= af_real_value x <= 1)%R ->
    m_c' x <= 2 * m_c x -> (0 <= af_real_value (prod_R x) <= 1)%R.
intros [a a' b b' c c' v1 v2]; unfold positive_coefficients, af_real_value,
   prod_R, m_a, m_a',  m_b, m_b', m_c, m_c', m_v1, m_v2.
intros H Hint Hc; decompose [and] H.
repeat rewrite (mult_IZR 2).
split.
repeat (apply Rplus_le_le_0_compat||apply Rmult_le_pos||
        apply Rdiv_le_0_compat|| apply (IZR_le 0)|| apply (IZR_lt 0));
  fourier || (try omega).
generalize (real_value_in01 v1);intros [H' H''];assumption.
generalize (real_value_in01 v2);intros [H' H''];assumption.
match goal with  id : (0 <= ?x <= 1)%R |- (?y <= 1)%R =>
    replace y with (2*x-1)%R
end.
apply prod_R_in_bound_sup;auto.
rewrite <- Z_R_minus.
rewrite mult_IZR; simpl (IZR 2).
field.
repeat apply prod_neq_R0; apply Rlt_neq_0; apply (IZR_lt 0);assumption.
Qed.


Theorem represents_real_value :
 forall s, represents s (real_value s).
idtac "mark1".
cofix.
intros [[ | | ] s].
rewrite real_value_lift; unfold lift.
replace (real_value s / 2 + 1 / 2)%R with ((real_value s +1)/2)%R.
apply reprR.
apply represents_real_value.
apply real_value_in01.
field.
apply (not_O_IZR 2); omega.

rewrite real_value_lift; unfold lift.
apply reprL.
apply represents_real_value.
apply real_value_in01.

rewrite real_value_lift; unfold lift.
replace (real_value s / 2 + 1 / 4)%R with ((2*real_value s +1)/4)%R.
apply reprC.
apply represents_real_value.
apply real_value_in01.
field.
replace (4*2)%R with (IZR 8).
apply not_O_IZR; omega.
simpl;field.
Qed.

Theorem smaller_half_is_zero : forall a, (0 < a -> a <= a /2 -> a = 0)%R.
intros; fourier.
Qed.

Theorem represents_equal :
  forall s r, represents s r -> real_value s = r.
intros s r Hrep.
assert (Happrox : (forall v, 0 < v -> exists n, 1/(2^n) <= v/2)%R).
intros v Hvpos;exists (Zabs_nat (up (1/v))+1)%nat.
destruct (archimed (1/v)) as [Hb1 _].
rewrite pow_add.
replace (1 / (2 ^ Zabs_nat (up (1 / v))*2^1))%R
 with ((1 / (2 ^ Zabs_nat (up (1 / v))))*/2)%R.
replace (v/2)%R with (v*/2)%R.
apply Rmult_le_compat_r.
assert (0 < /2)%R.
apply Rinv_0_lt_compat;auto.
fourier.
fourier.
apply pow_half_smaller.
fourier.
field.
apply (not_O_IZR 2); omega.
replace (2^1)%R with 2%R.
field.
assert (2<>0)%R.
apply (not_O_IZR 2); omega.
assert (2^Zabs_nat(up(1*/v)) <> 0)%R.
apply pow_nonzero;auto.
repeat apply prod_neq_R0;auto.
simpl; field.
destruct (total_order_T (real_value s) r) as [[Hlt1 | Heq] | Hlt2].
assert (r - real_value s <= (r - real_value s)/2)%R.
destruct (Happrox (r - real_value s)%R) as [n Hle].
fourier.
assert (-1/(2^n) <= r - real_value s <= 1/(2^n))%R.
apply represent_diff_2pow_n with s.
apply Hrep.
apply represents_real_value.
decompose [and] H.
apply Rle_trans with (1/(2^n))%R;auto.
replace r with ((r - real_value s)+real_value s)%R.
rewrite (smaller_half_is_zero (r - real_value s)); auto.
field.
fourier.
field.
assumption.

assert (real_value s -r <= (real_value s - r)/2)%R.
destruct (Happrox (real_value s - r)%R) as [n Hle].
fourier.
assert (-1/(2^n) <= real_value s - r<= 1/(2^n))%R.
apply represent_diff_2pow_n with s.
apply represents_real_value.
apply Hrep.
decompose [and] H.
apply Rle_trans with (1/(2^n))%R;auto.
replace (real_value s) with ((real_value s - r) + r)%R.
rewrite (smaller_half_is_zero (real_value s - r)); auto.
field.
fourier.
field.
Qed.

Theorem axbyc_consume_decrease :
     forall x, positive_coefficients x ->
         m_a' x < 8*m_a x \/ m_b' x < 8*m_b x -> 
         order (axbyc_consume x) x.
intros x; case x.
unfold positive_coefficients, axbyc_consume, order,
  affine_measure,
  m_a, m_a', m_b, m_b', m_c, m_c', m_v1, m_v2.
intros a a' b b' c c' [[ | | ] v1'] [[ | | ] v2'] 
   [Ha [Ha' [Hb [Hb' [Hc Hc']]]]] [Hlt | Hlt];
   case (Z_le_gt_dec (8*a-a')(8*b-b'));
   intros H1;
   case (Z_le_gt_dec (8*a-2*a')(8*b-2*b'));
   intros H2;
  ((rewrite (Zmax_le1 (8*a-a') (8*b-b'));[idtac | omega])||
   (try (rewrite (Zmax_le2 (8*a-a')(8*b-b'));[idtac | omega])));
  ((rewrite (Zmax_le1 (8*a-2*a') (8*b-2*b'));[idtac | omega])||
   (try (rewrite (Zmax_le2 (8*a-2*a')(8*b-2*b'));[idtac | omega])));
  match goal with
    |- Zwf 0 (Zmax 0 ?v) _ =>
       case (Z_le_gt_dec 0 v); intros H3;
       (rewrite (Zmax_le1 0 v); [idtac | omega]) ||
       (rewrite (Zmax_le2 0 v); [idtac | omega])
  end;
  (rewrite Zmax_le1;[idtac | omega]) ||
  (rewrite Zmax_le2;[idtac | omega]);
  unfold Zwf; omega.
Qed.

Theorem formulas_condition1 :
    forall a a' b b' c c',
      0 <= a -> 0 < a' -> 0 <= b -> 0 < b' ->
      0 <= c -> 0 < c' -> 8*a <= a' -> 8*b <= b' ->
      4*c < c' -> 
      2*(a*b'*c'+b*a'*c'+a'*b'*c) <= a'*b'*c'.
intros a a' b b' c c' Ha Ha' Hb Hb' Hc Hc' H8a H8b H4c.
assert (4*(a'*b'*c) <= a'*b'*c').
replace (4*(a'*b'*c)) with (a'*b'*(4*c)). 
auto with zarith.
ring.
assert (8*(a*b'*c') <= a'*b'*c').
replace (8*(a*b'*c')) with (8*a*b'*c').
auto with zarith.
ring.
assert (8*(b*a'*c') <= a'*b'*c').
replace (a'*b'*c') with (b'*a'*c').
replace (8*(b*a'*c')) with (8*b*a'*c').
auto with zarith.
ring.
ring.
assert (4*(a*b'*c'+b*a'*c') <= a'*b'*c').
omega.
omega.
Qed.

Theorem formulas_condition2 :
    forall a a' b b' c c',
      0 <= a -> 0 < a' -> 0 <= b -> 0 < b' ->
      0 <= c -> 0 < c' -> 8*a <= a' -> 8*b <= b' ->
      4*(a*b'*c'+b*a'*c'+a'*b'*c) > 3*a'*b'*c' ->
      c' <= 2*c.
intros a a' b b' c c' Ha Ha' Hb Hb' Hc Hc' H8a H8b Hbig.
case (Z_le_gt_dec c' (2*c)).
auto.
intros Hc'gt2c.
assert (2*a'*b'*c <= a'*b'*c').
replace (2*a'*b'*c) with ((a'*b')*(2*c)).
auto with zarith.
ring.
assert (8*(a*b'*c') <= a'*b'*c').
replace (8*(a*b'*c')) with (8*a*b'*c').
auto with zarith.
ring.
assert (8*(b*a'*c') <= a'*b'*c').
replace (8*(b*a'*c')) with (8*b*a'*c').
replace (a'*b'*c') with (b'*a'*c').
auto with zarith.
ring.
ring.
assert (4*a*b'*c'+4*b*a'*c' <= a'*b'*c').
replace (4*a*b'*c') with (4*(a*b'*c')).
replace (4*b*a'*c') with (4*(b*a'*c')).
omega.
ring.
ring.
assert (4*(a*b'*c'+b*a'*c'+a'*b'*c) <= 3*a'*b'*c').
replace (4*(a*b'*c'+b*a'*c'+a'*b'*c)) with
   ((4*(a*b'*c')+4*(b*a'*c'))+2*(2*a'*b'*c)).
replace (3*a'*b'*c') with (a'*b'*c'+2*(a'*b'*c')).
omega.
ring.
ring.
omega.
Qed.

Definition axbyc_test (x:affine_data) :
   positive_coefficients x ->
   {m_c' x <= 2*m_c x}+
   {2*(m_a x*m_b' x*m_c' x +
        m_b x*m_a' x*m_c' x + m_a' x*m_b' x*m_c x)<=
     m_a' x*m_b' x*m_c' x}+
   {4*(m_a x*m_b' x*m_c' x +
        m_b x*m_a' x*m_c' x + m_a' x*m_b' x*m_c x)<=
     3*m_a' x*m_b' x*m_c' x /\ m_c' x <= 4*m_c x}+
   {m_a' x < 8*m_a x \/ m_b' x < 8*m_b x}.
intros x; case x.
unfold m_a, m_a', m_b, m_b', m_c, m_c'.
intros a a' b b' c c' v1 v2
   [Ha [Ha' [Hb [Hb' [Hc Hc']]]]].
case (Z_le_gt_dec c' (2*c)).
intros; left; left; left;assumption.
intros H2cltc'.
case (Z_le_gt_dec (2*(a*b'*c'+ b*a'*c'+a'*b'*c))
         (a'*b'*c')).
intros; left; left; right; assumption.
intros Hdenlt2num.
case (Z_le_gt_dec (4*(a*b'*c'+b*a'*c'+a'*b'*c))
          (3*a'*b'*c')).
intros H4numle3den.
case (Z_le_gt_dec c' (4*c)).
intros Hc'le4c;left;right;split;assumption.
intros H4cltc';right.
case (Z_le_gt_dec (8*a) a').
intros H8alea'.
case (Z_le_gt_dec (8*b) b').
intros H8bleb'.
assert (2*(a*b'*c'+b*a'*c'+a'*b'*c) <= a'*b'*c').
apply formulas_condition1; auto with zarith.
omega.
intros H8bgtb';right;omega.
intros H8agta';left;omega.
intros H3denlt4num;right.
case (Z_le_gt_dec (8*a) a').
intros H8alea'.
case (Z_le_gt_dec (8*b) b').
intros H8bleb'.
assert (c' <= 2*c).
apply formulas_condition2 with (9:=H3denlt4num);auto with zarith.
omega.
intros;omega.
intros;omega.
Qed.

Theorem order_wf : well_founded order.
exact (wf_inverse_image affine_data Z (Zwf 0) affine_measure
        (Zwf_well_founded 0)).
Qed.

Theorem prod_R_pos : 
  forall x, positive_coefficients x ->
      m_c' x <= 2*m_c x -> positive_coefficients (prod_R x).
intros x; case x.
unfold positive_coefficients, prod_R, m_a, m_a', m_b, m_b', m_c, m_c'.
intros a a' b b' c c' v1 v2 [Ha [Ha' [Hb [Hb' [Hc Hc']]]]] Hle.
omega.
Qed.

Theorem prod_L_pos :
     forall x, positive_coefficients x ->
        positive_coefficients (prod_L x).
intros x; case x.
unfold positive_coefficients, prod_L, m_a, m_a', m_b, m_b', m_c, m_c'.
intros a a' b b' c c' [Ha [Ha' [Hb [Hb' [Hc Hc']]]]] Hle.
omega.
Qed.

Theorem prod_C_pos : 
  forall x, positive_coefficients x ->
      m_c' x <= 4*m_c x -> positive_coefficients (prod_C x).
intros x; case x.
unfold positive_coefficients, prod_C, m_a, m_a', m_b, m_b', m_c, m_c'.
intros a a' b b' c c' v1 v2 [Ha [Ha' [Hb [Hb' [Hc Hc']]]]] Hle.
omega.
Qed.

Theorem axbyc_consume_pos :
   forall x, positive_coefficients x ->
       positive_coefficients (axbyc_consume x).
intros x; case x.
unfold positive_coefficients, axbyc_consume, m_a, m_a', m_b, m_b', m_c, m_c',
   m_v1, m_v2.
intros a a' b b' c c' [[ | | ] v1'] [[ | | ] v2'] Hpos.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
repeat (apply conj);try omega.
auto with zarith;fail.
repeat (apply Zmult_lt_O_compat); auto with zarith;fail.
Qed.

Theorem axbyc_consume_same_value :
  forall x, positive_coefficients x -> same_value x (axbyc_consume x).
idtac "starting last proof".
assert (H2nonzero : (2 <> 0)%R).
apply (not_O_IZR 2); omega.
intros [a a' b b' c c' [[ | | ] x] [[ | | ] y]];
unfold positive_coefficients, same_value, af_real_value, axbyc_consume,
  m_a, m_a', m_b, m_b', m_c, m_c', m_v1, m_v2; 
repeat rewrite real_value_lift; unfold lift;
intros Hpos; decompose [and] Hpos;
repeat (rewrite mult_IZR || rewrite plus_IZR);
replace (IZR 2) with 2%R;try assumption;
replace (IZR 4) with 4%R;try apply IZR4;
field;
repeat apply prod_neq_R0; (apply not_O_IZR; omega) || (try assumption).
Qed.

End done_proofs.