I'm working on problem #69 from the Project Euler website. It regards Euler's Totient function. The definition of this function can be found here. Anyway, the problem asks to find the maximum value of n / Phi(n) , where Phi is the Totient function, over the integers 2 to 1,000,000. I know that my code works, as it correctly finds n=6 as a maximum when the search interval is 2 to 10. But, it is horrendously slow - using an upper limit of 1,000 extends the computation time to roughly 30 seconds, meaning that using an upper limit of 1,000,000 would take roughly at least 8 hours! On the Project Euler website it states that no single computation should take more than a minute on a moderately powered computer. My computer is plenty powerful so there is no deficiency in that regard. Perhaps the Jupyter Notebooks IDE I'm using has a particularly inefficient compiler? I'm not really sure. Help on this matter would be greatly appreciated. Below is my code:
def FactorsOf(n,include_ends=False): #Returns the factors of a positive integer in an array.
factors = [] #The include_ends param can be set to True if 1 & n are to be
a=0 #included in the returned array.
b=0
if (include_ends):
a=1
b=n+1
else:
a=2
b=n
for k in range(a,b):
if (n%k==0):
factors.append(k)
return factors
def AreRelativelyPrime(a,b):
a_factors=FactorsOf(a,include_ends=True)
b_factors=FactorsOf(b,include_ends=True)
for i in range(1,len(a_factors)): #Searches through both factor arrays to see if there
for j in range(1,len(b_factors)): #are any elements in common. Of course the first element,
if (a_factors[i]==b_factors[j]): # 1, is excluded.
return False
return True
def Totient(n):
totient=1 #The Totient function's minimum value is 1.
n_factors = FactorsOf(n)
for m in range(2,n):
if(AreRelativelyPrime(n,m)): #Increments the Totient function every time we find a
totient+=1 # number relatively prime to n.
return totient
n_over_phi_MAX = 2
maxes = []
for n in range(2,1001):
n_over_phi = n/Totient(n)
if(n_over_phi > n_over_phi_MAX):
n_over_phi_MAX = n_over_phi
maxes.append(n)
print("The maxiumum value of n/phi(n) is " + str(n_over_phi_MAX) + " at a value of " + str(maxes[-1]))
Factorizing every integer from 2 to 1M will take a long time in any programming language no matter how you program it. You need to use a different algorithm altogether. A more efficient approach could be to leverage the fact that the phi function is multiplicative. Given that phi of the power of a prime is easy to compute: p**k - p**(k-1), you can generate primes (using the sieve of Eratosthenes) and just run through all multiples of the powers of these primes. Any multiplier of these prime powers (that is not itself a multiple of said prime) will have a gcd of 1 with the prime power. The multiplicative property of the phi function then allows computation of the phi value for that multiple based on previously calculated phi values (of the multiplier and of the prime power).
This runs in 0.5 second on my laptop.
Here's an example of the algorithm:
N = 1000000
phi = [0]+[1]*N
done = [False]*(N+1)
prime = [True]*(N+1)
for p in range(2,N+1):
if not prime[p]: continue
prime[p*p::p] = [False]*len(range(p*p,N+1,p)) # sieve of Eratosthenes
n = p
while n < N: # n is every power of prime p (in range of N)
phi[n] = n - n//p # phi(n) for a power of a prime
done[n] = True
for m in range(2,N//n+1): # Multipliers of n will have gcd(m,n) == 1
if m%p == 0: continue # if m not divisible by p
if not done[m]: continue # Expand from known phi(m)
nm = n*m
phi[nm] = phi[n]*phi[m] # totient function is multiplicative
done[nm] = True
n *= p
# once you have the phi values for all numbers in range,
# you can get the maximum ratio of n over Phi(n)
maxN2Phi = max(range(2,N),key=lambda n:n/phi[n])
output:
# print(f"Max n/phi(n) from 1 to {N}: {maxN2Phi} / {phi[maxN2Phi]} ({maxN2Phi/phi[maxN2Phi]})")
print("\nFIRST 143 phi(n):")
for b in range(0,144,12):
print(" ".join(f"{n or '':4}" for n in phi[b:b+12]))
FIRST 143 phi(n):
1 1 2 2 4 2 6 4 6 4 10
4 12 6 8 8 16 6 18 8 12 10 22
8 20 12 18 12 28 8 30 16 20 16 24
12 36 18 24 16 40 12 42 20 24 22 46
16 42 20 32 24 52 18 40 24 36 28 58
16 60 30 36 32 48 20 66 32 44 24 70
24 72 36 40 36 60 24 78 32 54 40 82
24 64 42 56 40 88 24 72 44 60 46 72
32 96 42 60 40 100 32 102 48 48 52 106
36 108 40 72 48 112 36 88 56 72 58 96
32 110 60 80 60 100 36 126 64 84 48 130
40 108 66 72 64 136 44 138 48 92 70 120
Not printing the answer here in accordance with project Euler's rules
Another alternative is to leverage the phi(n) = n * (1-1/p) * (1-1/p) * ... property. This is much simpler to implement but runs slower in my tests (still under 1 second though):
N = 1000000
phi = list(range(N+1))
prime = [True]*(N+1)
for p in range(2,N+1):
if not prime[p]: continue
prime[p*p::p] = [False]*len(range(p*p,N+1,p)) # sieve of Eratosthenes
phi[p::p] = [ n - n//p for n in phi[p::p] ] # phi(n) = n*(1-1/p)*(1-1/p)...
maxN2Phi = max(range(2,N),key=lambda n:n/phi[n])
[EDIT] A much much faster solution
By focusing closer to the objective (i.e. n/phi(n)) rather than producing the phi values, we can strategize that the largest ratio will be one where n is as large as possible and phi(n) is as small as possible.
Given that phi(n) = n*(1-1/p)*(1-1/p)... each prime in n reduces the value of phi(n). Moreover, the smaller primes reduce it more than the larger ones. So, we need to have the greatest possible number of the smallest possible primes as factors of n. This points to selecting all the first primes until their product becomes greater than 1,000,000.
Also, since we want n to be as large as possible, once we reach the maximum number of primes that will fit, we can further multiply the product of these primes by 2 or 3, or 4 ... as long as n remains below 1,000,000.
This approach gives the solution directly with only a prime number generation up to a very small fraction of 1,000,000.
The algorithm produces the solution in 0.00005 second (50 microseconds)
here's the code:
N = 1000000
prime = [True]*int(N**0.5) # largest prime used will be smaller than square root of N
n = 1
for p in range(2,len(prime)):
if not prime[p]: continue
prime[p*p::p] = [False]*len(prime[p*p::p]) # Eratosthenes
if n*p < N: n *= p # product of first primes
else: break # while product fits within N
n = n*(N//n) # multiply to maximize n within N
phi = n # compute phi(n)
for f in range(2,p):
if prime[f]: phi -= phi//f # n*(1-1/p)(1-1/p) ...
if printIt:
print(f"Max n/phi(n) from 1 to {N}: n={n} phi(n)={phi} ({n/phi})")