(**************************************************************)
(*   Copyright Dominique Larchey-Wendling [*]                 *)
(*                                                            *)
(*                             [*] Affiliation LORIA -- CNRS  *)
(**************************************************************)
(*      This file is distributed under the terms of the       *)
(*         CeCILL v2 FREE SOFTWARE LICENSE AGREEMENT          *)
(**************************************************************)

(** ** Binomial theorem *)

Require Import Arith Omega.

Require Import utils_tac gcd.

Set Implicit Arguments.

Section factorial.

  Fixpoint fact n := match n with 0 => 1 | S n => (S n) * fact n end.

  Fact fact_0 : fact 0 = 1.
(* Goal: @eq nat (fact (S n)) (Init.Nat.mul (S n) (fact n)) *)
  Proof. trivial. Qed.

  Fact fact_S n : fact (S n) = (S n)*fact n.
(* Goal: @eq nat (fact (S n)) (Init.Nat.mul (S n) (fact n)) *)
  Proof. trivial. Qed.

  Fact fact_gt_0 n : 0 < fact n.
  Proof.
(* Goal: lt O (fact n) *)
    unfold lt; simpl.
(* Goal: le (S O) (fact n) *)
    induction n as [ | n IHn ]; simpl; auto.
(* Goal: le (S O) (Init.Nat.add (fact n) (Init.Nat.mul n (fact n))) *)
    generalize (n*fact n); intros; omega.
  Qed.

End factorial.

Section binomial.

  Infix "<d" := divides (at level 70, no associativity).

  Hint Resolve divides_refl.

  Let fact_neq_0 n : fact n <> 0.
  Proof. generalize (fact_gt_0 n); omega. Qed.

  Fixpoint binomial n p :=
    match n, p with
      | n, 0     => 1
      | 0, S _   => 0
      | S n, S p => binomial n p + binomial n (S p)
    end.

  Fact binomial_n0 n : binomial n 0 = 1.
  Proof. destruct n; auto. Qed.

  Fact binomial_SS n p : binomial (S n) (S p) = binomial n p + binomial n (S p).
  Proof. auto. Qed.

  Fact binomial_n1 n : 1 <= n -> binomial n 1 = n.
  Proof.
(* Goal: forall _ : le (S O) n, @eq nat (binomial n (S O)) n *)
    destruct n as [ | n ]; try omega; intros _.
(* Goal: @eq nat (binomial (S n) (S p)) (Init.Nat.add (binomial n p) (binomial n (S p))) *)
    induction n as [ | n IHn ]; auto.
(* Goal: @eq nat (binomial (S (S n)) (S O)) (S (S n)) *)
    rewrite binomial_SS, IHn, binomial_n0; omega.
  Qed.

  Fact binomial_gt n : forall p, n < p -> binomial n p = 0.
  Proof.
(* Goal: forall (p : nat) (_ : lt n p), @eq nat (binomial n p) O *)
    induction n as [ | n IHn ]; intros [|] ?; simpl; auto; try omega.
(* Goal: @eq nat (Init.Nat.add (binomial n n0) (binomial n (S n0))) O *)
    do 2 (rewrite IHn; try omega).
  Qed.

  Fact binomial_nn n : binomial n n = 1.
  Proof.
(* Goal: @eq nat (binomial n n) (S O) *)
    induction n; auto; rewrite binomial_SS, binomial_gt with (p := S _); omega.
  Qed.

  Theorem binomial_thm n p : p <= n -> fact n = binomial n p * fact p * fact (n-p).
  Proof.
(* Goal: forall _ : le p n, @eq nat (binomial n p) (div (fact n) (Init.Nat.mul (fact p) (fact (Init.Nat.sub n p)))) *)
    intros H.
(* Goal: @eq nat (fact n) (Init.Nat.mul (Init.Nat.mul (binomial n p) (fact p)) (fact (Init.Nat.sub n p))) *)
    replace n with (n-p+p) at 1 2 by omega.
(* Goal: @eq nat (fact (Init.Nat.add (Init.Nat.sub n p) p)) (Init.Nat.mul (Init.Nat.mul (binomial (Init.Nat.add (Init.Nat.sub n p) p) p) (fact p)) (fact (Init.Nat.sub n p))) *)
    generalize (n-p); clear n H; intros n.
(* Goal: @eq nat (fact (Init.Nat.add n p)) (Init.Nat.mul (Init.Nat.mul (binomial (Init.Nat.add n p) p) (fact p)) (fact n)) *)
    induction on n p as IH with measure (n+p).
(* Goal: @eq nat (binomial (S n) (S p)) (Init.Nat.add (binomial n p) (binomial n (S p))) *)
    revert n p IH; intros [ | n ] [ | p ] IH; simpl plus; auto.
    + rewrite binomial_nn; simpl; omega.
    + rewrite Nat.add_0_r, binomial_n0; simpl; omega.
    + rewrite fact_S, binomial_SS.
      replace (S (n+S p)) with (S p+S n) by omega.
(* Goal: @eq nat (Init.Nat.mul (Init.Nat.add (S p) (S n)) (fact (Init.Nat.add n (S p)))) (Init.Nat.mul (Init.Nat.mul (Init.Nat.add (binomial (Init.Nat.add n (S p)) p) (binomial (Init.Nat.add n (S p)) (S p))) (fact (S p))) (fact (S n))) *)
      do 3 rewrite Nat.mul_add_distr_r; f_equal.
      * replace (n+S p) with (S n+p) by omega.
(* Goal: @eq nat (Nat.mul (S p) (fact (Init.Nat.add (S n) p))) (Nat.mul (Nat.mul (binomial (Init.Nat.add (S n) p) p) (fact (S p))) (fact (S n))) *)
        rewrite fact_S, IH; try omega; ring.
      * rewrite (fact_S n), IH; try omega; ring.
  Qed.

  Fact binomial_le n p : p <= n -> binomial n p = div (fact n) (fact p * fact (n-p)).
  Proof.
(* Goal: forall _ : le p n, @eq nat (binomial n p) (div (fact n) (Init.Nat.mul (fact p) (fact (Init.Nat.sub n p)))) *)
    intros H.
(* Goal: @eq nat (binomial n p) (div (fact n) (Init.Nat.mul (fact p) (fact (Init.Nat.sub n p)))) *)
    symmetry; apply div_prop with (r := 0).
    + rewrite binomial_thm with (p := p); auto; ring.
    + red; change 1 with (1*1); apply mult_le_compat; apply fact_gt_0.
  Qed.

  Fact binomial_sym n p : p <= n -> binomial n p = binomial n (n-p).
  Proof.
(* Goal: forall _ : le p n, @eq nat (binomial n p) (binomial n (Init.Nat.sub n p)) *)
    intros H; do 2 (rewrite binomial_le; try omega).
(* Goal: @eq nat (div (fact n) (Init.Nat.mul (fact p) (fact (Init.Nat.sub n p)))) (div (fact n) (Init.Nat.mul (fact (Init.Nat.sub n p)) (fact (Init.Nat.sub n (Init.Nat.sub n p))))) *)
    rewrite mult_comm; do 3 f_equal; omega.
  Qed.

  Fact binomial_spec n p : p <= n -> fact n = binomial n p * fact p * fact (n-p).
  Proof. apply binomial_thm. Qed.

  Fact binomial_0n n : 0 < n -> binomial 0 n = 0.
  Proof. intros; rewrite binomial_gt; auto; simpl. Qed.

  Theorem binomial_pascal n p : binomial (S n) (S p) = binomial n p + binomial n (S p).
  Proof. auto. Qed.

End binomial.
