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.
Haskell
import qualified Data.Set as Set
minus :: [Int] -> [Int] -> [Int]
@(x:xt) ys@(y:yt) = case compare x y of
minus xsLT -> x : minus xt ys
EQ -> minus xt yt
GT -> minus xs yt
= a
minus a _
union :: [Int] -> [Int] -> [Int]
@(x:xt) ys@(y:yt) = case compare x y of
union xsLT -> x : union xt ys
EQ -> x : union xt yt
GT -> y : union xs yt
= a
union a [] = b
union [] b
uniq :: Ord a => [a] -> [a]
= uniq' Set.empty xs where
uniq xs = []
uniq' _ [] :ys) | Set.member y set = uniq' set ys
uniq' set (y| otherwise = y : uniq' (Set.insert y set) xs
primes :: [Int]
= 2:3:5:7: gaps 11 wheel (fold3t $ roll 11 wheel primes_)
primes where
= 11: gaps 13 (tail wheel) (fold3t $ roll 11 wheel primes_) -- separate feed
primes_ :xs): ~(ys:zs:t)) = x : union xs (union ys zs)
fold3t ((x`union` fold3t (pairs t) -- fold3t: 5% ~ 10% speedup
:xs):ys:t) = (x : union xs ys) : pairs t
pairs ((x= 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:
wheel 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
@(w:t) cs@ ~(c:u) | k==c = gaps (k+w) t u -- (* better fold, w/ Wheel! *)
gaps k ws| True = k : gaps (k+w) t cs
@(w:t) ps@ ~(p:u) | k==p = scanl (\c d->c+p*d) (p*p) ws
roll k ws: roll (k+w) t u
| True = roll (k+w) t ps
factorize :: Int -> [Int]
= primeFactors n primes where
factorize n 1 _ = []
primeFactors = []
primeFactors _ [] :ps) | m < p * p = [m]
primeFactors m (p| r == 0 = p : primeFactors q (p:ps)
| otherwise = primeFactors m ps
where (q, r) = quotRem m p
totient :: Int -> Double
1 = 1.0
totient = (fromIntegral n) * product [1.0 - (1.0 / (fromIntegral p)) | p <- uniq $ factorize n]
totient n
main :: IO ()
= print $ snd $ maximum [((fromIntegral n) / (totient n), n) | n <- [1..1000000]] main
$ 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:
2)
res.append(//= 2
n # try odd numbers up to sqrt(n)
= math.sqrt(n+1)
limit = 3
i while i <= limit:
if n % i == 0:
res.append(i)//= i
n = math.sqrt(n+i)
limit else:
+= 2
i 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 .default = 0 # return 0 instead nil if key not found in hash
factors= orig
n = 2
i = 4 # square of i
sqi <= n do
while sqi .modulo(i) == 0 do
while n/= i
n [i] += 1
factorsend
# we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
+= 2 * i + 1
sqi += 1
i end
if (n != 1) && (n != orig)
[n] += 1
factorsend
factorsend
def totient(n)
* factorize(n).each_key.map { |p| 1 - 1.0 / p }.reduce(1, :*)
n 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