公式解説では、（使える数字の種類を \(\sigma\) 種として） \(O(N\sigma)\) 時間のDPで解く方針が解説されていますが、この問題は \(N\) がより大きくても高速に解くことができます。 \(N\) がより大きい場合に用いることができるアルゴリズムを2つ解説します。

この解説では、説明に線形代数や形式的冪級数を用います。

公式解説で述べられている DP は \(N\times\sigma\) の DP テーブルを用いて行われていて、\((\operatorname{dp}[n][1],\ldots,\operatorname{dp}[n][9])\) から \((\operatorname{dp}[n+1][1],\ldots,\operatorname{dp}[n+1][9])\) への遷移の様子は大小比較などを伴わず「定数倍して足す」だけで、\(n\) によらない形になっています。 このような遷移をする DP はある行列を \(1\) つのベクトルに繰り返しかけることに言い換えられます。 実際、

\[ \begin{pmatrix}\operatorname{dp}[n+1][1]\\\operatorname{dp}[n+1][2]\\\operatorname{dp}[n+1][3]\\\operatorname{dp}[n+1][4]\\\operatorname{dp}[n+1][5]\\\operatorname{dp}[n+1][6]\\\operatorname{dp}[n+1][7]\\\operatorname{dp}[n+1][8]\\\operatorname{dp}[n+1][9]\end{pmatrix}=\begin{pmatrix}1&1&0&0&0&0&0&0&0\\1&1&1&0&0&0&0&0&0\\0&1&1&1&0&0&0&0&0\\0&0&1&1&1&0&0&0&0\\0&0&0&1&1&1&0&0&0\\0&0&0&0&1&1&1&0&0\\0&0&0&0&0&1&1&1&0\\0&0&0&0&0&0&1&1&1\\0&0&0&0&0&0&0&1&1\end{pmatrix}\begin{pmatrix}\operatorname{dp}[n][1]\\\operatorname{dp}[n][2]\\\operatorname{dp}[n][3]\\\operatorname{dp}[n][4]\\\operatorname{dp}[n][5]\\\operatorname{dp}[n][6]\\\operatorname{dp}[n][7]\\\operatorname{dp}[n][8]\\\operatorname{dp}[n][9]\end{pmatrix} \]

が成り立ちます。

\(\operatorname{dp}[1][1]=\cdots=\operatorname{dp}[1][9]=1\) なので、

\[ \begin{pmatrix}\operatorname{dp}[n][1]\\\operatorname{dp}[n][2]\\\operatorname{dp}[n][3]\\\operatorname{dp}[n][4]\\\operatorname{dp}[n][5]\\\operatorname{dp}[n][6]\\\operatorname{dp}[n][7]\\\operatorname{dp}[n][8]\\\operatorname{dp}[n][9]\end{pmatrix}=\begin{pmatrix}1&1&0&0&0&0&0&0&0\\1&1&1&0&0&0&0&0&0\\0&1&1&1&0&0&0&0&0\\0&0&1&1&1&0&0&0&0\\0&0&0&1&1&1&0&0&0\\0&0&0&0&1&1&1&0&0\\0&0&0&0&0&1&1&1&0\\0&0&0&0&0&0&1&1&1\\0&0&0&0&0&0&0&1&1\end{pmatrix}^{n-1}\begin{pmatrix}1\\1\\1\\1\\1\\1\\1\\1\\1\end{pmatrix} \]

が成り立ちます（帰納法で示せます）。

サイズ \(\sigma\) の正方行列の積は \(O(\sigma^3)\) 回の要素の積で計算できるので、繰り返し二乗法によってこの問題を \(O(\sigma^3\log N)\) 時間で解くことができます。

さらに高速に解くこともできます。

\(a_n\) を \(\operatorname{dp}[n][1]+\cdots+\operatorname{dp}[n][9]\) とします。つまり、

\[ a_n=\begin{pmatrix}1&1&1&1&1&1&1&1&1\end{pmatrix}\begin{pmatrix}1&1&0&0&0&0&0&0&0\\1&1&1&0&0&0&0&0&0\\0&1&1&1&0&0&0&0&0\\0&0&1&1&1&0&0&0&0\\0&0&0&1&1&1&0&0&0\\0&0&0&0&1&1&1&0&0\\0&0&0&0&0&1&1&1&0\\0&0&0&0&0&0&1&1&1\\0&0&0&0&0&0&0&1&1\end{pmatrix}^{n-1}\begin{pmatrix}1\\1\\1\\1\\1\\1\\1\\1\\1\end{pmatrix} \]

です。 簡単のため、\(a_n=\bm{e}^\top A^{n-1}\bm{e}\) と書くことにします。

形式的冪級数 \(F(x)\coloneqq a_0+a_1x+a_2x^2\cdots\) とします。

いま、ある \(c_0,\ldots,c_d\) に対して \(c_dA^d+\cdots+c_0A^0=O\) が成り立っているとし、\(P(x)\coloneqq c_d+\cdots+c_0x^d\) とします。

ここで、

\[ \begin{aligned} F(x)P(x)=\;&a_0c_d+(a_1c_d+a_0c_{d-1})x+\cdots+(a_{d-1}c_d+\cdots+a_0c_1)x^{d-1}+\\ &(a_dc_d+\cdots+a_0c_0)x^d+\cdots+(a_{d+k}c_d+\cdots+a_kc_0)x^{d+k}+\cdots \end{aligned} \]

となります。\(d+k\) 次の項を考えると（\(k\geq0\)）

\[ \begin{aligned} a_{d+k}c_d+\cdots+a_kc_0&=\bm{e}^\top c_dA^{d+k-1}\bm{e}+\cdots+\bm{e}^\top c_0A^{k-1}\bm{e}\\ &=\bm{e}^\top (c_dA^d+\cdots+c_0A_0)A^{k-1}\bm{e}\\ &=\bm{e}^\top OA^{k-1}\bm{e}=0 \end{aligned} \]

となるため、\(F(x)P(x)\) はたかだか \(d-1\) 次の式になることがわかります。

Cayley-Hamilton の定理より、\(\det(A-xI)\) は \(P(x)\) の条件を満たします。

事前に \(P(x)\) および \(F(x)P(x)\) を求めておくことで Bostan-Mori 法などで \(O(\sigma\log\sigma\log N)\) 時間で答えを求めることができます。

また、\(a_n=\bm{x}^\top Y^n\bm{z}\) の形で書けるとわかった時点で（愚直などで）小さいところの値を計算し、 Berlekamp–Massey algorithm を用いることでも \(P(x)\) および \(F(x)P(x)\) を求めることができます。