需要帮助加快我的 Project Euler #650 解决方案

Question

所以我一直在研究欧拉问题#650 项目，它为较小的值生成正确的答案，但他们正在寻找高达 20,000 的值的解决方案。 按照我的代码当前运行的速度，大约需要十亿年。 任何改进此代码并使其更高效/更快的提示将不胜感激。 这是问题的链接：https://projecteuler.net/problem=650

和代码

from scipy.special import comb
import math
import time

t0 = time.time()

def B(x):
    #takes product of binomial coefficients
    product = 1
    for i in range(x + 1):
        product *= comb(x, i, exact=True)
    return product



def D(y):
    #Sums all divisors of the product of binom. coefficients
    x = B(y)
    summation = []
    for i in range(1, int(math.sqrt(x)) + 1):
        if x % i == 0:
            summation.append(i)
            summation.append(x // i)
    summation = list(dict.fromkeys(summation))
    return sum(summation)


def S(z):
    """sums up all the sums of divisors of the product of the binomial 
    coefficients"""
    sigma = []
    for i in range(1, z + 1):
        sigma.append(D(i))
    return sum(sigma)

print(S(20000))

t1 = time.time()

total = t1 - t0


print(total)

Answer 1

我也是 但你可以试试这个

from math import factorial

def n_choose_k(n, k):
    if k == 0:
        return 1
    else:
        return factorial(n)//((factorial(k))*factorial(n-k))

def B(n):
    prods = []
    for k in range(0, n):
        prods.append(n_choose_k(n, k))

Answer 2

import copy
import primefac
import collections
import tqdm

primes = list(primefac.primegen(20000))
prime_factors = [{}]
prime_factors_factorials = [collections.defaultdict(int)] # prime factor of n!
for n in tqdm.tqdm(range(2,20001)):
    cur = copy.copy(prime_factors_factorials[-1])
    prime_factor = collections.Counter(primefac.primefac(n))
    for p,v in prime_factor.items():
        cur[p] += v
    prime_factors.append(prime_factor)
    prime_factors_factorials.append(cur)

def sum_factors(prime_fac_dict,mod):
    ans = 1
    for p,v in prime_fac_dict.items():
        ans *= (pow(p, v+1, mod) - 1)*pow(p-1, -1, mod)
        ans %= mod
    return ans

mod = 10**9 + 7
n = 20000
cur = collections.defaultdict(int)
ans = 1
for i in tqdm.tqdm(range(2,n+1)):
    add_factors = prime_factors[i - 1]
    for p,v in add_factors.items():
        cur[p] += (i - 1)*v
    minus_factors = prime_factors_factorials[i - 2]
    for p,v in minus_factors.items():
        cur[p] -= v
    ans += sum_factors(cur, mod)
    ans %= mod

print(ans)

时间复杂度为O(n^2/log(n))，运行大约20s即可得到答案。