-- Problem 35: Determine the prime factors of a given positive integer.
-- Trial division will do here, I didn't want to go into implementing a sieve if I didn't
-- have to (joke's on me). If we considered candidate factors in an ascending order, we
-- wouldn't need to make sure they're prime:
primefactors :: Integer -> [Integer]
primefactors n = primefactors' n 2
where
primefactors' 1 _ = []
primefactors' n f
-- Does f divide n? Good, try n/f now:
| n `mod` f == 0 = f : primefactors' (n `div` f) f
-- If not, just move on to the next number:
| otherwise = primefactors' n (f + 1)
-- Problem 36: Determine the prime factors of a given positive integer.
-- To encode the multiplicity of prime factors, we use our knowledge of `group` from the
-- time of Problem 10:
prime_factors_group :: Integer -> [(Integer, Int)]
prime_factors_group = map encode . group . primefactors
where encode xs = (head xs, length xs)
-- Problem 37: Calculate Euler's totient function phi(m).
-- Since the formula is given to us, we see that it's a simple application of foldl (aka
-- reduce) on (p-1)*p^(m-1) (mapped onto each prime):
phi :: Integer -> Integer
phi n = foldl (*) 1 [(p - 1) * p ^ (m - 1) | (p, m) <- prime_factors_group n]
Problem 39 is a handful. Here I was thinking I wouldn’t have to implement a prime sieve, but the problem asks for one. I decided to go for an infinite one, then drop all unnecessary primes:
primesR a b = takeWhile (<= b) $ dropWhile (< a) primes
But how will we implement a sieve? We will build one based on the idea of a prime wheel. Read on it. We’ll construct a pretty bigger wheel than the P6 candidate prime generator we used earlier. First, we’ll define a wheel by its circumference and the “spiked” positions:
data Wheel = Wheel Integer [Integer]
-- We generate our candidates K*n+R via an infinite list:
roll (Wheel n rs) = [n*k+r | k <- [0..], r <- rs]
-- This is the unit wheel, which yields all numbers:
w0 = Wheel 1 [1]
-- We can produce any wheel from w0 and a "spike" p by excluding multiples of p:
nextSize (Wheel n rs) p = Wheel (p*n) [r2 | k <- [0..(p-1)]
, r <- rs
, let r2 = n*k+r, r2 `mod` p /= 0]
-- We can use nextSize on any list of numbers to "avoid". Multiples of 2 would be the
-- simplest example, and we can extend that to however many we want.
mkWheel ds = foldl nextSize w0 ds
-- Now what do we do with the wheel? Starting with a list of primes (let's call it "small")
-- we generate a wheel, roll it and start testing for primality...but, against what? Well,
-- the "small" list will do! (remember prime factorisation?) However, once a number is
-- proven prime, we can add it to this list and use it for future candidates.
primes = small ++ large
where 1:p:candidates = roll $ mkWheel small
small = [2,3,5,7]
large = p : filter isPrime candidates
isPrime n = all (not . divides n)
$ takeWhile (\p -> p*p <= n) large
divides n p = n `mod` p == 0
This sieve took me a while to get working correctly. The most difficult part was
constructing the wheel in a generic fashion (from the unit wheel), and it was made easier
only because I’d read a few things in the past about wheel factorisation1. Notice that it
highlights a particularity of big wheels: they’re not good generators of primes. Since the
density of primes decreases (prime number theorem) with increasing N, we remove less and
less of the composites as we increase the wheel size. Notice that the wheel generated from
[2,3] removes 4/6 = 0.66 of the composites. The wheel generated from [2,3,5,7] removes
162/210 = 0.7714. To remove 90% of the composites, we must use the primes up to 251 (54
in number), and so forth.
-
and I’m still not sure I can properly explain it in terms of map/reduce. ↩