99 problems in Haskell, 35-39

2015 Jul 29

-- 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
    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.

  1. and I’m still not sure I can properly explain it in terms of map/reduce. 

« Past Future »


comments powered by Disqus