Problem: You are given positive integers N and M.  How many sequences a of length N consisting of positive integers satisfy a_1 \times a_2 \times ... \times a_N = M? Find the count modulo 10^9+7.  Here, two sequences a' and a'' are considered different when there exists some i such that a_i' \neq a_i''.

Submission attempts (in order):
from collections import deque

N, M = map(int,input().split())


def make_divisors(n):
    divisors = []
    for i in range(1, int(n**0.5)+1):
        if n % i == 0:
            divisors.append(i)
            if i != n // i:
                divisors.append(n//i)
    return sorted(divisors,reverse=True)


def count(nlist):
    nset = list(set(nlist))
    denom = 1
    for n in nset:
        cnt = nlist.count(n)
        prod = 1
        for i in range(1,cnt+1):
            prod *= i
        denom *= prod
    anom = 1
    for i in range(1,len(nlist)+1):
        anom *= i
    return(anom//denom)




divs = make_divisors(M)
divdict = {str(div):i for i, div in enumerate(divs)} 
divlist = [] 
stack = deque([])
リスト0のところを追加しておこう
for div in divs:
    stack.append([div,M//div])

ans = 0
if N == 2:
    for div in divs:
        if div < M//div:
            ans += 2
        elif div == M//div:
            ans += 1
else:
    
    while len(stack) > 0:
        nlist = stack.pop()
        div = nlist[-2]
        divided = nlist[-1]
        
        ind = divdict[str(div)]
        for div in divs[ind:]:
            if len(nlist) < N and divided%div == 0:
                new = nlist[:-1]+[div, divided//div]
                stack.append(new)
                
                if len(new) == N:
                    prod = nlist[0]
                    for i in range(1,N):
                        prod *= new[i]
                    if new[-1] <= new[-2] and prod==M: 
                        ans += count(new)


print(ans%(10**9+7))


"""
8:30:00
"""

from collections import defaultdict
def getprime(n):
    if not isinstance(n, int):
        raise TypeError("Input int")
    if n < 2:
        raise ValueError("N >= 2")
    prime = []
    
    data = [i+1 for i in range(1,n)]
    while True:
        p = data[0]
        if p >= int(n**0.5):
            return prime+data
        prime.append(p)
        
        data = [d for d in data if d%p != 0]

def factorization(n):
    factors = defaultdict(int)
    primes = getprime(int(n**0.5))
    for prime in primes:
        while n % prime == 0:
            factors[prime] += 1
            n //= prime
            # print(factors,n)
    if n != 1:
        factors[n] += 1
    return factors


def modconb(n,k,mod):
    
    fac = [0]*(n+1)
    finv = [0]*(n+1)
    inv = [0]*(n+1)  

    fac[0] = fac[1] = 1
    finv[0] = finv[1] = 1
    inv[1] = 1

    for i in range(2,n+1):
        fac[i] = fac[i-1]*i % mod
        inv[i] = mod - inv[mod%i] * (mod//i) % mod
        finv[i] = finv[i-1] * inv[i] % mod

    if n<k: return 0
    if n<0 or k<0: return 0
    return(fac[n]*(finv[k]*finv[n-k]%mod)%mod)



mod = 10**9 + 7
n, m = map(int, input().split())
if m == 1:
    print(1)
    exit()
factors = factorization(m)
ans = 1
for k in factors.keys():
    x = factors[k]
    ans *= modconb(x+n-1,x,mod)
    ans %= mod
print(ans)


Verdicts: Wrong Answer

Accepted