## Question

Euler’s Totient function, φ(n) [sometimes called the phi function], is used to determine the number of numbers less than n which are relatively prime to n. For example, as 1, 2, 4, 5, 7, and 8, are all less than nine and relatively prime to nine, φ(9)=6.

 n Relatively Prime φ(n) n/φ(n) 2 1 1 2 3 1,2 2 1.5 4 1,3 2 2 5 1,2,3,4 4 1.25 6 1,5 2 3 7 1,2,3,4,5,6 6 1.1666… 8 1,3,5,7 4 2 9 1,2,4,5,7,8 6 1.5 10 1,3,7,9 4 2.5

It can be seen that n=6 produces a maximum n/φ(n) for n ≤ 10.

Find the value of n ≤ 1,000,000 for which n/φ(n) is a maximum.

import qualified Data.Set as Set

minus :: [Int] -> [Int] -> [Int]
minus xs@(x:xt) ys@(y:yt) = case compare x y of
LT -> x : minus xt ys
EQ ->     minus xt yt
GT ->     minus xs yt
minus a         _         = a

union :: [Int] -> [Int] -> [Int]
union xs@(x:xt) ys@(y:yt) = case compare x y of
LT -> x : union xt ys
EQ -> x : union xt yt
GT -> y : union xs yt
union a         []        = a
union []        b         = b

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

primes :: [Int]
primes = 2:3:5:7: gaps 11 wheel (fold3t $roll 11 wheel primes_) where primes_ = 11: gaps 13 (tail wheel) (fold3t$ roll 11 wheel primes_)     -- separate feed
fold3t ((x:xs): ~(ys:zs:t)) = x : union xs (union ys zs)
union fold3t (pairs t)              -- fold3t: 5% ~ 10% speedup
pairs ((x:xs):ys:t)         = (x : union xs ys) : pairs t
wheel = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:
4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel
gaps k ws@(w:t) cs@ ~(c:u) | k==c  = gaps (k+w) t u              -- (*  better fold, w/ Wheel!   *)
| True  = k : gaps (k+w) t cs
roll k ws@(w:t) ps@ ~(p:u) | k==p  = scanl (\c d->c+p*d) (p*p) ws
: roll (k+w) t u
| True  = roll (k+w) t ps

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

totient :: Int -> Double
totient 1 = 1.0
totient n = (fromIntegral n) * product [1.0 - (1.0 / (fromIntegral p)) | p <- uniq $factorize n] main :: IO () main = print$ snd $maximum [((fromIntegral n) / (totient n), n) | n <- [1..1000000]] $ ghc -O2 -o totient totient.hs
$time ./totient real 0m0.855s user 0m0.848s sys 0m0.004s ## Python #!/usr/bin/env python import math from operator import mul from functools import reduce def factorize(n): res = [] # iterate over all even numbers first. while n % 2 == 0: res.append(2) n //= 2 # try odd numbers up to sqrt(n) limit = math.sqrt(n+1) i = 3 while i <= limit: if n % i == 0: res.append(i) n //= i limit = math.sqrt(n+i) else: i += 2 if n != 1: res.append(n) return res def totient(n): return n * reduce(mul, [1 - (1.0 / p) for p in set(factorize(n))]) def main(): print(max((n / totient(n), n) for n in range(2, 1000001))[1]) if __name__ == "__main__": main() $ time python3 totient.py
real   0m20.170s
user   0m20.036s
sys    0m0.004s

## Ruby

#!/usr/bin/env ruby
require 'set'

def factorize(orig)
factors = {}
n = orig
i = 2
sqi = 4 # square of i
while sqi <= n do
while n.modulo(i) == 0 do
n /= i
factors[i] += 1
end
# we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
sqi += 2 * i + 1
i += 1
end

if (n != 1) && (n != orig)
factors[n] += 1
end
factors
end

def totient(n)
n * factorize(n).each_key.map { |p| 1 - 1.0 / p }.reduce(1, :*)
end

puts (2..1000000).map { |n| [n / totient(n), n] }.max[1]
\$ time ruby totient.rb
real   0m17.427s
user   0m17.284s
sys    0m0.036s