ヒント

ヒント1

約数の個数を求める公式を思い出しましょう。もし知らなければ、調べてみましょう。

ヒント2

約数の公式は何かの積で表せます。積の和典型を適用することを考えましょう。

ヒント3

いくつかの状態をもつDPを何らかの方法で高速化するのがよいでしょう。どのような状態を持てばよいでしょうか。

ヒント3の補足

「何らかの方法」としては、行列累乗などが考えられます。もう少し高速な方法もありますが、本問題の制約下では持つ情報を誤らなければ行列累乗で十分です。

ヒント4

$M$ 以下の素数の個数はとても少ないです。そう、本当にとても少ないです。

ヒント4の補足

より本質的なことを書くと、$M$ 以下の素数の個数は、$\pi(x)$ を $x$ 以下の素数の個数とすると、$O(8^{\pi(M)}\log N)$ が高速なほどに少ないです。

ヒント5

つまり、行列累乗を使えば、$2^{\pi(M)}$ 個の状態を持つDPすら間に合います。

解説

前提として、正整数 \(x\) についての \(x\) の約数の個数の求め方を知っている必要があります。\(x\) の素因数分解された形を \(p_1^{a_1}p_2^{a_2}\dots p_k^{a_k}\) とすると、\(x\) の約数の個数は \((a_1+1)(a_2+1)\dots (a_k+1)\) と表されます。

ところで、\(M\) が最大の場合でも \(M\) 以下の素数は \(2,3,5,7,11,13\) の \(6\) 個のみです。なので、この問題では \(x\) の約数の個数を以下のように表すことができます。

• \(x\) が \(2^{a_1}3^{a_2}5^{a_3}7^{a_4}11^{a_5} 13^{a_{6}}\) と表せるときの、 \((a_1+1)(a_2+1)(a_3+1)(a_4+1)(a_5+1)(a_6+1)\) の値

すると、問題を以下のように言い換えることができます。

積み木のタワーが \(6\) 個ある。はじめ、どのタワーにもちょうど \(1\) 個の積み木がある。このタワーたちに、\(1\) 回以上 \(N\) 回以下、以下に示す操作を行う。

• \(1\leq m\leq M\) なる整数 \(m\) を選び、以下のようにタワーに積み木を積み上げる。
  • \(m\) が \(2^{a_1}3^{a_2}5^{a_3}7^{a_4}11^{a_5}13^{a_6}\) と表せるとき、\(i\) 番目のタワーに \(a_i\) 個の積み木を積む。

すべての操作列に対する以下の値の総和を \(998244353\) で割った余りを求めよ。

• \(i\) 番目のタワーに積まれている積み木の数を \(b_i\) としたときの、\(\prod_{i=1}^{6}b_i\) の値

この言い換えがもとの問題と対応していることは容易に示せます。

ところで、操作列ごとに求める値は、「各タワーからちょうどひとつの積み木を選ぶ通り数」と言い換えることができます。この積み木を選ぶ行為を操作の途中に行うことで、この問題を以下のようなbitDPで解くことができます。

• \(DP[i][S]=i\) 番目の操作まで行ったとき、積み木をすでに選んだタワーの番号の集合が \(S\) と一致するような通り数
• \(DP_{ans}[i]=i\) 番目の操作まで行ったときの答え

この解法の計算量は、\(\pi(x)\) を \(x\) 以下の素数の個数とすると \(O(2^{\pi(M)}N)\) です。もちろん、これは到底間に合いません。

しかし、このDPの遷移は容易に行列累乗にすることができます。すると計算量は \(O(8^{\pi(M)}\log N)\) となり、これは十分高速です。

想定解 (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!: \(M\leq 30\) の場合でも高速に解ける方法を考えてみましょう。形式的冪級数を利用して考察する方法があるようですが、行列累乗をある高速な方法で代替することでも達成できます。(実装例)