-- Perfect numbers: σ(n) = 2n. The divisors of 2^(k+1) are 1 and the doubles
-- of the divisors of 2^k, so σ(2^k) = 2^(k+1) − 1; Euclid's theorem that
-- 2^(p−1) · (2^p − 1) is perfect when 2^p − 1 is prime, and Euler's
-- converse that every even perfect number has this form
-- Part of Cumulus, the Grace17 mathematics library.
-- Copyright (c) 2026 Grant Passmore.
-- Released under the Apache License 2.0. See LICENSE in the repository root.
-- Authors: Grant Passmore

module Grace.Cumulus.NumberTheory.PerfectNumbers
import DivisorSum
import Nat.ParityLemmas
import Nat.SubtractionLemmas

-- The doubles of a set of naturals, and doubling as a bijection onto them.
def doubles (S : Set) : Set := { n ∈ ℕ | ∃ d ∈ S, n = mul(2, d) }

defgraph doubleMap (S : Set) (d ∈ S) : doubles(S) := mul(2, d)
requires naturals : S ⊆ ℕ
proof
  have dIn := naturals with d, member
  have dNat := nat.omega.backward with dIn
  have doubleNat : Nat(mul(2, d)) by typecheck
  have doubleIn := nat.omega.forward with doubleNat
  have same : mul(2, d) = mul(2, d) by rfl
  exact doubles.membership.backward with ⟨doubleIn, d, member, same⟩

theorem doubleMap.injective (S : Set)
  requires naturals : S ⊆ ℕ
  shows InjectiveBetween(doubleMap(S), S, doubles(S))
proof
  have function := doubleMap.function at S with naturals
  refine ⟨function, _⟩
  intros x hx y hy z hz relatesX relatesY
  have xNat := nat.omega.backward with (naturals with x, hx)
  have yNat := nat.omega.backward with (naturals with y, hy)
  have valueX := FunctionBetween.apply_value with function, hx, relatesX
  have valueY := FunctionBetween.apply_value with function, hy, relatesY
  have doubleX := doubleMap.value at S, x with naturals, hx
  have doubleY := doubleMap.value at S, y with naturals, hy
  have same : mul(2, x) = mul(2, y) by cc
  have twoNonzero : 2 ≠ 0 by arith
  exact mul.left_cancel_positive with twoNonzero, same

theorem doubleMap.surjective (S : Set)
  requires naturals : S ⊆ ℕ
  shows SurjectiveBetween(doubleMap(S), S, doubles(S))
proof
  have function := doubleMap.function at S with naturals
  refine ⟨function, _⟩
  intros z hz
  have parts := doubles.membership.forward with hz
  obtain ⟨d, hd, shape⟩ := parts.right
  have relates := FunctionBetween.relates_apply with function, hd
  have value := doubleMap.value at S, d with naturals, hd
  rw [value, ← shape] at relates
  exact ⟨d, hd, relates⟩

theorem doubleMap.bijective (S : Set)
  requires naturals : S ⊆ ℕ
  shows BijectiveBetween(doubleMap(S), S, doubles(S))
proof
  have injective := doubleMap.injective with naturals
  have surjective := doubleMap.surjective with naturals
  exact ⟨injective, surjective⟩

theorem doubles.finite (S : Set)
  requires finite : Finite(S)
  requires naturals : S ⊆ ℕ
  shows Finite(doubles(S))
proof
  have bijective := doubleMap.bijective with naturals
  have equipotent : Equipotent(S, doubles(S)) := ⟨doubleMap(S), bijective⟩
  exact Finite.of_equipotent with finite, equipotent

theorem doubles.weights (S : Set)
  requires naturals : S ⊆ ℕ
  shows ∀ t ∈ doubles(S), identity(ℕ)(t) ∈ ℕ
proof
  intros t ht
  have parts := doubles.membership.forward with ht
  have tNat := nat.omega.backward with parts.left
  exact identityWeight.natural at t

-- The doubles sum to twice the sum.
theorem natSum.doubles (S : Set)
  requires finite : Finite(S)
  requires naturals : S ⊆ ℕ
  shows natSum(doubles(S), identity(ℕ)) =
    add(natSum(S, identity(ℕ)), natSum(S, identity(ℕ)))
proof
  have weights := doubles.weights with naturals
  have weightsS : ∀ s ∈ S, identity(ℕ)(s) ∈ ℕ
    intros s hs
    have sNat := nat.omega.backward with (naturals with s, hs)
    exact identityWeight.natural at s
  have function := doubleMap.function at S with naturals
  have bijective := doubleMap.bijective with naturals
  have values : ∀ s ∈ S, doubleMap(S)(s) = identity(ℕ)(doubleMap(S)(s))
    intros s hs
    have image := FunctionBetween.apply_closed with function, hs
    have parts := doubles.membership.forward with image
    symm
    exact identity.apply at doubleMap(S)(s), ℕ with parts.left
  have reindexed := natSum.reindex with finite, bijective, weights, values
  have agree : ∀ s ∈ S,
      doubleMap(S)(s) = add(identity(ℕ)(s), identity(ℕ)(s))
    intros s hs
    have sNat := nat.omega.backward with (naturals with s, hs)
    have value := doubleMap.value at S, s with naturals, hs
    have idValue := identityWeight.value at s
    rw [value, idValue]
    exact mul.two_left at s
  have correct := cardinality.correct with finite
  have cardNat := correct.left
  have size : cardinality(S) = cardinality(S) by rfl
  have pointwise := natSum.add_pointwise_of_size at doubleMap(S), identity(ℕ),
    identity(ℕ), cardinality(S)
  have added := pointwise with S, finite, size, weightsS, weightsS, agree
  rw [reindexed, added]
  rfl

-- The divisors of 2^(k+1): 1 and the doubles of the divisors of 2^k.
theorem divisors.two_pow_succ (k : Nat)
  shows divisors(2 ^ (k + 1)) = insert(1, doubles(divisors(2 ^ k)))
proof
  have twoPositive : 0 < 2 by arith
  have twoNonzero : 2 ≠ 0 by arith
  have powPositive := pow.positive at 2, k with twoPositive
  have powNonzero := Less.ne_zero_of_positive with powPositive
  have succPositive := pow.positive at 2, k + 1 with twoPositive
  have succNonzero := Less.ne_zero_of_positive with succPositive
  have expanded := pow.add_one at 2, k
  have two := Prime.two
  ext d
  constructor
  . intro member
    have dNat := divisors.natural with member
    have factor := divisors.divides with member
    apply insert.membership.backward
    have parity := parity.exhaustive at d
    cases parity with even | odd
    . right
      have characterization := Even.iff_two_mul at d
      obtain ⟨c, cNat, shape⟩ := characterization.forward with even
      rw [shape, expanded] at factor
      have cancel := divides.mul_cancel_left at 2, c, 2 ^ k with twoNonzero
      have cDivides := cancel.forward with factor
      have cBelow := divides.le with cDivides, powNonzero
      have cCharacterization := divisors.membership_le at 2 ^ k, c
      have cMember := cCharacterization.backward with ⟨cBelow, cDivides⟩
      have dIn := nat.omega.forward with dNat
      have same : 2 * c = mul(2, c) by rfl
      rw shape at dIn
      rw shape
      exact doubles.membership.backward with ⟨dIn, c, cMember, same⟩
    . left
      have oddCharacterization := divides.not_two_iff_odd at d
      have notTwo := oddCharacterization.backward with odd
      have coprimeTwo := Prime.gcd_of_not_divides with two, notTwo
      have coprimePow := coprime.pow_left at 2, d, k + 1 with coprimeTwo
      have characterization := gcd.eq_one_iff at 2 ^ (k + 1), d
      have all := characterization.forward with coprimePow
      have selfDivides := divides.reflexive at d
      exact all with d, dNat, factor, selfDivides
  . intro member
    have split := insert.membership.forward with member
    cases split with one | double
    . rw one
      exact divisors.one_member at 2 ^ (k + 1) with succNonzero
    . have parts := doubles.membership.forward with double
      obtain ⟨c, hc, shape⟩ := parts.right
      have cNat := divisors.natural with hc
      have cDivides := divisors.divides with hc
      have twoSelf := divides.reflexive at 2
      have scaled := divides.mul_mono with twoSelf, cDivides
      rw ← expanded at scaled
      have below := divides.le with scaled, succNonzero
      have characterization := divisors.membership_le at 2 ^ (k + 1), 2 * c
      rw shape
      exact characterization.backward with ⟨below, scaled⟩

-- σ(2^(k+1)) = 2 σ(2^k) + 1, hence σ(2^k) = 2^(k+1) − 1.
theorem sigma.two_pow_succ (k : Nat)
  shows sigma(2 ^ (k + 1)) = 2 * sigma(2 ^ k) + 1
proof
  have split := divisors.two_pow_succ at k
  have finite := divisors.finite at 2 ^ k
  have naturals := divisors.subset_naturals at 2 ^ k
  have finiteDoubles := doubles.finite with finite, naturals
  have weights := doubles.weights with naturals
  have doubled := natSum.doubles with finite, naturals
  rw ← sigma.definition at doubled
  have weightOne := identityWeight.natural at 1
  have valueOne := identityWeight.value at 1
  have fresh : ¬1 ∈ doubles(divisors(2 ^ k))
    intro inside
    have parts := doubles.membership.forward with inside
    obtain ⟨c, hc, shape⟩ := parts.right
    have cNat := divisors.natural with hc
    have even := Even.two_mul at c
    rw ← shape at even
    exact Even.not_one with even
  have inserted := natSum.insert with finiteDoubles, weights, weightOne, fresh
  have twice := mul.two_left at sigma(2 ^ k)
  rw [sigma.definition, split, inserted, doubled, valueOne, twice]
  rfl

theorem sigma.two_pow_add_one (k : Nat)
  shows sigma(2 ^ k) + 1 = 2 ^ (k + 1)
proof
  induction k
  | zero =>
    rw [pow.zero, sigma.one, pow.add_one, pow.zero]
    ring
  | succ(j) using ih =>
    have step := sigma.two_pow_succ at j
    have asSucc := add.right_one at j
    rw asSucc at step
    rw [step, ← asSucc, pow.add_one, ← ih]
    ring

theorem sigma.two_pow (k : Nat)
  shows sigma(2 ^ k) = 2 ^ (k + 1) - 1
proof
  have sum := sigma.two_pow_add_one at k
  symm
  exact sub.of_add_eq with sum

-- A natural number is perfect when it is the sum of its proper divisors:
-- σ(n) = 2n.
predicate PerfectNumber (n : Nat) := 0 < n ∧ sigma(n) = 2 * n

-- Euclid: 2^(p−1) · (2^p − 1) is perfect when 2^p − 1 is prime.
theorem perfectNumber.of_mersenne (p : Nat)
  requires prime : Prime(2 ^ p - 1)
  shows PerfectNumber(2 ^ (p - 1) * (2 ^ p - 1))
proof
  have qNonzero := Prime.nonzero with prime
  have pNonzero : p ≠ 0
    intro zero
    rw [zero, pow.zero, sub.self] at qNonzero
    apply qNonzero
    rfl
  have oneLeP : 1 ≤ p
    have characterization := LessEqual.one_iff_ne_zero at p
    exact characterization.backward with pNonzero
  have restored := sub.add_cancel at p, 1 with oneLeP
  have sum := sigma.two_pow_add_one at p - 1
  rw restored at sum
  have sigmaValue := sub.of_add_eq with sum
  have expanded := pow.add_one at 2, p - 1
  rw restored at expanded
  have oneLeTwo : 1 ≤ 2 by arith
  have powGeOne := pow.ge_one at 2, p with oneLeTwo
  have powRestored := sub.add_cancel at 2 ^ p, 1 with powGeOne
  have notTwo : 2 ^ p - 1 ≠ 2
    intro two
    have even := Even.two_mul at 2 ^ (p - 1)
    rw [← expanded, ← powRestored, two] at even
    have characterization := Even.iff_two_mul at 2 + 1
    obtain ⟨c, cNat, shape⟩ := characterization.forward with even
    have cBelow : c < 2 by arith
    have alternatives := natural.below_two with cBelow
    cases alternatives with zero | one
    . rw zero at shape
      arith
    . rw one at shape
      arith
  have odd := Prime.odd_of_ne_two with prime, notTwo
  have oddCharacterization := divides.not_two_iff_odd at 2 ^ p - 1
  have notDivides := oddCharacterization.backward with odd
  have two := Prime.two
  have coprimeTwo := Prime.gcd_of_not_divides with two, notDivides
  have coprime := coprime.pow_left at 2, 2 ^ p - 1, p - 1 with coprimeTwo
  have multiplicative := sigma.multiplicative with coprime
  have sigmaPrime := sigma.prime with prime
  have twoPositive : 0 < 2 by arith
  have powPositive := pow.positive at 2, p - 1 with twoPositive
  have powNonzero := Less.ne_zero_of_positive with powPositive
  have productNonzero := mul.nonzero with powNonzero, qNonzero
  have positive := Less.positive_of_ne_zero with productNonzero
  refine ⟨positive, _⟩
  rw [multiplicative, sigmaPrime, powRestored, sigmaValue, expanded]
  ring

-- Every positive natural is a power of two times an odd number.
theorem natural.two_pow_odd_below (B : Nat)
  shows ∀ n, Nat(n) → 0 < n → n ≤ B →
    ∃ k ∈ ℕ, ∃ m ∈ ℕ, Odd(m) ∧ n = 2 ^ k * m
proof
  induction B
  | zero =>
    intros n nNat positive below
    exfalso
    arith
  | succ(b) using ih =>
    intros n nNat positive below
    have alternatives := LessEqual.cases at n, succ(b) with below
    cases alternatives with strict | same
    . have belowB : n ≤ b by arith
      exact ih with n, nNat, positive, belowB
    . have parity := parity.exhaustive at n
      cases parity with even | odd
      . have characterization := Even.iff_two_mul at n
        obtain ⟨c, cNat, shape⟩ := characterization.forward with even
        have cPositive : 0 < c by arith
        have cBelow : c ≤ b by arith
        obtain ⟨k, kIn, m, mIn, oddM, form⟩ :=
          ih with c, cNat, cPositive, cBelow
        have kNat := nat.omega.backward with kIn
        have mNat := nat.omega.backward with mIn
        have succNat : Nat(k + 1) by typecheck
        have succIn := nat.omega.forward with succNat
        have expanded := pow.add_one at 2, k
        have shapeN : n = 2 ^ (k + 1) * m
          rw [shape, form, expanded]
          ring
        exact ⟨k + 1, succIn, m, mIn, oddM, shapeN⟩
      . have zeroNat : Nat(0) by typecheck
        have zeroIn := nat.omega.forward with zeroNat
        have nIn := nat.omega.forward with nNat
        have shapeN : n = 2 ^ 0 * n
          rw [pow.zero, mul.left_unit]
          rfl
        exact ⟨0, zeroIn, n, nIn, odd, shapeN⟩

theorem natural.two_pow_odd (n : Nat)
  requires positive : 0 < n
  shows ∃ k ∈ ℕ, ∃ m ∈ ℕ, Odd(m) ∧ n = 2 ^ k * m
proof
  have nNat : Nat(n) by typecheck
  have below := LessEqual.reflexive at n
  have all := natural.two_pow_odd_below at n
  exact all with n, nNat, positive, below

-- Euler: an even perfect number is 2^(p−1) · (2^p − 1) with 2^p − 1
-- prime.
theorem perfectNumber.even_form (n : Nat)
  requires perfect : PerfectNumber(n)
  requires even : Even(n)
  shows ∃ p ∈ ℕ, Prime(2 ^ p - 1) ∧ n = 2 ^ (p - 1) * (2 ^ p - 1)
proof
  rcases perfect with ⟨positive, sum⟩
  obtain ⟨k, kIn, m, mIn, odd, form⟩ := natural.two_pow_odd with positive
  have kNat := nat.omega.backward with kIn
  have mNat := nat.omega.backward with mIn
  have mNonzero : m ≠ 0
    intro zero
    rw [zero, mul.right_zero] at form
    rw form at positive
    exact Less.irreflexive at 0 with positive
  have kNonzero : k ≠ 0
    intro zero
    rw [zero, pow.zero, mul.left_unit] at form
    rw form at even
    have notEven := Odd.not_even with odd
    exact notEven with even
  have two := Prime.two
  have oddCharacterization := divides.not_two_iff_odd at m
  have notTwo := oddCharacterization.backward with odd
  have coprimeTwo := Prime.gcd_of_not_divides with two, notTwo
  have coprime := coprime.pow_left at 2, m, k with coprimeTwo
  have multiplicative := sigma.multiplicative with coprime
  have powSum := sigma.two_pow_add_one at k
  have expanded := pow.add_one at 2, k
  rw [form, multiplicative] at sum
  -- With L = σ(2^k): L · σ(m) = (L + 1) · m.
  have key : sigma(2 ^ k) * sigma(m) = (sigma(2 ^ k) + 1) * m
    rw [powSum, expanded, sum]
    ring
  have kPositive := Less.positive_of_ne_zero with kNonzero
  have oneLeK : 1 ≤ k by arith
  have oneLeTwo : 1 ≤ 2 by arith
  have powMonotone := pow.exponent_monotone at 2, 1, k with oneLeTwo, oneLeK
  rw pow.one_right at powMonotone
  have threeLe : 3 ≤ sigma(2 ^ k)
    rw expanded at powSum
    arith
  have lPositive : 0 < sigma(2 ^ k) by arith
  have lNonzero := Less.ne_zero_of_positive with lPositive
  -- L divides m, since L is coprime to L + 1.
  have consecutive := gcd.consecutive at sigma(2 ^ k)
  have coprimeSucc := coprime.symmetric with consecutive
  have dividesProduct : sigma(2 ^ k) ∣ m * (sigma(2 ^ k) + 1)
    have commuted : m * (sigma(2 ^ k) + 1) = sigma(2 ^ k) * sigma(m)
      rw key
      ring
    rw commuted
    exact divides.mul_left_self at sigma(2 ^ k), sigma(m)
  have dividesM := coprime.cancel_right with coprimeSucc, dividesProduct
  obtain ⟨t, tNat, mForm⟩ := dividesM
  -- Then σ(m) = m + t, with t a divisor of m below m.
  have sigmaM : sigma(m) = m + t
    have scaled : sigma(2 ^ k) * sigma(m) = sigma(2 ^ k) * (m + t)
      rw [key, ← mForm]
      ring
    exact mul.left_cancel_positive with lNonzero, scaled
  have tNonzero : t ≠ 0
    intro zero
    rw [zero, mul.right_zero] at mForm
    apply mNonzero
    symm
    exact mForm
  have tPositive := Less.positive_of_ne_zero with tNonzero
  have tLess : t < m
    have scaled := LessEqual.mul_right at 3, sigma(2 ^ k), t with threeLe
    rw mForm at scaled
    arith
  have tNotSelf := Less.ne with tLess
  have tDivides : t ∣ m
    have commuted : t * sigma(2 ^ k) = m
      rw ← mForm
      ring
    exact divides.intro with commuted
  have tCases : t = 1 ∨ t ≠ 1 by tauto
  cases tCases with one | notOne
  . -- t = 1: m = L = 2^(k+1) − 1 and σ(m) = m + 1, so m is prime.
    rw [one, mul.right_unit] at mForm
    rw one at sigmaM
    have aboveOne : 1 < m by arith
    have primeM := Prime.of_sigma with aboveOne, sigmaM
    have succNat : Nat(k + 1) by typecheck
    have succIn := nat.omega.forward with succNat
    have predecessor := sub.add_one_right at k
    have mersenne : 2 ^ (k + 1) - 1 = m
      have value := sub.of_add_eq with powSum
      rw value
      exact mForm
    have primeMersenne : Prime(2 ^ (k + 1) - 1)
      rw mersenne
      exact primeM
    have shapeN : n = 2 ^ (k + 1 - 1) * (2 ^ (k + 1) - 1)
      rw [mersenne, predecessor]
      exact form
    exact ⟨k + 1, succIn, primeMersenne, shapeN⟩
  . -- t ≠ 1: t is a divisor strictly between, and σ(m) ≥ m + 1 + t.
    exfalso
    have member := innerDivisors.member with mNonzero, tDivides, notOne,
      tNotSelf
    have bound := innerDivisors.le_sum with member
    have aboveOne : 1 < m by arith
    have split := sigma.split with aboveOne
    rw sigmaM at split
    arith

-- The even perfect numbers are exactly the Euclid numbers.
theorem perfectNumber.even_iff (n : Nat)
  requires even : Even(n)
  shows PerfectNumber(n) ↔
    ∃ p ∈ ℕ, Prime(2 ^ p - 1) ∧ n = 2 ^ (p - 1) * (2 ^ p - 1)
proof
  constructor
  . intro perfect
    exact perfectNumber.even_form with perfect, even
  . intro form
    obtain ⟨p, pIn, prime, shape⟩ := form
    have pNat := nat.omega.backward with pIn
    rw shape
    exact perfectNumber.of_mersenne with prime