Project Euler Problem 243 Solution

Question

A positive fraction whose numerator is less than its denominator is called a proper fraction.

For any denominator, d, there will be d - 1 proper fractions; for example, with d = 12:

\frac{1}{12}, \frac{2}{12}, \frac{3}{12}, \frac{4}{12}, \frac{5}{12}, \frac{6}{12}, \frac{7}{12}, \frac{8}{12}, \frac{9}{12}, \frac{10}{12}, \frac{11}{12}.

We shall call a fraction that cannot be cancelled down a resilient fraction.

Furthermore we shall define the resilience of a denominator, R(d), to be the ratio of its proper fractions that are resilient; for example, R(12) = \frac{4}{11}.

In fact, d = 12 is the smallest denominator having a resilience R(d) < \frac{4}{10}.

Find the smallest denominator d, having a resilience R(d) < \frac{15499}{94744}.

Haskell

import Data.List (union)
import qualified Data.Set as Set

pairwise :: (a -> a -> a) -> [a] -> [a]
pairwise f (xs:ys:t) = f xs ys : pairwise f t
pairwise _ t = t

primes :: [Int]
primes = 2 : _Y ((3 :) . gaps 5 . _U . map (\p-> [p*p, p*p+2*p..]))
    where
        _Y g = g (_Y g)                      -- recursion, Y combinator
        _U ((x:xs):t) = x : (union xs . _U . pairwise union) t   -- ~= nub.sort.concat
        gaps k s@(x:xs) 
            | k < x     = k : gaps (k+2) s    -- ~= [k,k+2..]\\s, when
            | otherwise =     gaps (k+2) xs   --  k <= head s && null(s\\[k,k+2..])

factorize :: Int -> [Int]
factorize n = primeFactors n primes where
    primeFactors 1 _ = []
    primeFactors _ [] = []
    primeFactors m (p:ps) | m < p * p = [m]
                          | r == 0 = p : primeFactors q (p:ps)
                          | otherwise = primeFactors m ps
                          where (q, r) = quotRem m p

uniq :: Ord a => [a] -> [a]
uniq xs = uniq' Set.empty xs where
    uniq' _ [] = []
    uniq' set (y:ys) | Set.member y set = uniq' set ys
                     | otherwise = y : uniq' (Set.insert y set) xs

totient :: Int -> Double
totient 1 = 1.0
totient n = (fromIntegral n) * product [1.0 - (1.0 / (fromIntegral p)) | p <- uniq $ factorize n]

resilience :: Int -> Double
resilience d = (totient d) / (fromIntegral (d - 1))

primorials :: [Int]
primorials = scanl1 (*) primes

candidates :: [Int]
candidates = expand 1 primorials where
    expand m ps@(x:y:_) | m * x < y = m * x : expand (m+1) ps
                        | otherwise = expand 1 (tail ps)

main :: IO ()
main = print $ head [d | d <- candidates, resilience d < 15499 / 94744]

$ ghc -O2 -o resilience resilience.hs
$ time ./resilience
real   0m0.002s
user   0m0.000s
sys    0m0.002s