Speeding up my solution for Project Euler #650

Question

I've been working on project Euler problem #650, and it generates the correct answer for smaller values, but they are looking for the solution to values up to 20,000. At the speed my code is currently running that will take like a billion years. How can this code be improved and made more efficient/faster?

Here is a link to the problem: https://projecteuler.net/problem=650

And the code:

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

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)

The time complexity is O(n^2/log(n)) and it runs about 20s to get the answer.

Answer 2

me too but you can try this

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