別解です．終盤公式解説に合流します．

[1] モーメント母関数とその性質

確率変数 \(X\) に対して，その期待値を \(E[X]\) と書きます．また，次で定まる形式的べき級数 \(M_X\) を \(X\) の モーメント母関数 といいます：

\[M_X(x) = \sum_{i=0}^{\infty} \frac{E[X^i]}{i!}x^i.\]

確率変数の期待値は（収束の問題があるため）すべての確率変数に対して存在するわけではないですが，以下では特に断らずに，\(M_X\) が定義できる確率変数のみを扱うこととします．

\(X, Y\) を独立な確率変数とするとき，モーメント母関数は次の性質を持ちます：

\[M_{X+Y} = M_X\cdot M_Y\]

これは，\(E[(X+Y)^n] = \sum_{i}\binom{n}{i}E[X^iY^{n-i}] = \sum_{i}\binom{n}{i}E[X^i]E[Y^{n-i}]\) から分かります（後者の等号で \(X, Y\) の独立性を用いる）．

[2] 本問題の解法

サイコロを振る回数を \(X\) とすると，本問は \(M_X\) を \(M\) 次まで求める問題であるといえます．

出た目の種類が \(i\) 種類になってから \(i+1\) 種類になるまでにサイコロを振る回数を \(X_i\) とすると，\(X = \sum_{i=0}^{N-1}X_i\) と分解できて，[1] で述べたことより \(M_X = \prod_{i=0}^{N-1} M_{X_i}\) です．

\(X_i\) は，成功確率 \(p_i := \dfrac{N-i}{N}\) の幾何分布で，そのモーメント母関数は \(\dfrac{p_ie^x}{1-(1-p_i)e^x}\) です．よって，本問は次で定まる形式的べき級数を計算する問題に帰着されました：

\[\prod_{i=0}^{N-1}\dfrac{p_ie^x}{1-(1-p_i)e^x}\]

これは \(N\) 次多項式 \(P, Q\) を用いて \(\dfrac{P(e^x)}{Q(e^x)}\) と書けます．\(P, Q\) は \(1\) 次多項式 \(N\) 個の総積なので \(O(N\log^2N)\) 時間で計算できます．なお \(p_i\) が等差数列であることから，第1種スターリング数の計算と同じように \(O(N\log N)\) 時間で計算することもできます．

\(N\) 次多項式 \(P(t), Q(t)\) に対して \(P(e^x), Q(e^x)\) を \(M\) 次まで計算することは公式解説のように \(O(N\log^2N + M\log M)\) 時間でできます（参考）．

以上により本問題は \(O(N\log^2N + M\log M)\) 時間で解けます．