Hints

Hint 1

Recall the formula for finding the number of divisors. If you don't know it, look it up.

Hint 2

The formula for the number of divisors can be expressed as a product. Consider applying a general technique for sums of products.

Hint 3

It would be good to speed up a DP with several states in some way. What kind of states should it have?

Hint 3 Supplement

One possible method is matrix exponentiation. Although there are faster methods, matrix exponentiation with the correct states is sufficient under the constraints of this problem.

Hint 4

The number of prime numbers not greater than $M$ is very small. Yes, really very small.

Hint 4 Supplement

To be more precise, the number of prime numbers not greater than $M$ is so small that $O(8^{\pi(M)}\log N)$ is fast enough, where $\pi(x)$ is the number of primes not greater than $x$.

Hint 5

In other words, using matrix exponentiation, even a DP with $2^{\pi(M)}$ states is feasible.

Solution

As a prerequisite, you need to know how to find the number of divisors of a positive integer \(x\). If \(x\) is factorized as \(p_1^{a_1}p_2^{a_2}\dots p_k^{a_k}\), the number of divisors of \(x\) is given by \((a_1+1)(a_2+1)\dots (a_k+1)\).

By the way, even for the maximum \(M\), there are only six primes not greater than \(M\): \(2, 3, 5, 7, 11, 13\). Therefore, in this problem, the number of divisors of \(x\) can be expressed as follows:

• \((a_1+1)(a_2+1)(a_3+1)(a_4+1)(a_5+1)(a_6+1)\), where \(x\) can be expressed as \(2^{a_1}3^{a_2}5^{a_3}7^{a_4}11^{a_5}13^{a_6}\).

Then, the problem can be rephrased as follows:

There are six towers of blocks. Initially, each tower has exactly one block. Perform the following operation between \(1\) and \(N\) times:

• Choose an integer \(m\) such that \(1 \leq m \leq M\), and stack blocks on the towers as follows:
  • When \(m\) can be expressed as \(2^{a_1}3^{a_2}5^{a_3}7^{a_4}11^{a_5}13^{a_6}\), stack \(a_i\) blocks on the \(i\)-th tower.

Find the sum of the following values over all sequences of operations, modulo \(998244353\):

• \(\prod_{i=1}^{6}b_i\), where \(b_i\) is the number of blocks stacked on the \(i\)-th tower.

It is easy to show that this rephrasing corresponds to the original problem.

By the way, the value to be found for each sequence of operations can be rephrased as “the number of ways to choose exactly one block from each tower.” By performing this block selection during the operations, this problem can be solved using the following bit DP:

• \(DP[i][S] = \) the number of ways such that the set of towers from which blocks have already been chosen matches \(S\) after the \(i\)-th operation.
• \(DP_{ans}[i] = \) the answer after the \(i\)-th operation.

The complexity of this solution is \(O(2^{\pi(M)}N)\), where \(\pi(x)\) is the number of primes not greater than \(x\). Of course, this is far from feasible.

However, the transitions of this DP can easily be converted into matrix exponentiation. Then, the complexity becomes \(O(8^{\pi(M)}\log N)\), which is fast enough.

Judge’s Solution (PyPy3)

from collections import defaultdict

MOD=998244353
MAX_PRIME=6
PRIMES=[2,3,5,7,11,13]

def matrix_pow(a,x):
    mat_size=len(a)
    ret=[[0]*mat_size for _ in range(mat_size)]
    for i in range(mat_size):
        ret[i][i]=1
    while x>0:
        if x&1:
            nret=[[0]*mat_size for _ in range(mat_size)]
            for i in range(mat_size):
                for k in range(mat_size):
                    for j in range(mat_size):
                        nret[i][j]+=ret[i][k]*a[k][j]
                        nret[i][j]%=MOD
            ret=nret
        na=[[0]*mat_size for _ in range(mat_size)]
        for i in range(mat_size):
            for k in range(mat_size):
                for j in range(mat_size):
                    na[i][j]+=a[i][k]*a[k][j]
                    na[i][j]%=MOD
        a=na
        x>>=1
    return ret

N,M=map(int,input().split())
mat=[[0]*((1<<MAX_PRIME)+1) for _ in range((1<<MAX_PRIME)+1)]
nxt=defaultdict(int)
for i in range(1<<MAX_PRIME):
    for j in range(1,M+1):
        add=1
        for k in range(MAX_PRIME):
            if not (i>>k)&1:
                continue
            cnt=0
            while j%PRIMES[k]==0:
                cnt+=1
                j//=PRIMES[k]
            add*=cnt
        nxt[i]+=add
for i in range(1<<MAX_PRIME):
    for j in nxt:
        if (i&j)==0:
            mat[i][i|j]+=nxt[j]
mat[(1<<MAX_PRIME)-1][1<<MAX_PRIME]=1
mat[1<<MAX_PRIME][1<<MAX_PRIME]=1
mat=matrix_pow(mat,N+1)
ans=0
for i in range(1<<MAX_PRIME):
    ans+=mat[i][1<<MAX_PRIME]
print((ans-1)%MOD)

Bonus!: Consider how to solve the problem quickly even when \(M \leq 30\). There seems to be a method using formal power series, but it can also be achieved by replacing matrix exponentiation with a faster method. (Sample Implementation)