(* ========================================================================= *)
(* Perfect number theorems.                                                  *)
(* ========================================================================= *)

needs "Library/prime.ml";;

prioritize_num();;

(* ------------------------------------------------------------------------- *)
(* The sum-of-divisors function.                                             *)
(* ------------------------------------------------------------------------- *)


let sigma = new_definition
 `sigma(n) = if n = 0 then 0 else nsum {d | d divides n} (\i. i)`;;

(* ------------------------------------------------------------------------- *)
(* Definition of perfection.                                                 *)
(* ------------------------------------------------------------------------- *)


let perfect = new_definition
 `perfect n <=> ~(n = 0) /\ sigma(n) = 2 * n`;;

(* ------------------------------------------------------------------------- *)
(* Various number-theoretic lemmas.                                          *)
(* ------------------------------------------------------------------------- *)


let ODD_POW2_MINUS1 = 
prove
 (`!k. ~(k = 0) ==> ODD(2 EXP k - 1)`,

  REPEAT STRIP_TAC THEN
  SUBGOAL_THEN `EVEN(2 EXP k) <=> EVEN((2 EXP k - 1) + 1)` MP_TAC THENL
   [AP_TERM_TAC THEN REWRITE_TAC[ARITH_RULE `k = k - 1 + 1 <=> ~(k = 0)`] THEN
    REWRITE_TAC[EXP_EQ_0; ARITH];
    ASM_REWRITE_TAC[GSYM NOT_EVEN; EVEN_ADD; EVEN_EXP; ARITH]]);;


let EVEN_ODD_DECOMP = 
prove
 (`!n. ~(n = 0) ==> ?r s. ODD s /\ n = 2 EXP r * s`,

  MATCH_MP_TAC num_WF THEN X_GEN_TAC `n:num` THEN
  MP_TAC(SPEC `n:num` EVEN_OR_ODD) THEN
  REWRITE_TAC[EVEN_EXISTS; ODD_EXISTS] THEN
  DISCH_THEN(DISJ_CASES_THEN (X_CHOOSE_THEN `m:num` SUBST_ALL_TAC)) THENL
   [DISCH_THEN(MP_TAC o SPEC `m:num`) THEN
    REWRITE_TAC[MULT_EQ_0; ARITH; ARITH_RULE `m < 2 * m <=> ~(m = 0)`] THEN
    ASM_CASES_TAC `m = 0` THEN ASM_REWRITE_TAC[] THEN
    ONCE_REWRITE_TAC[SWAP_EXISTS_THM] THEN MATCH_MP_TAC MONO_EXISTS THEN
    X_GEN_TAC `s:num` THEN DISCH_THEN(X_CHOOSE_TAC `r:num`) THEN
    EXISTS_TAC `SUC r` THEN ASM_REWRITE_TAC[EXP; GSYM MULT_ASSOC];
    REPEAT(DISCH_THEN(K ALL_TAC)) THEN EXISTS_TAC `0` THEN
    REWRITE_TAC[EXP; MULT_CLAUSES] THEN MESON_TAC[]]);;


let FINITE_DIVISORS = 
prove
 (`!n. ~(n = 0) ==> FINITE {d | d divides n}`,

  REPEAT STRIP_TAC THEN MATCH_MP_TAC FINITE_SUBSET THEN
  EXISTS_TAC `{d | d <= n}` THEN REWRITE_TAC[FINITE_NUMSEG_LE] THEN
  REWRITE_TAC[SUBSET; IN_ELIM_THM] THEN ASM_MESON_TAC[DIVIDES_LE]);;


let MULT_EQ_COPRIME = 
prove
 (`!a b x y. a * b = x * y /\ coprime(a,x)
             ==> ?d. y = a * d /\ b = x * d`,

  REPEAT STRIP_TAC THEN
  MP_TAC(SPECL [`a:num`; `x:num`; `y:num`] COPRIME_DIVPROD) THEN
  MP_TAC(SPECL [`x:num`; `a:num`; `b:num`] COPRIME_DIVPROD) THEN
  REPEAT(ANTS_TAC THENL
   [ASM_MESON_TAC[DIVIDES_REFL; DIVIDES_RMUL; COPRIME_SYM];
    REWRITE_TAC[divides] THEN STRIP_TAC]) THEN
  UNDISCH_TAC `a * b = x * y` THEN ASM_REWRITE_TAC[] THEN
  REWRITE_TAC[ARITH_RULE
    `(a * x * u = x * a * v) <=> (a * x) * u = (a * x) * v`] THEN
  REWRITE_TAC[EQ_MULT_LCANCEL; MULT_EQ_0] THEN ASM_MESON_TAC[]);;


let COPRIME_ODD_POW2 = 
prove
 (`!k n. ODD(n) ==> coprime(2 EXP k,n)`,

  SIMP_TAC[coprime; PRIME_2; DIVIDES_PRIMEPOW] THEN REWRITE_TAC[divides] THEN
  REPEAT STRIP_TAC THEN UNDISCH_TAC `ODD n` THEN ASM_REWRITE_TAC[] THEN
  SIMP_TAC[ODD_MULT; ODD_EXP; ARITH]);;


let MULT_NSUM = 
prove
 (`!s t. FINITE s /\ FINITE t
         ==> nsum s f * nsum t g =
             nsum {(x:A,y:B) | x IN s /\ y IN t} (\(x,y). f(x) * g(y))`,

(* ------------------------------------------------------------------------- *)
(* Some elementary properties of the sigma function.                         *)
(* ------------------------------------------------------------------------- *)


let SIGMA_0 = 
prove
 (`sigma 0 = 0`,

  REWRITE_TAC[sigma]);;


let SIGMA_1 = 
prove
 (`sigma(1) = 1`,

  REWRITE_TAC[sigma; DIVIDES_ONE; SET_RULE `{d | d = 1} = {1}`] THEN
  REWRITE_TAC[ARITH; NSUM_SING]);;


let SIGMA_LBOUND = 
prove
 (`!n. 1 < n ==> n + 1 <= sigma(n)`,

  REPEAT STRIP_TAC THEN
  FIRST_ASSUM(ASSUME_TAC o MATCH_MP (ARITH_RULE `1 < n ==> ~(n = 0)`)) THEN
  ASM_REWRITE_TAC[sigma] THEN MATCH_MP_TAC LE_TRANS THEN
  EXISTS_TAC `nsum {1,n} (\i. i)` THEN CONJ_TAC THENL
   [SIMP_TAC[NSUM_CLAUSES; FINITE_RULES; IN_SING; NOT_IN_EMPTY] THEN
    ASM_ARITH_TAC;
    MATCH_MP_TAC NSUM_SUBSET_SIMPLE THEN ASM_SIMP_TAC[FINITE_DIVISORS] THEN
    REWRITE_TAC[SUBSET; IN_ELIM_THM; NOT_IN_EMPTY; IN_INSERT] THEN
    MESON_TAC[DIVIDES_1; DIVIDES_REFL]]);;


let SIGMA_MULT = 
prove
 (`!a b. 1 < a /\ 1 < b ==> 1 + b + a * b <= sigma(a * b)`,

  REPEAT STRIP_TAC THEN
  EVERY_ASSUM(ASSUME_TAC o MATCH_MP (ARITH_RULE `1 < n ==> ~(n = 0)`)) THEN
  ASM_REWRITE_TAC[sigma] THEN MATCH_MP_TAC LE_TRANS THEN
  EXISTS_TAC `nsum {1,b,a*b} (\i. i)` THEN CONJ_TAC THENL
   [SIMP_TAC[NSUM_CLAUSES; FINITE_RULES; IN_INSERT; NOT_IN_EMPTY] THEN
    ONCE_REWRITE_TAC[ARITH_RULE `x = a * b <=> a * b = 1 * x`] THEN
    ASM_REWRITE_TAC[EQ_MULT_RCANCEL] THEN
    REWRITE_TAC[MULT_CLAUSES; MULT_EQ_1] THEN
    ASM_ARITH_TAC;
    ASM_REWRITE_TAC[MULT_EQ_0] THEN
    MATCH_MP_TAC NSUM_SUBSET_SIMPLE THEN
    ASM_SIMP_TAC[FINITE_DIVISORS; MULT_EQ_0] THEN
    REWRITE_TAC[SUBSET; IN_ELIM_THM; NOT_IN_EMPTY; IN_INSERT] THEN
    MESON_TAC[DIVIDES_1; DIVIDES_REFL; DIVIDES_LMUL]]);;


let SIGMA_PRIME = 
prove
 (`!p. prime(p) ==> sigma(p) = p + 1`,

  GEN_TAC THEN
  ASM_CASES_TAC `p = 0` THEN ASM_REWRITE_TAC[PRIME_0; SIGMA_0; ARITH] THEN
  ASM_CASES_TAC `p = 1` THEN ASM_REWRITE_TAC[PRIME_1; SIGMA_1; ARITH] THEN
  DISCH_TAC THEN ASM_REWRITE_TAC[sigma] THEN
  MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC `nsum {1,p} (\i. i)` THEN
  CONJ_TAC THENL
   [AP_THM_TAC THEN AP_TERM_TAC THEN
    REWRITE_TAC[EXTENSION; IN_ELIM_THM; IN_INSERT; NOT_IN_EMPTY] THEN
    ASM_MESON_TAC[prime; DIVIDES_1; DIVIDES_REFL];
    ASM_SIMP_TAC[NSUM_CLAUSES; IN_SING; FINITE_RULES; NOT_IN_EMPTY] THEN
    ARITH_TAC]);;


let SIGMA_PRIME_EQ = 
prove
 (`!p. prime(p) <=> sigma(p) = p + 1`,

  GEN_TAC THEN EQ_TAC THEN REWRITE_TAC[SIGMA_PRIME] THEN
  ONCE_REWRITE_TAC[GSYM CONTRAPOS_THM] THEN
  REWRITE_TAC[prime; DE_MORGAN_THM] THEN
  ASM_CASES_TAC `p = 1` THEN ASM_REWRITE_TAC[SIGMA_1; ARITH] THEN
  REWRITE_TAC[NOT_FORALL_THM; NOT_IMP; divides; DE_MORGAN_THM] THEN
  DISCH_THEN(X_CHOOSE_THEN `a:num` (CONJUNCTS_THEN2 MP_TAC ASSUME_TAC)) THEN
  DISCH_THEN(X_CHOOSE_THEN `b:num` SUBST_ALL_TAC) THEN
  MP_TAC(SPECL [`a:num`; `b:num`] SIGMA_MULT) THEN
  ASM_CASES_TAC `a = 0` THEN ASM_REWRITE_TAC[MULT_CLAUSES; SIGMA_0; ARITH] THEN
  ASM_CASES_TAC `b = 0` THEN ASM_REWRITE_TAC[MULT_CLAUSES; SIGMA_0; ARITH] THEN
  REPEAT(POP_ASSUM MP_TAC) THEN REWRITE_TAC[MULT_EQ_1] THEN
  ONCE_REWRITE_TAC[ARITH_RULE `a = a * b <=> a * b = a * 1`] THEN
  REWRITE_TAC[EQ_MULT_LCANCEL] THEN ARITH_TAC);;


let SIGMA_POW2 = 
prove
 (`!k. sigma(2 EXP k) = 2 EXP (k + 1) - 1`,

  GEN_TAC THEN REWRITE_TAC[sigma; EXP_EQ_0; ARITH] THEN
  MATCH_MP_TAC EQ_TRANS THEN
  EXISTS_TAC `nsum {2 EXP i | i <= k} (\i. i)` THEN CONJ_TAC THENL
   [AP_THM_TAC THEN AP_TERM_TAC THEN
    SIMP_TAC[DIVIDES_PRIMEPOW; PRIME_2; EXTENSION; IN_ELIM_THM];
    ALL_TAC] THEN
  MATCH_MP_TAC(ARITH_RULE `x + 1 = y ==> x = y - 1`) THEN
  SPEC_TAC(`k:num`,`k:num`) THEN INDUCT_TAC THEN REWRITE_TAC[LE] THENL
   [REWRITE_TAC[SET_RULE `{2 EXP i | i = 0} = {2 EXP 0}`] THEN
    REWRITE_TAC[ARITH; NSUM_SING];
    ALL_TAC] THEN
  REWRITE_TAC[SET_RULE
   `{2 EXP i | i = SUC k \/ i <= k} =
    (2 EXP (SUC k)) INSERT {2 EXP i | i <= k}`] THEN
  POP_ASSUM MP_TAC THEN
  REWRITE_TAC[SET_RULE
   `{2 EXP i | i <= k} = IMAGE (\i. 2 EXP i) {i | i <= k}`] THEN
  SIMP_TAC[NSUM_CLAUSES; FINITE_IMAGE; FINITE_NUMSEG_LE] THEN
  REWRITE_TAC[IN_IMAGE; GSYM LE_ANTISYM; LE_EXP; ARITH] THEN
  REWRITE_TAC[LE_ANTISYM; IN_ELIM_THM; UNWIND_THM1] THEN
  REWRITE_TAC[ARITH_RULE `~(SUC k <= k)`] THEN
  DISCH_TAC THEN ASM_REWRITE_TAC[GSYM ADD_ASSOC] THEN
  REWRITE_TAC[EXP; EXP_ADD; ARITH] THEN ARITH_TAC);;

(* ------------------------------------------------------------------------- *)
(* Multiplicativity of sigma, the most interesting property.                 *)
(* ------------------------------------------------------------------------- *)


let SIGMA_MULTIPLICATIVE = 
prove
 (`!a b. coprime(a,b) ==> sigma(a * b) = sigma(a) * sigma(b)`,

  REPEAT GEN_TAC THEN
  ASM_CASES_TAC `a = 0` THEN ASM_REWRITE_TAC[SIGMA_0; MULT_CLAUSES] THEN
  ASM_CASES_TAC `b = 0` THEN ASM_REWRITE_TAC[SIGMA_0; MULT_CLAUSES] THEN
  DISCH_TAC THEN ASM_REWRITE_TAC[sigma; MULT_EQ_0] THEN
  ASM_SIMP_TAC[FINITE_DIVISORS; MULT_NSUM] THEN
  REWRITE_TAC[IN_ELIM_THM] THEN MATCH_MP_TAC EQ_TRANS THEN
  EXISTS_TAC
   `nsum (IMAGE (\(x,y). x * y)
         {x,y | x divides a /\ y divides b}) (\i. i)` THEN
  CONJ_TAC THENL
   [AP_THM_TAC THEN AP_TERM_TAC THEN
    REWRITE_TAC[EXTENSION; IN_IMAGE; IN_ELIM_THM; EXISTS_PAIR_THM] THEN
    REWRITE_TAC[PAIR_EQ] THEN
    ONCE_REWRITE_TAC[TAUT `(a /\ b) /\ c <=> c /\ a /\ b`] THEN
    REWRITE_TAC[GSYM CONJ_ASSOC; RIGHT_EXISTS_AND_THM; UNWIND_THM1] THEN
    X_GEN_TAC `n:num` THEN EQ_TAC THEN REWRITE_TAC[DIVISION_DECOMP] THEN
    REWRITE_TAC[divides] THEN REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THEN
    MESON_TAC[MULT_AC];
    ALL_TAC] THEN
  W(fun (asl,w) -> MP_TAC(PART_MATCH (lhs o rand) NSUM_IMAGE (lhand w))) THEN
  REWRITE_TAC[o_DEF; ETA_AX] THEN DISCH_THEN MATCH_MP_TAC THEN
  REWRITE_TAC[FORALL_PAIR_THM; IN_ELIM_THM] THEN
  MAP_EVERY X_GEN_TAC [`w:num`; `x:num`; `y:num`; `z:num`] THEN
  REWRITE_TAC[PAIR_EQ] THEN STRIP_TAC THEN
  REPEAT(FIRST_X_ASSUM(SUBST_ALL_TAC o SYM o
    check (is_var o rand o concl))) THEN
  REWRITE_TAC[GSYM DIVIDES_ANTISYM] THEN
  ASM_MESON_TAC[COPRIME_DIVISORS; COPRIME_SYM; COPRIME_DIVPROD;
                DIVIDES_RMUL; DIVIDES_REFL; MULT_SYM]);;

(* ------------------------------------------------------------------------- *)
(* Hence the main theorems.                                                  *)
(* ------------------------------------------------------------------------- *)


let PERFECT_EUCLID = 
prove
 (`!k. prime(2 EXP k - 1) ==> perfect(2 EXP (k - 1) * (2 EXP k - 1))`,

  GEN_TAC THEN ASM_CASES_TAC `k = 0` THEN ASM_REWRITE_TAC[ARITH; PRIME_0] THEN
  DISCH_TAC THEN
  SUBGOAL_THEN `coprime(2 EXP (k - 1),2 EXP k - 1)` ASSUME_TAC THENL
   [MATCH_MP_TAC COPRIME_ODD_POW2 THEN ASM_SIMP_TAC[ODD_POW2_MINUS1];
    ALL_TAC] THEN
  ASM_SIMP_TAC[perfect; SIGMA_MULTIPLICATIVE; SIGMA_PRIME; SIGMA_POW2] THEN
  ASM_SIMP_TAC[ARITH_RULE `~(k = 0) ==> k - 1 + 1 = k`; EXP_EQ_0;
               MULT_EQ_0; ARITH] THEN
  CONJ_TAC THENL [ASM_MESON_TAC[PRIME_0]; ALL_TAC] THEN
  REWRITE_TAC[MULT_ASSOC] THEN GEN_REWRITE_TAC RAND_CONV [MULT_SYM] THEN
  AP_TERM_TAC THEN REWRITE_TAC[GSYM(CONJUNCT2 EXP)] THEN
  AP_TERM_TAC THEN UNDISCH_TAC `~(k = 0)` THEN ARITH_TAC);;


let PERFECT_EULER = 
prove
 (`!n. EVEN(n) /\ perfect(n)
       ==> ?k. prime(2 EXP k - 1) /\ n = 2 EXP (k - 1) * (2 EXP k - 1)`,

  GEN_TAC THEN MP_TAC(SPEC `n:num` EVEN_ODD_DECOMP) THEN
  ASM_CASES_TAC `n = 0` THENL
   [ASM_REWRITE_TAC[perfect]; ASM_REWRITE_TAC[]] THEN
  REWRITE_TAC[LEFT_IMP_EXISTS_THM; GSYM NOT_EVEN] THEN
  MAP_EVERY X_GEN_TAC [`r:num`; `s:num`] THEN
  DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC SUBST_ALL_TAC) THEN
  ASM_REWRITE_TAC[EVEN_EXP; EVEN_MULT; ARITH] THEN
  DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
  REWRITE_TAC[perfect] THEN
  ASM_SIMP_TAC[SIGMA_MULTIPLICATIVE; SIGMA_POW2;
               COPRIME_ODD_POW2; GSYM NOT_EVEN] THEN
  DISCH_TAC THEN EXISTS_TAC `r + 1` THEN
  REWRITE_TAC[ADD_SUB; EQ_MULT_LCANCEL] THEN REWRITE_TAC[EXP_EQ_0; ARITH] THEN
  FIRST_X_ASSUM(MP_TAC o check(is_eq o concl)) THEN
  REWRITE_TAC[MULT_ASSOC; GSYM(CONJUNCT2 EXP); ADD1] THEN
  DISCH_THEN(MP_TAC o MATCH_MP
    (REWRITE_RULE[IMP_CONJ] MULT_EQ_COPRIME)) THEN
  ANTS_TAC THENL
   [ONCE_REWRITE_TAC[COPRIME_SYM] THEN MATCH_MP_TAC COPRIME_ODD_POW2 THEN
    SIMP_TAC[ODD_POW2_MINUS1; ADD_EQ_0; ARITH_EQ];
    ALL_TAC] THEN
  DISCH_THEN(X_CHOOSE_THEN `d:num` MP_TAC) THEN
  ASM_CASES_TAC `d = 0` THEN ASM_REWRITE_TAC[MULT_CLAUSES] THENL
   [ASM_MESON_TAC[EVEN]; ALL_TAC] THEN
  ASM_CASES_TAC `d = 1` THENL
   [ASM_REWRITE_TAC[MULT_CLAUSES; SIGMA_PRIME_EQ] THEN
    DISCH_THEN(CONJUNCTS_THEN2 (ASSUME_TAC o SYM) ASSUME_TAC) THEN
    ASM_REWRITE_TAC[] THEN EXPAND_TAC "s" THEN
    MATCH_MP_TAC(GSYM SUB_ADD) THEN
    REWRITE_TAC[ARITH_RULE `1 <= n <=> ~(n = 0)`; EXP_EQ_0; ARITH];
    ALL_TAC] THEN
  DISCH_THEN(CONJUNCTS_THEN2 (ASSUME_TAC o SYM) ASSUME_TAC) THEN
  MP_TAC(SPECL [`2 EXP (r + 1) - 1`; `d:num`] SIGMA_MULT) THEN
  ASM_REWRITE_TAC[] THEN
  MATCH_MP_TAC(TAUT `a /\ ~b ==> (a ==> b) ==> c`) THEN
  REPEAT CONJ_TAC THENL
   [MATCH_MP_TAC(ARITH_RULE `2 EXP 1 < a ==> 1 < a - 1`) THEN
    REWRITE_TAC[LT_EXP; ARITH] THEN
    UNDISCH_TAC `~(r = 0)` THEN ARITH_TAC;
    MAP_EVERY UNDISCH_TAC [`~(d = 0)`; `~(d = 1)`] THEN ARITH_TAC;
    REWRITE_TAC[NOT_LE] THEN EXPAND_TAC "s" THEN
    REWRITE_TAC[RIGHT_SUB_DISTRIB; MULT_CLAUSES] THEN
    MATCH_MP_TAC(ARITH_RULE `1 * d < x * d ==> x * d < 1 + d + x * d - d`) THEN
    ASM_REWRITE_TAC[LT_MULT_RCANCEL] THEN
    MATCH_MP_TAC(ARITH_RULE `2 EXP 0 < a ==> 1 < a`) THEN
    REWRITE_TAC[LT_EXP] THEN UNDISCH_TAC `~(r = 0)` THEN ARITH_TAC]);;