There are infinite primes

Prime.agda

module Prime where

open import Coinduction

open import Data.Empty
open import Data.Nat
open import Data.Nat.Properties
open import Data.Nat.Divisibility
open import Data.Fin hiding (pred; _+_; _<_; _≤_; compare)
open import Data.Fin.Props hiding (_≟_)
open import Data.Sum
open import Data.Product
open import Data.Stream

open import Function

open import Relation.Nullary
open import Relation.Binary
open import Relation.Binary.PropositionalEquality

open import Induction.WellFounded as WF
import Induction.Nat as NI

open import Algebra

open DecTotalOrder decTotalOrder using () renaming (refl to ≤-refl; trans to ≤-trans)
open CommutativeSemiring commutativeSemiring using (*-assoc; *-comm; *-identity; +-assoc; +-comm; +-identity)

-- We can negate a decision on a type A
not : {A : Set} → Dec A → Dec (¬ A)
not (yes p) = no (λ z → z p)
not (no ¬p) = yes ¬p

-- Decisions exactly the law of excluded middle for a specific type, which gives us double negation elimination on that type
DNE : {A : Set} → Dec A → ¬ ¬ A → A
DNE (yes p) f = p
DNE (no ¬p) f = ⊥-elim (f ¬p)

-- A decision procedure for an arbitrary predicate P
Decide : {A : Set} (P : A → Set) → Set
Decide P = ∀ i → Dec (P i)

module Finite where
  -- Either the decidable predicate is true over its entire (finite) domain, or there exists an input for which it is false.
  -- Thanks to Ulf Norell for this!
  LPO : ∀ {n} {P : Fin n → Set} (p? : Decide P) → (∀ i → P i) ⊎ ∃ (¬_ ∘ P)
  LPO {zero}  p? = inj₁ (λ ())
  LPO {suc n} p? with LPO (p? ∘ suc)
  LPO {suc n} p? | inj₁ ok  with p? zero
  LPO {suc n} p? | inj₁ ok | yes p = inj₁ (λ { zero → p ; (suc i) → ok i })
  LPO {suc n} p? | inj₁ ok | no ¬p = inj₂ (zero , ¬p)
  LPO {suc n} p? | inj₂ (i , ¬pi) = inj₂ (suc i , ¬pi)

  -- A simple negation of the above, using DNE defined earlier, because this is the one we're usually interested in
  ¬LPO : ∀ {n} {P : Fin n → Set} (p? : Decide P) → ∃ P ⊎ (∀ i → ¬ P i)
  ¬LPO p? with LPO (not ∘ p?)
  ¬LPO p? | inj₁ x = inj₂ x
  ¬LPO p? | inj₂ (n , pf) = inj₁ (n , DNE (p? n) pf)

-- We actually care about the LPO on naturals < k, instead of finite sets, so let's convert the above definition to work with them
LPO : ∀ n {P : ℕ → Set} (p? : ∀ i → i < n → Dec (P i)) → (∀ i → i < n → P i) ⊎ ∃ (λ i → i < n × ¬ P i)
LPO _ p? with Finite.LPO (λ i → p? _ (bounded i))
LPO _ {P} p? | inj₁ x = inj₁ (λ i i<n → subst P (toℕ-fromℕ≤ i<n) (x _))
LPO _ p? | inj₂ (i , ¬P[i]) = inj₂ (toℕ i , bounded i , ¬P[i])

-- As above
¬LPO : ∀ n {P : ℕ → Set} (p? : ∀ i → i < n → Dec (P i)) → ∃ (λ i → i < n × P i) ⊎ (∀ i → i < n → ¬ P i) 
¬LPO n p? with LPO n (λ i i<n → not (p? i i<n))
¬LPO n p? | inj₁ x = inj₂ x
¬LPO n p? | inj₂ (i , i<n , pf) = inj₁ (i , i<n , DNE (p? i i<n) pf)

-- A couple of simple lemmas for later
unsplit-≤ : ∀ {n p} → p ≢ suc n → p ≤ suc n → p ≤ n
unsplit-≤ pf z≤n = z≤n
unsplit-≤ {zero} pf (s≤s z≤n) = ⊥-elim (pf refl)
unsplit-≤ {suc n} pf (s≤s m≤n) = s≤s (unsplit-≤ (λ x → pf (cong suc x)) m≤n)

*-cong′ : ∀ {x y} k → x ∣ y → x ∣ y * k
*-cong′ {x} n (divides q refl) rewrite *-comm q x | *-assoc x q n | *-comm x (q * n) = ∣-* (q * n)

-- I have a couple of other definitions of Prime lying around but this proved easier to deal with
record Prime (p : ℕ) : Set where
  constructor prime
  field
    -- We don't like primes smaller than 2
    p>1 : p > 1
    -- If you have something that divides our p, it's either 1 or p
    pf  : ∀ {d} → d ∣ p → d ≡ 1 ⊎ d ≡ p

-- A composite number has a prime divisor smaller than it
record Composite (c : ℕ) : Set where
  constructor composite
  field
    {d}  : ℕ
    P[d] : Prime d
    d<c  : d < c
    pf   : d ∣ c

-- Numbers aren't both prime and composite
¬Prime×Composite : ∀ {p} → Prime p → Composite p → ⊥
¬Prime×Composite (prime p>1 pf₁) (composite d>1 d<n pf₂) with pf₁ pf₂
¬Prime×Composite (prime p>1 pf₁) (composite (prime (s≤s ()) _) d<n pf₂) | inj₁ refl
¬Prime×Composite (prime p>1 pf₁) (composite d>1 d<n pf₂) | inj₂ refl = 1+n≰n d<n

-- This well-founded recursion is needed to appease the termination checker for the definition of primality below.
-- For some reason, the standard library includes a proof that _<′_ is well-founded but leaves out _<_ (probably because
-- it's so easy to do what I'm doing here). This incantation gives me well-founded recursion over _<_.
open Subrelation {_<₁_ = λ i j → i < j} ≤⇒≤′
open WF.All (well-founded NI.<-well-founded) renaming (wfRec to <-rec)

-- The meat of this module: All numbers above 1 are either prime or composite.
--
-- This decision procedure proceeds by using the ¬LPO machinery above to find a number between 2 and p - 1 that divides
-- p. If one exists, then the number is composite, and if it doesn't, the number is prime.
primality : ∀ p → p > 1 → Prime p ⊎ Composite p
primality 0 ()
primality 1 (s≤s ())
primality (suc (suc p)) p>1 = <-rec _ helper p
  where
  module ∣-Poset = Poset poset
  helper : ∀ p → (∀ i → i < p → Prime (2 + i) ⊎ Composite (2 + i)) → Prime (2 + p) ⊎ Composite (2 + p)
  helper p f with ¬LPO p (λ i i<p → (2 + i) ∣? (2 + p))
  -- There exists a number i, smaller than p, that divides p. This is the tricky case.
  -- Most of the trickiness comes from the fact that in the composite case, the ¬LPO gives us a number i smaller than p that
  -- divides p, but says nothing about whether that number is prime or not. Obviously, we'd like to ask whether that number
  -- is itself prime or not, but unfortunately it's not obviously structurally smaller than p. If we just naively call
  -- primality on it, Agda will complain about non-obvious termination. This is where the well-founded recursion comes
  -- into play. We know that i < p, and we showed above that _<_ is well founded, so we use the well-founded recursion voodoo
  -- and call ourselves recursively (f is effectively the primality test hidden behind the safe well-founded recursion stuff).
  helper p f | inj₁ (i , i<p , 2+i∣2+p) with f i i<p
  -- Aha! i is prime, so we're good and can return that proof in our composite proof
  helper p f | inj₁ (i , i<p , 2+i∣2+p) | inj₁ P[i] = inj₂ (composite P[i] (s≤s (s≤s i<p)) 2+i∣2+p)
  -- i is composite, but that means it has a prime divisor, which we can show transitively divides p! sweet!
  helper p f | inj₁ (i , i<p , 2+i∣2+p) | inj₂ (composite P[d] d<c pf) = inj₂ (composite P[d] (<-trans d<c (s≤s (s≤s i<p))) (∣-Poset.trans pf 2+i∣2+p))

  -- The boring case: there exist no numbers between 2 and p - 1 that divide p, so transmogrify that proof into our primality
  -- statement.
  helper p f | inj₂ y = inj₁ (prime (s≤s (s≤s z≤n)) lemma₁)
    where
    lemma₀ : ∀ {n p} → 2 + n ≤ 2 + p → 1 + n ≤ 1 + p
    lemma₀ (s≤s 1+n≤1+p) = 1+n≤1+p

    lemma₁ : {d : ℕ} → d ∣ 2 + p → d ≡ 1 ⊎ d ≡ 2 + p
    lemma₁ {d} d∣2+p with d ≟ 1 | d ≟ 2 + p
    lemma₁ d∣2+p | yes refl | yes ()
    lemma₁ d∣2+p | yes d≡1 | no  d≢2+p = inj₁ d≡1
    lemma₁ d∣2+p | no  d≢1 | yes d≡2+p = inj₂ d≡2+p
    lemma₁ {zero} d∣2+p | no d≢1 | no d≢2+p rewrite 0∣⇒≡0 d∣2+p = ⊥-elim (d≢2+p refl)
    lemma₁ {suc zero} d∣2+p | no d≢1 | no d≢2+p = ⊥-elim (d≢1 refl)
    lemma₁ {suc (suc n)} d∣2+p | no d≢1 | no d≢2+p = ⊥-elim (y n (unsplit-≤ (d≢2+p ∘ cong suc) (lemma₀ (∣⇒≤ d∣2+p))) d∣2+p)

-- The factorial function!
_! : ℕ → ℕ
zero ! = 1
suc n ! = suc n * n !

-- All positive numbers less than n divide n!
!-divides : ∀ n {m} → suc m ≤ n → suc m ∣ n !
!-divides zero ()
!-divides (suc n) {m} m≤n rewrite *-comm (suc n) (n !) with suc m ≟ suc n
!-divides (suc n) m≤n | yes refl = ∣-* (n !)
!-divides (suc n) m≤n | no ¬p = *-cong′ (suc n) (!-divides n (unsplit-≤ ¬p m≤n))

!>0 : ∀ n → n ! > 0
!>0 zero = s≤s z≤n
!>0 (suc n) = _*-mono_ {1} {suc n} (s≤s z≤n) (!>0 n)

n≤n! : ∀ n → n ≤ n !
n≤n! zero = z≤n
n≤n! (suc n) rewrite sym (proj₂ *-identity (suc n)) = _*-mono_ {suc n} ≤-refl (!>0 n)

∣-∸′ : ∀ {i n} m → i ∣ m + n → i ∣ n → i ∣ m
∣-∸′ {i} {n} m pf₁ pf₂ rewrite +-comm m n = ∣-∸ pf₁ pf₂

¬≤⇒> : ∀ {x y} → ¬ x ≤ y → y < x
¬≤⇒> {x} {y} pf with compare x y
¬≤⇒> pf | less x k = ⊥-elim (pf (≤-step (m≤m+n x k)))
¬≤⇒> pf | equal x = ⊥-elim (pf ≤-refl)
¬≤⇒> pf | greater x k = s≤s (m≤m+n x k)

generate : ∀ n → ∃ λ i → i > n × Prime i
generate n with primality (1 + n !) (s≤s (!>0 n))
generate n | inj₁ x = , (s≤s (n≤n! n) , x)
generate n | inj₂ (composite {d} P[d] d<c pf) with suc n ≤? d
generate n | inj₂ (composite P[d] d<c pf) | yes p = , (p , P[d])
generate n | inj₂ (composite {0} (prime () _) d<c pf) | no ¬p
generate n | inj₂ (composite {suc d} P[d] d<c pf) | no ¬p with ∣1⇒≡1 (∣-∸′ 1 pf (!-divides n (¬≤⇒> (¬p ∘ s≤s))))
generate n | inj₂ (composite {suc .0} (prime (s≤s ()) _) d<c pf) | no ¬p | refl

data Primes (p : ℕ) : Set where
  _[_]∷_ : ∀ {p′} (P[p] : Prime p) (p<p′ : p < p′) (Ps : ∞ (Primes p′)) → Primes p

toStream : ∀ {p} → Primes p → Stream ℕ
toStream {p} (P[p] [ p<p′ ]∷ Ps) = p ∷ ♯ toStream (♭ Ps)

primes′ : ∀ {p} → Prime p → Primes p
primes′ {p} P[p] with generate p
primes′ {p} P[p] | p′ , p′>p , P[p′] = P[p] [ p′>p ]∷ ♯ (primes′ P[p′])

P[2] : Prime 2
P[2] with primality 2 (s≤s (s≤s z≤n))
P[2] | inj₁ x = x
P[2] | inj₂ (composite {0} (prime () _) (s≤s m≤n) pf)
P[2] | inj₂ (composite {1} (prime (s≤s ()) _) (s≤s m≤n) pf)
P[2] | inj₂ (composite {suc (suc n)} P[d] (s≤s (s≤s ())) pf)

primes : Primes 2
primes = primes′ P[2]

I don't understand the language, but this is almost like reading a math book
LPO : ∀ n {P : ℕ → Set} (p? : ∀ i → i < n → Dec (P i)) → (∀ i → i < n → P i) ⊎ ∃ (λ i → i < n × ¬ P i)
very cool