(* begin hide *)
Require Import Classical_Pred_Type ClassicalChoice ClassicalEpsilon.
Require Import ssreflect ssrbool eqtype ssrnat ssrfun seq fintype bigop.
Require Import ZArith_ext.
(* end hide *)

(** SSReflect formalization of the Uniform Candy Distribution %\cite{ucd}%
   puzzle following %\cite{ucdsolution}%. Using the SSReflect
   extension %\cite{ssr}% of the Coq proof-assistant %\cite{coq}%. *)

(* begin hide *)

(** %\section{Preliminaries}% *)

Lemma filter_iota : forall n i m, (m <= i < m + n)%nat ->
  filter (xpred1 i) (iota m n) = [:: i].
Proof.
elim => [i m | n IH /= i m].
  rewrite addn0; case/andP => H1 H2.
  move: (leq_ltn_trans H1 H2); by rewrite ltnn.
case/andP => Hm Hi; case: ifP => [|m_i].
  move/eqP => ?; subst i.
  f_equal.
  transitivity (filter pred0 (iota m.+1 n)).
    apply eq_in_filter => j.
    rewrite mem_iota.
    case/andP.
    rewrite ltn_neqAle .
    case/andP.
    move/negbTE.
    by rewrite eq_sym => ->.
  by rewrite filter_pred0.
destruct i as [|i]; first by destruct m.
apply IH.
rewrite addSnnS Hi andbC /=.
by rewrite leq_eqVlt m_i /= in Hm.
Qed.

Notation "x `_ i" := (nth 0 x i).

Open Local Scope zarith_ext_scope.

Lemma iter_Zplus : forall n z, ssrnat.iter n (Zplus z) 0 = Z_of_nat n * z.
Proof. elim=> // n IH z. rewrite iterS IH Z_S; ring. Qed.

(** Properties of sutraction and division w.r.t.%\% even integers: *)

Lemma even_lt a b : Zeven_bool a -> Zeven_bool b -> a < b -> a <= b - 2.
Proof.
case/Zeven_boolP/Zeven_ex => pa Ha; subst a.
case/Zeven_boolP/Zeven_ex => pb Hb; subst b; omega.
Qed.

Lemma even_div2_mul2 z : Zeven_bool z -> z / 2 * 2 = z.
Proof.
case/Zeven_boolP/Zeven_ex => x ->.
rewrite Zmult_comm; f_equal; by rewrite Zmult_comm Z_div_mult_full.
Qed.

Lemma Zeven_bool_div2 z : Zeven_bool z -> z / 2 + z / 2 = z.
Proof.
move/Zeven_boolP/Zeven_div2 => ->.
rewrite Zmult_comm Z_div_mult_full //; ring.
Qed.

Canonical Structure Zplus_monoid := Monoid.Law Zplus_assoc Zplus_0_l Zplus_0_r.
Canonical Structure Zmult_muloid := Monoid.MulLaw Zmult_0_l Zmult_0_r.
Canonical Structure Zplus_comoid := Monoid.ComLaw Zplus_comm.
Canonical Structure Zplus_addoid := Monoid.AddLaw Zmult_plus_distr_l Zmult_plus_distr_r.

Notation "\zsum_ ( m <= i < n ) F" := (\big[Zplus/Z0]_( m <= i < n ) F) 
  (at level 41, F at level 41, i, m, n at level 50, 
    format "'[' \zsum_ ( m  <=  i  <  n ) '/  '  F ']'") : zarith_ext_scope.

Lemma big_Zle0_compat : forall b f, (forall i, (0 <= i < b)%nat -> 0 <= f i) ->
  0 <= \zsum_(0 <= i < b) f i.
Proof.
elim=> [f _| n IH f H]; first by rewrite big_nil.
rewrite big_nat_recr /=; apply Zplus_le_0_compat.
apply IH => i; case/andP => Hi Hi'.
apply H; by rewrite Hi /= ltnS leq_eqVlt Hi' orbC /=.
apply H => //; by rewrite ltnSn leq0n.
Qed.

(** If a sum of positive integers is null, then all values are null: *)

Lemma sum_pos_0_inv : forall n f, \zsum_(0 <= i < n) f i = 0 ->
  (forall i, (0 <= i < n)%nat -> 0 <= f i) -> forall i, (0 <= i < n)%nat -> f i = 0.
Proof.
elim=> // n IH f /= H H1 i Hi.
rewrite big_nat_recr in H.
apply Zplus0_inv in H; last 2 first.
  apply big_Zle0_compat => j Hj.
  apply H1.
  case/andP : Hj => _ Hj.
  by apply ltn_trans with n.
  apply H1; by rewrite ltnSn.
case : H => H2 H3.
case/orP : (orbN (i == n)) => X.
  move/eqP : X => X; by subst i.
apply IH.
- tauto.
- move=> j; case/andP => H4 h5.
  by apply H1, ltn_trans with n.
- rewrite ltnS in Hi.
  by rewrite leq0n /= ltn_neqAle X /= Hi.
Qed.

Lemma diff_0_equal : forall l, (forall i, (0 <= i < (size l).-1)%nat -> (l`_ i - l`_i.+1) ^^ 2 = 0) ->
  forall i, (0 <= i < size l)%nat -> l`_i = l`_0 .
Proof.
case=> h t // H; elim=> // i /= IH.
rewrite ltnS => Hi.
rewrite -IH; last by apply ltn_trans with (size t).
rewrite ( _ : t`_i = (h :: t)`_i.+1) //.
move/Zsquare_0_inv: (H _ Hi) => ?; omega.
Qed.

Lemma big_Zle_compat' : forall (c : seq nat) f g, (forall i, i \in c -> f i <= g i) ->
  \big[Zplus/Z0]_(i <- c) f i <= \big[Zplus/Z0]_(i <- c) g i.
Proof.
elim=> [ * | n c IH f g H]; first by rewrite !big_nil; apply Zle_refl.
rewrite !big_cons; apply Zplus_le_compat.
apply H; by rewrite in_cons eqxx.
apply IH => i Hi; apply H; by rewrite in_cons Hi orbC.
Qed.

(** The iterated sum is compatible with comparison: *)

Lemma big_Zle_compat : forall n f g, (forall i, (i < n)%nat -> f i <= g i) ->
  \zsum_(0 <= i < n) f i <= \zsum_(0 <= i < n) g i.
Proof.
move=> n f g H; rewrite /index_iota subn0.
apply big_Zle_compat' => i Hi; apply H.
by rewrite mem_iota add0n in Hi.
Qed.

(** The lists of positive integers: *)

Definition poslst := all (Zle_bool 0).

(** The maxium of a list of positive integers: *)

Notation "\zmax_ ( m <= i < n ) F" := (\big[Zmax/Z0]_( m <= i < n ) F) 
  (at level 41, F at level 41, i, m, n at level 50, 
    format "'[' \zmax_ ( m  <=  i  <  n ) '/  '  F ']'") : zarith_ext_scope.

Definition Z_max l := \zmax_(0 <= i < size l) l`_i.

Lemma Z_max_l : forall l i, (i < size l)%nat -> l`_i <= Z_max l.
Proof.
elim=> // h t IH [_ | n] /=.
rewrite /Z_max big_cons /=; by apply Zle_max_l.
rewrite ltnS => H.
rewrite /Z_max big_cons /= -/(size t) -/(iota 1 (size t)).
move: {IH}(IH _ H) => IH.
set x := \big[_/_]_(j <- _) _.
have -> : x = \zmax_(1 <= i < (size t).+1) (h :: t)`_i.
  rewrite /x {x}.
  apply congr_big => //; by rewrite /index_iota subn1.
rewrite {x} big_add1.
set x := \zmax_(0 <= i < size t) _.
have -> : x = Z_max t by rewrite /x /Z_max; apply congr_big_nat.
eapply Zle_trans; by [apply IH | apply Zle_max_r].
Qed.
  
Lemma le_Z_max_compat' : forall l, 0 <= \zmax_(0 <= i < size l) l`_i.
Proof.
elim => [ | h t IH /= ]; first by rewrite big_nil; apply Zle_refl.
rewrite big_nat_recl /=; eapply Zle_trans; by [apply IH | apply Zle_max_r].
Qed.

Lemma le_Z_max_compat : forall k l, (forall i, (i < size k)%nat -> k`_i <= Z_max l) ->
  Z_max k <= Z_max l.
Proof.
elim => [l _ | h t IH l /= H].
rewrite /Z_max big_nil; by apply le_Z_max_compat'.
rewrite {1}/Z_max big_cons.
apply Zmax_lub; first by apply H.
rewrite -/(size t) -/(iota 1 (size t)).
set x := \big[_/_]_(j <- _) _.
have -> : x = \zmax_(1 <= i < (size t).+1) (h :: t)`_i.
  rewrite /x {x}; apply congr_big => //; by rewrite /index_iota subn1.
rewrite {x} big_add1.
set x := \big[_/_]_(0 <= i < size t) _.
have -> : x = Z_max t by rewrite /x /Z_max; apply congr_big_nat.
apply IH => i Hi; by apply (H i.+1).
Qed.

Local Close Scope Z_scope.
(* end hide *)

(** %\section{Problem Description and Basic Properties}% *)

(** The child on the right and related properties: *)

Definition right n i := if i < n.-1 then i.+1 else O.

Lemma right_O n : right n n.-1 = O. by rewrite /right ltnn. Qed.

Lemma lt_right n i : i < n -> right n i < n.
Proof.
move=> Hi; rewrite /right.
case/orP : (orbN (i < n.-1)) => [X|X].
by rewrite X -(addn1 i) addnC ltn_add_sub subn1.
rewrite (negbTE X); apply/ltP; move/ltP in Hi; omega.
Qed.

(* begin hide *)
Local Open Scope Z_scope.
(* end hide *)

Lemma sum_rotate_right : forall l,  
  \zsum_(0 <= i < size l) l`_(right (size l) i) ^^ 2 = \zsum_(0 <= i < size l) l`_i ^^ 2.
Proof.
case=> [| h t]; first by rewrite !big_nil.
rewrite /= big_nat_recr /= right_O /=.
symmetry.
rewrite big_cons /= Zplus_comm -/(iota 1 (size t)).
set x := \big[_/_]_(j <- iota 1 (size t)) _.
have -> : x = \zsum_(1 <= j < (size t).+1) ((h :: t)`_j * ((h :: t)`_j * 1)).
    apply congr_big => //; by rewrite /index_iota subn1.
rewrite big_add1.
f_equal.
apply congr_big_nat => // i.
case/andP => _; case/andP => _ H /=.
by rewrite !Zmult_1_r /right /= H.
Qed.

(** The ``Increment if odd'' operation: *)

Definition incr_odd x := if Zodd_bool x then x + 1 else x.

(* begin hide *)
Lemma incr_odd_pos x : 0 <= x -> 0 <= incr_odd x.
Proof. move=> ?; rewrite /incr_odd Zplus_comm; case: Zodd_bool => //; by apply Zle_plus_trans. Qed.

Lemma incr_odd_le : forall a b, a <= b - 1 -> incr_odd a <= b.
Proof. move=> a b Hab. rewrite /incr_odd. case:ifP => // _; omega. Qed.

Lemma even_incr_odd h : Zeven_bool (incr_odd h).
Proof.
rewrite /incr_odd; case: ifP.
move/Zodd_bool_Zeven_bool; by rewrite Zeven_bool_S => ->.
by move/negP/negP => ->.
Qed.

Lemma even_incr_odd_inv h : Zeven_bool h -> incr_odd h = h.
Proof.
rewrite /incr_odd; case: ifP.
by move/Zodd_bool_Zeven_bool/negbTE => ->.
by move/negP/negP => ->.
Qed.

Lemma leq_incr_odd x : x <= incr_odd x.
Proof.
rewrite /incr_odd; case/orP : (orbN (Zodd_bool x)) => [-> |]; first by omega.
move/negbTE => ->; by apply Zle_refl.
Qed.

Lemma all_even_incr_odd_inv : forall l, all Zeven_bool l -> map incr_odd l = l.
Proof.
elim=> // h t IH /=; case/andP => H1 h2.
by rewrite IH // (even_incr_odd_inv _ H1).
Qed.
(* end hide *)

(** ``Increment if odd'' always increases the number of candies%\ldots %*)

Lemma incr_odd_inc n l P : n = size l ->
  \big[Zplus/0]_(0 <= i < n | P i) l`_i <= \big[Zplus/0]_(0 <= i < n | P i) (map incr_odd l)`_i.
Proof.
move=> ->.
rewrite -big_filter.
set x := \big[_/_]_(i <- filter _ _ ) _ ; rewrite -big_filter /x {x}.
apply big_Zle_compat' => i Hi.
rewrite (nth_map 0); last first.
  rewrite mem_filter mem_index_iota in Hi.
  case/andP : Hi => _; by case/andP.
by apply leq_incr_odd.
Qed.

(** %\ldots% and increasing is strict when not all candies' numbers are even: *)

Lemma incr_odd_strict_inc l : ~~ all Zeven_bool l ->
  \zsum_(0 <= i < size l) l`_i < \zsum_(0 <= i < size l) (map incr_odd l)`_i.
Proof.
case/allPn => x; case/(nthP 0) => i i_l Hi Hx.
rewrite (bigID (xpred1 i)) /=.
set y := _ + _. rewrite (bigID (xpred1 i)) /= /y {y}.
apply Zplus_lt_le_compat; last by apply incr_odd_inc.
rewrite -big_filter /index_iota subn0 filter_iota //. set y := \big[_/_]_(i0 <- _) _.
rewrite -big_filter filter_iota // /y {y} !big_cons !big_nil !Zplus_0_r (nth_map 0) // /incr_odd Hi.
apply Zodd_bool_Zeven_bool in Hx; rewrite Hx; omega.
Qed.

(** Compute the next number of candy for the %i%the child: *)

Definition m_next i l := l`_i / 2  + l`_(right (size l) i) / 2.

(** Distribute the candies among the children: *)

Definition distribute l := map (m_next^~ l) (iota 0 (size l)).

(* begin hide *)
Lemma size_distribute l : size (distribute l) = size l.
Proof. by rewrite /distribute size_map size_iota. Qed.
(* end hide *)

(** Distribution preserves the total number of candies: *)

Lemma distribute_preserves_total : forall l, all Zeven_bool l ->
  \zsum_(0 <= i < size l) (distribute l)`_i = \zsum_(0 <= i < size l) l`_i.
Proof.
case/lastP => // h t H. set l := rcons h t.
rewrite /distribute /m_next. set x := \zsum_(0 <= i < size l) _.
have {x}-> : x = \zsum_(0 <= i < size l)
  ((map (fun x => l`_x / 2) (iota 0 (size l)))`_i +
   (map (fun x => l`_(right (size l) x) / 2) (iota 0 (size l)))`_i) .
  apply congr_big_nat => // i; case/andP => _. case/andP => _ ?.
  by rewrite !(nth_map O) // !size_iota.
rewrite big_split /=. set x := \zsum_(0 <= i < size l) _. set y := \zsum_(0 <= i < size l) _.
have -> : y = \zsum_(0 <= i < size l) (map (fun x => l`_x / 2) (iota 0 (size l)))`_i.
  rewrite {1}/l {}/y size_rcons big_nat_recr.
  symmetry; rewrite {}/x big_cons -/(iota 1 (size h)) Zplus_comm (_ : Monoid.operator _ = Zplus) //.
  f_equal; last first.
    rewrite !(nth_map O) ?size_iota // !nth_iota // !add0n (_ : size h = (size l).-1); by [rewrite right_O | rewrite /l size_rcons].
  set x := \big[_/_]_(j <- iota 1 (size h)) _.
  have -> : x = \big[Zplus/0]_(1 <= j < (size h).+1) (map (fun x => l`_x / 2) (iota 0 (size h).+1))`_j.
    apply congr_big => //; by rewrite /index_iota subn1.
  rewrite big_add1; apply congr_big_nat => // i; case/andP => _; case/andP => _ H1.
  rewrite (nth_map O); last by rewrite size_iota ltnS H1.
  rewrite (nth_map O); last by rewrite size_iota ltnS; apply ltnW.
  rewrite nth_iota; last by rewrite ltnS.
  rewrite add0n nth_iota; by [rewrite /right add0n H1 | apply ltnW; rewrite ltnS].
rewrite -big_split /=; apply congr_big_nat => // i; case/andP => _; case/andP => _ H1.
rewrite (nth_map O); last by rewrite size_iota.
rewrite nth_iota // add0n; apply Zeven_bool_div2.
move/allP : H; apply; apply/(nthP 0); by exists i.
Qed.

(** The candies are uniformly distributed: *)

Definition uniform l := all (fun x => x == l`_O) l.

(** The main algorithm: *)

Fixpoint ucd k l := if k is k'.+1 then
                      if uniform l then l else ucd k' (map incr_odd (distribute l))
                    else l. 

(* begin hide *)
Lemma size_ucd : forall k l, size (ucd k l) = size l.
Proof. 
elim=> // k IH /= l; case: ifP => // Hl.
eapply trans_eq; by [apply IH | rewrite size_map size_distribute].
Qed.

Lemma ucd_S : forall k l, ucd k.+1 l = ucd 1 (ucd k l).
Proof.
elim=> // k IH l; rewrite [k.+1]lock {1}[ucd _ _]/= -lock.
case: ifP => // Hl; by [rewrite /= !Hl | rewrite IH [ucd k.+1 l]/= Hl].
Qed.

Lemma ucd_add : forall a b l, ucd (a + b) l = ucd a (ucd b l).
Proof. elim=> // n IH b l; by rewrite addSn ucd_S IH (ucd_S n). Qed.
(* end hide *)

(** After each round, the list of numbers of candies is even: *)

Lemma ucd_even : forall k l, all Zeven_bool l -> all Zeven_bool (ucd k l).
Proof.
elim=> // n IH l Heven; rewrite ucd_S /=; case: ifP => // _; first by apply IH.
apply/allP => x; case/mapP => i [Hi ->]; by apply even_incr_odd.
Qed.

(* begin hide *)
Lemma ucd_uniform l : uniform l -> forall k, ucd k l = l.
Proof. move=> Hl; elim=> // n IH; by rewrite -addn1 ucd_add /= Hl. Qed.
(* end hide *)

(* begin hide *)
Lemma ucd_inc_S k l : all Zeven_bool l ->
  \zsum_(0 <= i < size l) (ucd k l)`_i <= \zsum_(0 <= i < size l) (ucd k.+1 l)`_i.
Proof.
move=> H; rewrite ucd_S /=.
case: ifP => // Huni; first by apply Zle_refl.
move: (distribute_preserves_total (ucd k l) (ucd_even k l H)).
rewrite size_ucd => <-; apply incr_odd_inc; by rewrite size_distribute size_ucd.
Qed.
(* end hide *)

(** The number of candies is (not strictly) increasing: *)

Lemma ucd_inc : forall d k l, all Zeven_bool l ->
  \zsum_(0 <= i < size l) (ucd k l)`_i <=  \zsum_(0 <= i < size l) (ucd (k + d) l)`_i.
Proof.
elim => [* | n IH k l Heven]; first by rewrite addn0; apply Zle_refl.
by rewrite addnS; eapply Zle_trans; [apply IH | apply ucd_inc_S].
Qed.

(* begin hide *)
Lemma poslst_ucd_1 l : poslst l -> poslst (ucd 1 l).
Proof.
move=> Hl /=; case: ifP => Huni //.
apply/allP => /= i; case/(nthP 0) => x.
rewrite 2!size_map size_iota => Hx <-.
rewrite (nth_map 0) ?size_map ?size_iota // (nth_map O) ?size_iota // /m_next nth_iota //.
apply/Zle_boolP; apply incr_odd_pos, Zplus_le_0_compat; 
  apply Z_div_pos => //; rewrite add0n; apply/Zle_boolP; move/allP : Hl; apply; apply/nthP.
by exists x.
exists (right (size l) x) => //; by apply lt_right.
Qed.

Lemma poslst_ucd : forall n l, poslst l -> poslst (ucd n l).
Proof. elim=> // ? IH ? ?. rewrite ucd_S. by apply poslst_ucd_1, IH. Qed.
(* end hide *)

(* begin hide *)
Lemma ucd_pos_S l : poslst l -> 0 <= \zsum_(0 <= i < size l) (ucd 1 l)`_i.
Proof.
move/poslst_ucd_1 => Hl.
apply big_Zle0_compat => i; case/andP => _ H.
apply/Zle_boolP; move/allP : Hl; apply.
apply/nthP; exists i => //; by rewrite size_ucd.
Qed.

Lemma ucd_pos : forall k l, poslst l -> 0 <= \zsum_(0 <= i < size l) (ucd k l)`_i.
Proof.
elim => [l /= | n IH l] Hl.
- apply big_Zle0_compat => i; case/andP => _ H1 /=.
  apply/Zle_boolP; move/allP : Hl; apply; apply/nthP; by exists i.
- rewrite ucd_S -(size_ucd n l); by apply ucd_pos_S, poslst_ucd.
Qed.
(* end hide *)

(** %\section{The Main Properties}% *)

(** %{\bf The Total Number of Candies is Bounded}% *)

(** After one round, the number of candies of one child is no
    more than what it or its neighbor had before: *)

Lemma one_round_one_child_bound l : all Zeven_bool l -> forall i, (i < size l)%nat ->
  (map incr_odd (distribute l))`_i <= Zmax l`_i l`_(right (size l) i).
Proof.
move=> Hl i Hi.
rewrite /m_next (nth_map 0); last by rewrite size_map size_iota.
rewrite (nth_map O); last by rewrite size_iota.
rewrite nth_iota // add0n.
have Hli : Zeven_bool l`_i by move/allP : Hl; apply; apply/nthP; exists i.
have Hlri : Zeven_bool l`_(right (size l) i).
  move/allP : Hl; apply; apply/nthP; exists (right (size l) i) => //; by apply lt_right.
case/orP : (orbN (l`_i == l`_(right (size l) i))) => X.
  move/eqP in X.
- rewrite /m_next -X Zeven_bool_div2 // even_incr_odd_inv // Zmax_l; by apply Zle_refl.
- have {X}[X|X] : l`_i < l`_(right (size l) i) \/ l`_(right (size l) i) < l`_i by move/eqP in X; omega.
  + have Y : l`_i <= l`_(right (size l) i) - 2 by apply even_lt.
    rewrite /m_next; set x := _ + _.
    have {Hli}Hli : x <= l`_(right (size l) i) - 1.
      rewrite /x; apply Zmult_lt_0_le_reg_r with 2 => //.
      rewrite Zmult_plus_distr_l !even_div2_mul2 //; omega.
    apply Zle_trans with l`_(right (size l) i); by [apply incr_odd_le | apply Zle_max_r].
  + have Y : l`_(right (size l) i) <= l`_i - 2  by apply even_lt.
    rewrite /m_next; set x := _ + _.
    have {Hli}Hli : x <= l`_i - 1.
      rewrite /x; apply Zmult_lt_0_le_reg_r with 2 => //.
      rewrite Zmult_plus_distr_l !even_div2_mul2 //; omega.
    apply Zle_trans with l`_i; by [apply incr_odd_le | apply Zle_max_l].
Qed.

(* begin hide *)
Lemma one_round_bound l : all Zeven_bool l -> Z_max (ucd 1 l) <= Z_max l.
Proof.
move=> even_l /=.
case: ifP => // Huni; first by apply Zle_refl.
apply le_Z_max_compat => // i.
rewrite 2!size_map size_iota // => Hi.
eapply Zle_trans; first by apply one_round_one_child_bound.
apply Zmax_lub; by [apply Z_max_l | apply Z_max_l, lt_right].
Qed.

Lemma one_child_bound_k : forall k l, all Zeven_bool l -> forall i, (i < size l)%nat ->
  (ucd k l)`_i <= Z_max l.
Proof.
elim=> [l Hl /= | n IH l Hl i Hi]; first by apply Z_max_l.
rewrite -addn1 ucd_add.
have H : all Zeven_bool (ucd 1 l) by apply ucd_even.
move: {IH}(IH _ H i); rewrite size_ucd; move/(_ Hi) => IH.
apply Zle_trans with (Z_max (ucd 1 l)) => //; by apply one_round_bound.
Qed.
(* end hide *)

(** The total sum of candies is never more than %$n$% times 
   the initial maximum for one child, where %$n$% is the total 
   number of children: *)

Lemma ucd_bound k l : all Zeven_bool l ->
  \zsum_(0 <= i < size l) (ucd k l)`_i <= Z_of_nat (size l) * Z_max l.
Proof.
move=> Hl; move: (@big_const_nat _ Z0 Zplus 0 (size l) (Z_max l)).
rewrite subn0 iter_Zplus => <-.
apply big_Zle_compat => *; by apply one_child_bound_k.
Qed.

(** %{\bf The Total Number of Candies Becomes Constant}% *)

(* begin hide *)
Local Open Scope nat_scope.
Lemma sum_stable_helper f g : (forall x, 0 <= f x) ->
  (forall x, 0 <= g x) -> (forall x, g x <= g x.+1) ->
  (forall x, g x < g (x + f x)) ->
  forall cst, exists i, cst < g (iter_nat i _ (fun x => x + f x) 0).
Proof.
move=> Hf Hg0 Hg g_f; elim => [| n [i IH]].
  exists 1%nat => /=; eapply leq_ltn_trans; by [apply Hg0 | apply g_f].
exists i.+1 => /=.
case/orP : (orbN (n.+1 == g (iter_nat i _ (fun x => x + f x) 0))) => [|X].
  move/eqP => ->; by apply g_f.
have {X}X : n.+1 < g (iter_nat i _ (fun x => x + f x) 0).
  move/ltP in IH. move/eqP in X.
  apply/ltP; omega.
eapply ltn_trans; by [apply X | apply g_f].
Qed.
Local Close Scope nat_scope.
(* end hide *)

(** After a while, the total number of candies cannot grow anymore
   (this proof uses the axiom of choice%---%a provable version in Coq) : *)

Lemma sum_stable l : poslst l -> all Zeven_bool l -> { k | forall d, 
  \zsum_(0 <= i < size l) (ucd (k + d) l)`_i = \zsum_(0 <= i < size l) (ucd k l)`_i}.
Proof.
move=> Hl Heven.
apply constructive_indefinite_description, NNPP.
move/not_ex_all_not => H.
have {H} : forall n, exists d, 
  \zsum_(0 <= i < size l) (ucd n l)`_i < \zsum_(0 <= i < size l) (ucd (n + d) l)`_i.
  move=> n; move: {H}(H n) => H.
  apply not_all_not_ex.
  contradict H => d; move: {H}(H d) (ucd_inc d n l Heven) => *; omega.
case/choice => f H.
have H4 : forall x, (Zabs_nat (\zsum_(0 <= i < size l) (ucd x l)`_i) < 
   Zabs_nat (\zsum_(0 <= i < size l) (ucd (x + f x) l)`_i))%nat.
  move=> x; apply/ltP/Zabs_nat_lt; split; by [apply ucd_pos | apply H].
have {H}H3 : (forall x, Zabs_nat (\zsum_(0 <= i < size l) (ucd x l)`_i) <=
  Zabs_nat (\zsum_(0 <= i < size l) (ucd x.+1 l)`_i))%nat.
  move=> x; apply/leP/Zabs_nat_le.
  split; by [apply ucd_pos | rewrite -addn1; apply ucd_inc].
have H2 : forall x, (0 <= Zabs_nat (\zsum_(0 <= i < size l) (ucd x l)`_i))%nat.
  move=> x; by rewrite (_ : O = Zabs_nat 0).
have H1 : forall x, (0 <= f x)%nat by move=> x; rewrite leq0n.
case: {H1 H2 H3 H4}(sum_stable_helper f 
  (fun x => Zabs_nat (\zsum_(0 <= i < size l) (ucd x l)`_i))
  H1 H2 H3 H4 (size l * Zabs_nat (Z_max l))%nat) => k Hk.
move: (ucd_bound (iter_nat k _ (fun x => x + f x)%nat O) l) => H.
set x := \big[_/_]_(0 <= i < size l) _ in H.
have X : 0 <=x by rewrite /x; apply ucd_pos.
move/Zabs_nat_le : {X H}(conj X (H Heven)).
rewrite Zabs_nat_mult Zabs_nat_Z_of_nat multE.
move/ltP in Hk; rewrite /x => ?; omega.
Qed.

(** Therefore, after a while, either the list is uniform, or
   the number of candies are all even after each new distribution step
   (in other words, no need for ``Increment if odd'' anymore): *)

Lemma even_distribute l : poslst l -> all Zeven_bool l -> { k | forall d, uniform (ucd (k + d) l) || 
  ~~ uniform (ucd (k + d) l) && all Zeven_bool (distribute (ucd (k + d) l)) }.
Proof.
move=> Hpos Heven; case: (sum_stable l Hpos Heven) => k Hk.
exists k => d.
apply/negPn/negP; rewrite negb_or; case/andP => H.
rewrite negb_and; case/orP => [|H1]; first by rewrite H.
move: {H1}(incr_odd_strict_inc _ H1); rewrite size_map size_iota => H1.
have : \zsum_(0 <= i < size l) (ucd (k + d) l)`_i < \zsum_(0 <= i < size l) (ucd (k + d.+1) l)`_i.
  rewrite size_ucd in H1.
  rewrite addnS ucd_S /= (negbTE H); eapply Zle_lt_trans; last by apply H1.
  move: (distribute_preserves_total _ (ucd_even (k + d) l Heven)).
  rewrite size_ucd => <-; by apply Zle_refl.
rewrite !Hk; by move/Zlt_irrefl.
Qed.

(** %{\bf After a While, the Sum of Squares Decreases at Each Round}% *)

(** Development of the sum of squares of differences
  %
  (i.e., $\sum_0^{n-2} (l_i - l_{i+1}) ^ 2 + (l_0 - l_{n-1}) ^ 2 
  =
  2 \, (\sum_0^{n-1} {l_i}^2 - \sum_0^{n-2} l_{i+1}  \, l_i - l_0  \, l_{n-1})
  $)%: *)

Lemma sum_squares_develop : forall n l, size l = n ->
  \zsum_(0 <= i < n.-1) (l`_i - l`_i.+1) ^^ 2 + (l`_O - l`_n.-1) ^^ 2 =
  2 * (\zsum_(0 <= i < n) (l`_i ^^ 2) - \zsum_(0 <= i < n.-1) (l`_i.+1 * l`_i) - l`_O * l`_n.-1).
Proof.
elim => [[] // _ |]; first by rewrite !big_nil.
case=> [_ [|h []] //= _ | n IH].
  by rewrite !big_cons !big_nil /= !Zmult_1_r Zplus_0_r Zminus_0_r !Zminus_diag.
case/lastP => // h t.
rewrite size_rcons; move=> [h_n].
move: {IH}(IH _ h_n) => IH.
have -> : \zsum_(0 <= i < n.+1)
  ((rcons h t)`_i - (rcons h t)`_i.+1) ^^ 2 + ((rcons h t)`_0 - (rcons h t)`_n.+1) ^^ 2 =
  \zsum_(0 <= i < n) (h`_i - h`_i.+1) ^^ 2 + (h`_0 - h`_n) ^^ 2 - (h`_0 - h`_n) ^^ 2 +
  ((rcons h t)`_n - (rcons h t)`_n.+1) ^^ 2 + ((rcons h t)`_0 - (rcons h t)`_n.+1) ^^ 2.
  rewrite big_nat_recr !nth_rcons !h_n ltnSn ltnn eqxx ltn0Sn; ring_simplify.
  have -> : \zsum_(0 <= i < n) ((rcons h t)`_i - (rcons h t)`_i.+1) ^^ 2 = 
      \zsum_(0 <= i < n) (h`_i - h`_i.+1) ^^ 2 .
    apply: congr_big_nat => // i; case/andP => _. case/andP => _ H.
    by rewrite !nth_rcons !h_n (ltn_trans H (ltnSn _)) ltnS H.
  rewrite /Zplus_monoid /=; ring.
rewrite {}IH !nth_rcons !h_n ltnn eqxx ltn0Sn ltnSn.
have -> : \zsum_(0 <= i < n.+2) (rcons h t)`_i ^^ 2 = \zsum_(0 <= i < n.+1) h`_i ^^ 2 + (rcons h t)`_n.+1 ^^ 2.
  rewrite big_nat_recr.
  have -> : \big[Zplus_monoid/0]_(0 <= i < n.+1) (rcons h t)`_i ^^ 2 =
    \big[Zplus/0]_(0 <= i < n.+1) h`_i ^^ 2.
    apply: congr_big_nat => // i; case/andP => _. case/andP => _ H.
    by rewrite !nth_rcons !h_n H.
  rewrite /Zplus_monoid /=; ring.
rewrite !nth_rcons !h_n ltnn eqxx.
have -> : \zsum_(0 <= i < n.+2.-1) ((rcons h t)`_i.+1 * (rcons h t)`_i) =
  \zsum_(0 <= i < n) (h`_i.+1 * h`_i) + (rcons h t)`_n.+1 * (rcons h t)`_ n.
  rewrite big_nat_recr.
  have -> :  \big[Zplus_monoid/0]_(0 <= i < n)((rcons h t)`_i.+1 * (rcons h t)`_i) = 
    \zsum_(0 <= i < n) (h`_i.+1 * h`_i).
    apply: congr_big_nat => // i; case/andP => _. case/andP => _ H.
    by rewrite !nth_rcons !h_n ltnS H ltnS leq_eqVlt H orbC.
  by rewrite !nth_rcons !h_n ltnn eqxx ltnSn.
rewrite !id_rem_minus -!Zpower_b_square; ring_simplify.
rewrite !nth_rcons !h_n ltnn eqxx ltnSn [(n.+1.-1)%nat]/=; ring.
Qed.

(** Therefore, after a while, the sum of squares of the number of candies
   strictly decreases at each distribution step: *)

Lemma one_iteration_decreases l : ~~ uniform l -> all Zeven_bool l -> 
  \zsum_(0 <= i < size l) (distribute l)`_i ^^ 2 < \zsum_(0 <= i < size l) l`_i ^^ 2.
Proof.
move=> Huni Heven.
apply Zmult_lt_reg_l with 4 => //; apply Zlt_0_minus_lt.
set l' := map incr_odd l.
  have -> : 4 * (\zsum_(0 <= i < size l) l`_i ^^ 2) - 4 * (\zsum_(0 <= i < size l) (distribute l)`_i ^^ 2) =
      2 * (\zsum_(0 <= i < size l) (l`_i ^^ 2) - \zsum_(0 <= i < (size l).-1) (l`_i.+1 * l`_i) - l`_O * l`_(size l).-1).
    set x1 := \zsum_(0 <= i < size l) _. set x2 := \zsum_(0 <= i < size l) _.
    have {x2}-> : x2 = \zsum_(0 <= i < size l) (l`_i / 2 + l`_(right (size l) i) / 2) ^^ 2.
      apply congr_big_nat => // i; case/andP => _ ; case/andP => _ ?.
      rewrite (nth_map O); by [rewrite nth_iota | rewrite size_iota].
    set x2 := \big[Zplus/0]_(0 <= i < size l) _.
    have {x2}-> : 4 * x2 = \zsum_(0 <= i < size l) (l`_i + l`_(right (size l) i)) ^^ 2.
      rewrite big_distrr.
      apply congr_big_nat => // i; case/andP=> _; case/andP => _ ?.
      rewrite 2!id_rem_plus.
      have [ni Hni] : exists ni, l`_i = 2 * ni.
        apply/Zeven_ex/Zeven_boolP; move/allP : Heven; apply; apply/nthP; by exists i.
      have [ni' Hni'] : exists ni, l`_(right (size l) i) = 2 * ni.
        apply/Zeven_ex/Zeven_boolP.
        move/allP : Heven; apply; apply/nthP; exists (right (size l) i) => //; by apply lt_right.
      rewrite Hni (Zmult_comm 2 ni) Hni' (Zmult_comm 2 ni') !Z_div_mult_full //.
      rewrite !Zpower_Zmult [2 ^^ 2]/= /Zmult_muloid (_: Monoid.mul_operator _ = Zmult) //; ring.
    set x2 := \zsum_(0 <= i < size l) _. have {x2}-> : x2 = \zsum_(0 <= i < size l) l`_i ^^ 2 +
      \zsum_(0 <= i < size l) l`_(right (size l) i) ^^ 2 +
      2 * \zsum_(0 <= i < size l) (l`_i * l`_(right (size l) i)).
      rewrite /x2 -big_split /Zplus_comoid (_ : Monoid.operator _ = Zplus) //.
      rewrite big_distrr -big_split /Zplus_comoid (_ : Monoid.operator _ = Zplus) //.
      apply congr_big_nat => // i; case/andP => _; case/andP => _ ?.
      rewrite id_rem_plus /Zmult_muloid (_: Monoid.mul_operator _ = Zmult) //; ring.
    rewrite /x1 !sum_rotate_right; ring_simplify.
    rewrite -Zminus_plus_distr; f_equal.
    destruct l as [|h t] => //.
    rewrite [size _]/= big_nat_recr right_O.
    have -> : \big[Zplus_monoid/0]_(0 <= i < size t) ((h :: t)`_i * (h :: t)`_(right (size (h :: t)) i)) = \zsum_(0 <= i < size t) ((h :: t)`_i.+1 * (h :: t)`_i).
      apply congr_big_nat => // i; case/andP => _; case/andP => _ H /=.
      by rewrite /right H /= Zmult_comm.
    rewrite /Zplus_monoid (_: Monoid.operator _ = Zplus) // [((size t).+1.-1)%nat]/=; ring.
  apply Zmult_lt_0_reg_r with 2 => //.
  rewrite Zmult_comm -sum_squares_develop //.
  apply Zlt_lt_double; split => //; apply Zlt0_pos.
  - apply big_Zle0_compat => i _; by apply Zsquare_pos.
  - by apply Zsquare_pos.
  - apply/negPn; move/eqP => H.
    apply Zplus0_inv in H; last 2 first.
      apply big_Zle0_compat => i _; by apply Zsquare_pos.
      by apply Zsquare_pos.
    case: H => H _.
    suff : uniform l by rewrite (negbTE Huni).
    move/sum_pos_0_inv in H.
    have : forall i, (0 <= i < (size l).-1)%nat -> 0 <= (l`_i - l`_i.+1) ^^ 2.
      move=> j _; by apply Zsquare_pos.
    move/H/(diff_0_equal l) => {H}H.
    apply/allP => x; move/nthP.
    case/(_ 0) => i Hi Hx; by rewrite -(H i) // Hx.
Qed.

(** %\section{The Termination Proof}% *)

(** The type of outputs of the ucd algorithm starting from %$l$% and after at least %$k_0$% steps: *)

Definition ucd_seq k0 l := sigT (fun s => { k | (k0 <= k)%nat * (ucd k l = s) } )%type.

(** The list inside an object of type %{\it ucd\_seq k0 l}% (projection): *)

Definition lst {k0 l} (ks : ucd_seq k0 l) := match ks with existT v _ => v end.
Notation "# l" := (lst l) (at level 50).

(* begin hide *)
Definition rank {k0 l} (ks : ucd_seq k0 l) := match ks with existT v (exist w _) => w end.
(* end hide *)

(* begin hide *)

(* the ucd_seq after %$k0$% steps: *) 

Definition ucd_seq_0 k0 l : ucd_seq k0 l. 
by exists (ucd k0 l), k0.
Defined.

(* the ucd_seq obtained after one more round: *)

Definition ucd_seq_next {k0 l} (L : ucd_seq k0 l) : ucd_seq k0 l.
case : L => x [k [H1 H2]].
exists (ucd 1 x), (k.+1)%nat; split; by [apply leq_trans with k | rewrite -H2 -ucd_S].
Defined.

(* end hide *)

(** Order relation between lists (parameterized by a minimum number of steps and the input of the ucd algorithm), defined by comparing the sum of squares: *)

Definition ss_lt k0 l (x y : ucd_seq k0 l ) :=
  \zsum_(0 <= i < size (# x)) (# x)`_i ^^ 2 < \zsum_(0 <= i < size (# y)) (# y)`_i ^^ 2.

(** This is a well-founded order: *)

Lemma well_founded_ss_lt k0 l : well_founded (ss_lt k0 l).
Proof.
apply well_founded_lt_compat with (fun l => Zabs_nat (\zsum_(0 <= i < size (# l)) (# l)`_i ^^ 2)) => x y x_y.
apply Zabs_nat_lt; split => //.
apply big_Zle0_compat => i _; by apply Zsquare_pos.
Qed.

(** Starting with a list of positive, even integers, the ucd
   algorithm converges (by well-founded induction using the preceding order): *)

Theorem ucd_terminates l : poslst l -> all Zeven_bool l -> 
  exists v, exists k1, forall k, (k1 <= k)%nat -> ucd k l = v.
Proof.
move=> Hpos Heven; case: (even_distribute l Hpos Heven) => k0 Hk0.
suff : exists v, exists k1, forall k, (k0 + k1 <= k)%nat -> ucd k l = v.
  case=> v [k1 Hv]; by exists v, (k0 + k1)%nat.
suff : { v : ucd_seq k0 l | exists k1, forall k, (k1 <= k)%nat -> ucd k (# ucd_seq_0 k0 l) = # v }.
  case=> v [k1 /= Hv]; exists (# v), k1 => k Hk.
  move: (Hv (k - k0)%nat); rewrite -(leq_add2l k0) subnKC ?(leq_trans (leq_addr _ _) Hk) //.
  move/(_ Hk); by rewrite -ucd_add subnK // -(leq_add2r k1) (leq_trans Hk) // leq_addr.
apply (well_founded_induction (well_founded_ss_lt k0 l)
  (fun x => { v | exists k1, forall k, (k1 <= k)%nat -> ucd k (# x) = # v })) => L HL.
apply constructive_indefinite_description.
case/orP : (orbN (uniform (ucd (rank L) l))) => Huni.
  exists L, k0 => k k0_k; apply ucd_uniform => /=.
  case: L HL Huni => L [kL [? /= H]] HL Huni; by rewrite /= -H.
have : ss_lt k0 l (ucd_seq_next L) L.
  rewrite /ss_lt; case: L HL Huni => L [k [H1 /= H2]] HL Huni.
  rewrite -H2 size_ucd (negbTE Huni) all_even_incr_odd_inv; last first.
    move: (Hk0 (k - k0)%nat); rewrite subnKC //.
    case/orP; by [rewrite (negbTE Huni) | case/andP].
  have -> : size l = size (ucd k l) by rewrite size_ucd.
  rewrite size_distribute; apply one_iteration_decreases => //; by apply ucd_even.
case/HL => v [k1 Hv] {HL}.
case: L Huni Hv => L [kL [? /= H]] Huni /= Hv.
exists v, (k1.+1)%nat; case => // k.
move/Hv => <-; by rewrite /= -H (negbTE Huni).
Qed.