We present a \(O\left( N(\log N+\log M) \right)\) solution.

In the original problem, groups must be non-empty and are indistinguishable. Let the answer be \(f(k)\) in such situation.

Consider solving the problem when groups can be empty and are distinguishable. Let the answer be \(g(k)\) in such situation.

\(g(k)\) can be counted by dividing sets of numbers that are the same modulo \(M\) one by one, in the following process.

For each remainder \(r=0,1,\dots,M-1\), let \(S_r=\{i|i\%M=r,1\leq i \leq N\}\). Divide \(S_r\) into \(k\) distinguishable groups, such that each group has at most one number.

Since there are \(M-N\%M\) sets with size \(|S_r|=\lfloor N/M \rfloor\), and \(N\%M\) sets with size \(|S_r|=\lfloor N/M \rfloor +1\), \(g(k)\) can be expressed as the following formula. For each \(k\), if factorials are pre-computed, \(g(k)\) can be evaluated in \(O(\log M)\).

\[ g(k)=\left( {}_kP_{\lfloor N/M \rfloor} \right)^{M-N\%M}\left( {}_kP_{\lfloor N/M \rfloor+1} \right)^{N\%M} \]

Note that \({}_nP_m=\frac{n!}{(n-m)!}\) is the permutation number, the number of ways of dividing \(m\) distinct objects into \(n\) distinct groups, such that each group has at most one object. For \(m>n\) we define \({}_nP_m=0\).

In addition, \(g(k)\) can also be expressed in terms of \(f(k)\). Consider assigning \(k\) distinct IDs to \(k'\leq k\) non-empty, indistinguishable groups, and the remaining \(k-k'\) IDs correspond to empty groups. We can express \(g(k)\) as the following formula.

\[ g(k)=\sum_{k'=0}^{k} {}_{k}P_{k'}f(k') \]

This can be transformed into the form of convolution.

\[ \frac{g(k)}{k!}=\sum_{k'=0}^{k} \frac{f(k')}{(k-k')!} \]

In other words, \([\frac{g(0)}{0!},\frac{g(1)}{1!},\frac{g(2)}{2!},\dots]\) equals to the convolution of \([f(0),f(1),f(2),\dots]\) and \([\frac{1}{0!},\frac{1}{1!},\frac{1}{2!},\dots]\).

By computing all \(g(k)\) in \(O(N\log M)\) and running an inverse convolution in \(O(N\log N)\), we can find all \(f(k)\). Alternatively, let \(c(x)=\frac{1}{0!}+\frac{x}{1!}+\frac{x^2}{2!}+\dots=e^x\). Noting that \(c(x)c(-x)=e^xe^{-x}=1\), we know that the convolutional inverse of \([\frac{1}{0!},\frac{1}{1!},\frac{1}{2!},\dots]\) is \([\frac{1}{0!},-\frac{1}{1!},\frac{1}{2!},\dots]\). Thus \(f(k)\) can be compute from \(g(k)\) by an ordinary convolution as well.

Either way, the overall time complexity is \(O\left( N(\log N+\log M) \right)\).

Sample implementation (C++, 9ms) https://atcoder.jp/contests/abc217/submissions/25619059