为什么这个 Python 代码这么慢？我怎样才能使它更有效率？

Question

我正在使用 Project Euler 网站上的 problem #69。它涉及欧拉的Totient函数。这个函数的定义可以在here找到。无论如何，问题要求在整数 2 到 1,000,000 上找到 n / Phi(n) 的最大值，其中 Phi 是 Totient 函数。我知道我的代码有效，因为当搜索间隔为 2 到 10 时，它正确地找到了 n=6 作为最大值。但是，它可怕慢 - 使用 1,000 的上限会延长计算时间到大约 30 秒，这意味着使用 1,000,000 的上限大约需要至少 8 小时！在 Project Euler 网站上，它指出在中等功率的计算机上任何一次计算都不应超过一分钟。我的电脑很强大，所以在这方面没有任何不足。也许我正在使用的 Jupyter Notebooks IDE 的编译器效率特别低？我不太确定。对此事的帮助将不胜感激。以下是我的代码：

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]))

Answer 1

在任何编程语言中，无论您如何编程，将每个整数从 2 分解到 1M 都需要很长时间。您需要完全使用不同的算法。更有效的方法可能是利用 phi 函数是乘法的这一事实。鉴于素数的幂的 phi 很容易计算：p**k - p**(k-1)，您可以生成素数（使用 Eratosthenes 的筛子）并且只需遍历这些素数的所有幂的倍数。这些素数幂的任何乘数（其本身不是所述素数的倍数）将与素数幂的 gcd 为 1。然后，phi 函数的乘法属性允许根据先前计算的 phi 值（乘数和素数幂）计算该倍数的 phi 值。

这在我的笔记本电脑上运行 0.5 秒。

这是一个算法示例：

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])

输出：

# 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

这里没有按照项目欧拉规则打印答案

另一种选择是利用 phi(n) = n * (1-1/p) * (1-1/p) * ... 属性。这实现起来要简单得多，但在我的测试中运行速度较慢（尽管仍不到 1 秒）：

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])

[编辑] 更快的解决方案

通过更接近目标（即 n/phi(n)）而不是产生 phi 值，我们可以制定策略，最大比率将是n 尽可能大而phi(n) 尽可能大的比率尽可能小。

假设 phi(n) = n*(1-1/p)*(1-1/p)...n 中的每个素数都会减小 phi(n) 的值。此外，较小的素数比较大的素数减少更多。因此，我们需要将尽可能多的最小可能素数作为 n 的因数。这指向选择所有第一个素数，直到它们的乘积大于 1,000,000。

此外，由于我们希望n 尽可能大，一旦我们达到适合的最大素数数量，我们可以进一步将这些素数的乘积乘以 2 或 3 或 4 ...只要n 保持在 1,000,000 以下。

这种方法直接给出了解决方案，只生成一个素数，最多为 1,000,000 的一小部分。

算法在 0.00005 秒（50 微秒）内产生解决方案

代码如下：

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})") 

Answer 2

这里不是 Python 的问题。您的算法效率不高。您需要使用欧拉 phi 函数的属性。其中之一是如果n = p1^a1 * p2^a2 * p3^a3 ... * pk^ak 其中p1, p2, ... ,pk 是素数而a1, a2, ..., ak 是整数，那么：

phi(n) = n * (1 - 1 / p1) * (1 - 1 / p2) * ... * (1 - 1 / pk)

见：https://en.wikipedia.org/wiki/Euler%27s_totient_function#Proof_of_Euler's_product_formula

下面是使用这一事实的更高效的代码：

import math

def Totient(n): # Returns the factors of a positive integer in an array.
    phi = n
    for i in range(2, math.ceil(math.sqrt(n)) + 1):
        if n % i == 0:
            while n % i == 0:
                n /= i
            phi -= phi / i
    if n > 1:
        phi -= phi / n
    return phi

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 maximum value of n/phi(n) is " + str(n_over_phi_MAX) + " at a value of " + str(maxes[-1]))

Answer 3

我测试过，如果我们缓存因子生成结果，整个过程会加快很多。

factors_cache = {}

def FactorsOf(n,include_ends=False):
    key = (n, include_ends)
    if key in factors_cache:
        return factors_cache[key]
    factors = []
    a=0
    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)
    factors_cache[key] = factors
    return factors