Require Import Arith Omega.

(* The subject of the exercise was to define a function that computes the sum
  of the first n natural numbers.  We don't give the solution to that exercises,
  but a solution to a similar exercise. *)

(* We use the Function command instead of Fixpoint, because this generates extra
 theorems for us.  *)

Function sum_of_squares (n:nat) : nat :=
  match n with
   0 => 0
  | S p => S p * S p + sum_of_squares p
  end.

(* The theorem sum_of_squares_equation is generated for us. *)
Check sum_of_squares_equation.

Lemma sum_of_squares_formula :
  forall n, 6*sum_of_squares n = n*(n+1)*(2*n+1).
induction n.
auto.
rewrite sum_of_squares_equation.
SearchRewrite (_ * (_ + _)).
rewrite mult_plus_distr_l.
rewrite IHn; ring.
Qed.

Function binomial (n m:nat) {struct n} : nat :=
   match n, m with
     S p, S q => binomial p q + binomial p (S q)
   | S p, 0 => 1 
   | 0, 0 => 1
   | 0, S q => 0
   end.

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

Lemma binomial_n_gt : forall n m, n < m -> binomial n m = 0.
double induction n m.
simpl.
omega.
simpl; auto.
simpl.
intros; omega.
clear n m; intros n IHn m IHm nm; simpl.
assert (nm' : m < n) by omega.
repeat rewrite IHm; try omega.
Qed.

Lemma minus_S_lt :
  forall n m, m < n -> n - m = S (n - S m).
double induction n m.
intros; omega.
intros; omega.
simpl. intros n0 IHn0 _; rewrite <- minus_n_O; trivial.
clear n m; intros n IHn m IHm .
simpl; intros nm.
apply IHm. omega.
Qed.

Lemma binomial_formula :
  forall n p, p <= n -> fact p * fact (n-p) * binomial n p = fact n.
induction n.
intros p h; inversion h.
auto.
intros [ | p].
intros _; simpl; ring.
intros pn.
Search le.
assert (pn' : p <= n) by (apply le_S_n; auto).
rewrite fact_equation.
simpl minus.
rewrite binomial_equation.
rewrite mult_plus_distr_l.
match goal with |- ?a + ?b = ?c => set (u:= b) end.
repeat rewrite <- mult_assoc.
rewrite (mult_assoc (fact p)).
rewrite IHn; trivial.
unfold u; clear u.
case (le_lt_dec n p).
intros np.
assert (pn2 : n = p) by omega.
rewrite binomial_n_gt; try omega.
simpl; subst p. ring.
clear pn' pn; intros pn.
assert (H :n-p = S (n - S p)).
apply minus_S_lt; auto.
assert (Hf: fact (n-p) = (n - p) * fact (n - S p)).
rewrite H; simpl; ring.
rewrite Hf.
match goal with |- ?a + ?b = ?c => replace b with
  ((n-p)* (fact (S p) * fact (n-S p) * binomial n (S p)));[idtac | simpl;ring;fail]
end.
rewrite IHn.
rewrite <- mult_plus_distr_r.
replace (S p + (n-p)) with (1 + (p + (n - p))).
rewrite le_plus_minus_r.
simpl; ring.
omega.
ring.
omega.
Qed.

(* The next part of the file is an experiment in trying to see if the computer
 proof can help us discover the parameters a b c d, such that 

 a*sum_of_squares n = n*n*n*b + n*n*c + n*d

*)

Definition equation (a b c d n : nat) :=
  a*sum_of_squares n = n*n*n*b + n*n*c + n*d.

(* First we should notice that if the equation holds for any n, it should at least
  hold for n = 1, n = 2, n = 3. *)

Definition three_equations (a b c d:nat) :=
   equation a b c d 1 /\ equation a b c d 2 /\ equation a b c d 3.

Lemma mult_reg_l :
  forall m n p, p <> 0 -> m * p = n * p -> m = n.
induction p.
intros H; case H; auto.
intros _ q; apply le_antisym.
Search le.
apply mult_S_le_reg_l with p.
rewrite mult_comm; rewrite q; rewrite mult_comm; auto.
apply mult_S_le_reg_l with p.
rewrite mult_comm; rewrite <- q; rewrite mult_comm; auto.
Qed.

Lemma sum_of_square_formula :
  forall a b c d n, three_equations a b c d -> equation a b c d n.
unfold three_equations, equation.
SearchRewrite (1*_).
repeat match goal with |- context[?a*S ?b] =>
   let v := eval compute in (a*S b) in replace (a*S b) with v by ring
end; simpl sum_of_squares; intros a b c d n; repeat rewrite mult_1_l;
rewrite mult_1_r.
intros H; elim H; clear H; intros Ha; rewrite Ha.
intros [H].
assert (Hc : c + 3*d = 3*b).
apply plus_reg_l with (4*c); apply plus_reg_l with (2*d); apply plus_reg_l with  (5*b).
match goal with H : ?a = _ |- ?b = _ => 
  replace b with a;[rewrite H; ring | ring]
end.
clear H; replace (27*b) with (9*(3*b)) by ring; rewrite <- Hc.
assert (H': (b+c+d)*14 = 4*(3*b)+(c+d)*14+2*b) by ring; rewrite H'; clear H'.
rewrite <- Hc.
intros H.
(* At this point, we would like ring_simplify to simplify the equation H, but
  ring_simplify does not work in hypotheses. *)
match goal with  H : ?a=?b |- _ => assert (H1: a = b) end.
ring_simplify.
repeat rewrite <- plus_assoc; apply f_equal with (f:= plus (18*c));
repeat rewrite plus_assoc.
replace (30*d) with (26*d + 4 * d) by ring; apply f_equal with (f:= plus (26*d)).
(* remember b=2*d (we actually see 2*b=4*d) *)
apply plus_reg_l with (18*c); apply plus_reg_l with (26*d).
match goal with H : ?a = _ |- ?b = _ => replace b with a; [rewrite H; ring | ring] end.
assert (Hb : b = 2*d).
apply mult_reg_l with 2; try discriminate.
apply plus_reg_l with (18*c); apply plus_reg_l with (26*d).
match goal with H : ?a = _ |- ?b = _ => replace b with a; [rewrite H; ring | ring] end.
clear H H1.
rewrite Hb in Hc.
rename Hc into Hc_old.
assert (Hc : c = 3*d).
(* remember c = 3*d *)
apply plus_reg_l with (3*d).
match goal with H : ?a = _ |- ?b = _ => replace b with a; [rewrite H; ring | ring] end.
clear Hc_old.
subst b; subst c; ring_simplify.
clear Ha.
induction n.
simpl; ring.
rewrite sum_of_squares_equation.
SearchRewrite (_ * (_ + _)).
rewrite mult_plus_distr_l; rewrite IHn; ring.
Qed.