> module Primes(primes)
>  where
>  import SequenceUtilities(takeUntil)


sequence of all prime numbers: [2, 3, 5, 7, 11, 13, 17, ...]

>  primes :: [Integer]
>  primes = (map head . iterate sieve) [2 ..]


Eratosthenes sieve

>  sieve :: [Integer] -> [Integer]
>  sieve(p : xs) = filter ((/= 0) . (`mod` p)) xs


prime factorization of a given positive integer
Note: if n::Integer, n > 0, and fs = factors n, then
      product fs = n, all elements of fs are prime, and
      fs is monotonically nondecreasing

>  factors :: Integer -> [Integer]
>  factors n =
>    (filter(/= 1) . map thd3 .
>       takeUntil((== 1) . fst3) . iterate factor) (n, primes, 1)


factor a given integer into two pieces, using given divisor if possible
discard proposed factor from sequence if it isn't a factor of the integer
Note: if n,p::Integer, n > 0, p > 0,
        f has no factors x with 1 < x < p, and
        factor(n, [p]++ps,  _) = (m, qs, f),
      then
        n = m*f
        f = p and qs = [p]++ps if p divides n
        f = n if p*p > n (which means that no q with p < q < n divides n)
        f = 1, p doesn't divide n, and qs = ps, otherwise

>  factor :: (Integer, [Integer], Integer) -> (Integer, [Integer], Integer)
>  factor (n, p : ps, _)
>    | p*p > n          = (1, ps, n)
>    | remainder == 0   = (quotient, [p] ++ ps, p)
>    | otherwise        = (n, ps, 1)
>    where
>    (quotient, remainder) = divMod n p


>  fst3(a, b, c) = a
>  thd3(a, b, c) = c


Theorem E
If [p, x1, x2, ...]::[Integer]
  (1) p < x1 < x2 < ...
  (2) p is a prime number
  (3) no number x with 1 < x < p divides any number in [x1, x2, ...]
  (4) all prime numbers exceeding p are in [x1, x2, ...]
and
  sieve[p, x1, x2, ...] = [q, y1, y2, ...] :: [Integer]
then
  (1) q < y1 < y2 < ...
  (2) q is a prime number
  (3) no number y with 1 < y < q divides any number in [y1, y2, ...]
  (4) all primes exceeding q are in [y1, y2, ...]
  (5) there is no prime number x with p < x < q
  (6) [q, y1, y2, ...] is a subsequence of [p, x1, x2, ...]
  (7) p is not in [q, y1, y2, ...]

proof:

 sieve[p, x1, x2, ...] = [q, y1, y2, ...]                  -- hypothesis
                       = filter d [x1, x2, ...]        -- def'n of sieve
                         where d(x) is True iff p divides x
Therefore,
 [q, y1, y2, ...] is a subsequence of [x1, x2, ...] -- property of filter
This proves conclusions (6) and (7).
Also,
 q < y1 < y2 < ...                   -- hypothesis (1) and conclusion (6)
This proves conclusion (1).

Also, 
 p doesn't divide any number in [q, y1, y2, ...] -- property of filter(d)
                                             -- call this Observation (P)
So,
 no number x with 1 < x <= p divides any number in [q, y1, y2, ...]
                            -- Observation (P) along with hypothesis (2)
Furthermore, because all numbers x with x < q failed to pass filter(d),
all such numbers must be divisible by p. Therefore,
 no number x < q in [x1, x2, ...] divides any number in [q, y1, y2, ...]
  -- otherwise p (which divides all such numbers x) would divide a number
     in [q, y1, y2, ...], which would violate Observation (P)
So, no number x with 1 < x < q divides any number in [q, y1, y2, ...].
  -- if it did, it would have a prime factor t (x = t*v) with t < x < q;
     t would have to be p or greater because, by hypothesis (3),
     no number less than p divides any number in [x1, x2, ...];
     t could not be p, because then q wouldn't have passed filter(d);
     so, t > p, which means (by hypothesis 4), that t is in [x1, x2, ...];
     and, since t < q, t would have to be one of the numbers dropped by
     filter(d) from [x1, x2, ...], which can't be the case since filter(d)
     only drops multiples of p (i.e., non-prime numbers)
This proves conclusions (2) and (3).
It also proves conclusion (5) because any prime between p and q would have
to be one of the numbers dropped from [x1, x2, ...] by filter(d).

Conclusion (4) follows from hypothesis (1) (all primes exceeding p are in
[x1, x2, ...]) and the overvation that filter(d) delivers all primes in
the sequence supplied as its argument.
                                                                     QED

Theorem P
All primes numbers are in the sequence primes and no non-prime numbers are.

proof

Let s(1) = [2 ..] and s(n) = sieve(s(n-1)) for each n > 0
Then, by definition, primes = [head(s(1)), head(s(2)), head(s(3)), ...]

Let P(n) stand for the following statement:
      (1) s(n) is strictly increasing
  and (2) the first element of s(n) is a prime number
  and (3) no x with 1 < x < head(s(n)) divides any number in tail(s(n))
  and (4) all prime numbers exceeding head(s(n)) are in tail(s(n))
  and (5) there are no prime numbers x with head(s(n)) < x < head(s(n+1))
  and (6) s(n+1) is a subsequence of s(n)
  and (7) head(s(n)) is not in s(n+1)

Observe that parts (1)-(4) of P(1) are true and that parts (5)-(7) of P(1)
are therefore true by Theorem E.

Induction hypothesis: P(n) is true

Then, by Theorem E, P(n+1) follows from the induction hypothesis.

Therefore, by the induction principle, P(n) is true for all n >= 0.

So, head(s(n)) is prime, for all n >= 0            -- part (2) of P(n)
This proves that no non-primes are in the sequence primes.

Suppose some prime number t is not in the sequence primes.

Observe that primes is a subsequence of [2 ..] (because sieve
eliminates the first element of its argument).
So, primes is a strictly increasing sequence with an infinite number of
elements (infinite because iterate delivers an infinite sequence).
Therefore, there must be two adjacent elements p and q in primes
such that p < t < q. But, this would violate conclusion (5) of Theorem E.

So, the supposition that the sequence primes failed to include some
prime number t must have been wrong.
This proves that the seqeunce primes contains all prime numbers.
                                                                    QED