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.

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.

Haskell

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.797s
user   0m0.796s
sys    0m0.000s

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   0m13.748s
user   0m13.731s
sys    0m0.016s

Ruby

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

def factorize(orig)
  factors = {}
  factors.default = 0 # return 0 instead nil if key not found in hash
  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   0m16.543s
user   0m16.503s
sys    0m0.040s