一眼非常典的题，但是冲了数个小时问题在于我的状态设计。感觉还是得加练一些思维，但是时间不多了，中考如果不爆炸的话高一加油吧。感觉讲题的 AC 自动机刻画状态的方法比较牛，但是我没听懂那个减状态的优化 /yun。

这篇题解前半部分纯改模拟赛题的，等价于单测的 \(k\le 200\) 做法，如果需要看生成函数部分可以直接跳过第一部分。

直接矩阵快速幂的做法

首先，注意到颜色数和要求连续段数都是 \(m\)，我们考虑怎么样才能运用这个性质。注意到这个 \(n\le 10^9\)，则这提示我们需要设计出一个简易的 DP 式子，然后使用矩阵快速幂优化 DP 的方式才能通过此题（？为什么 \(m\le 200\)，这很眉目你知道吗）。

我们考虑求出长度为 \(n\) 的不符合填充数量（也就是有至少一个子序列是一个 \(m\) 排列）。这个东西的状态设计比较难想（对我这种区来说，我想到的做法和 AC_lover 是一样的，但是没调出来），我们设 \(f_{i,j}\) 表示前 \(i\) 个元素，且后 \(j\) 个元素两两不同的方案数。

  • 首先考虑接一个出现在这里的元素。对于任意出现的元素，我们都可以找到前 \(j\) 个当中的一个位置，满足这个位置往后开始是符合条件的，也就是说 \(\forall 1\le k\le j,f_{i+1,k}\leftarrow f_{i+1,k}+f_{i,j}\)。
  • 你也可以选择加入一个新的元素，有 \(m-j\) 种选择，\(f_{i+1,j+1}\leftarrow f_{i+1,j+1}+(m-j)f_{i,j}\)。

由于这个东西的本质是每一位都计数全部 \(m\) 种情况，所以其是不重不漏的。那么我们就考虑直接对这个东西跑矩阵快速幂优化 DP。刷表法好想但是难转化，改一下：

\[f_{i,j}=(m-(j-1))f_{i-1,j-1}+\sum_{k=j}^{m-1} f_{i-1,k} \]

常规套路，使用矩阵快速幂的时候不考虑 \(i\) 维。然后我们看怎么填数：

\[f_j\leftarrow (m-(j-1)) pre_{j-1}+\sum_{k=j}^{m-1} pre_k \]

则考虑转移矩阵的填写方式。\(j-1\to j\) 的点显然都是满足 \((j-1,j)\) 的，此时这个点的贡献直接写就是 \((m-(j-1))\)。这里的系数是 \((m-j+1)\)。然后 \((k,j)\) 的点全部满足 \(j\le k\)，也就是在矩阵中满足 \(i\ge j\) 的点。这些点的贡献都是 \(1\)（一个经典套路是 \((i,j)\) 处的贡献对应 \(i\to j\) 的转移）。这样的矩阵量级就刚刚好能跑。

注意到我们并没有正难则反，所以最终跑容斥的时候要刻意避开 \(f_m\) 的转移去写，这样就能保证转移的正确性。最终时间复杂度 \(\mathcal O(Tm^3\log m)\)，T 了的话利用表格的 \(0\) 很多的性质加个常数优化就能过。

::::success[模拟赛能过的代码（很慢，1000ms）]

#define FileIO
#include
using namespace std;
#define ll long long
#define fi first
#define se second
#define il inline
bool stMem; double stT=clock();
const int N=1e5+10;
const ll P=998244353;
il ll qp(ll a,ll b=P-2){ll res=1;for(;b;b>>=1,a=(a*a)%P) if(b&1) res=(res*a)%P;return res;
}
// ---------------------------------------------------- //
const int M=205;
ll n,m;
struct matrix{ll a[M][M];il void init_empty(){for(int i=0;i=j) a[i][j]=1;else a[i][j]=0;}return;}
} trans, dp;
il matrix operator*(matrix a,matrix b){matrix c; c.init_empty();for(int k=1;k>=1,a=a*a) if(b&1) c=c*a;return c;
}
// ---------------------------------------------------- //
il void solve(){cin>>n>>m;if(m==1){ cout<<0<<"\n"; return; }dp.init_empty(); dp.a[1][1]=m; // 这里是 f[1]，矩阵快速幂只需要 n-1 层trans.init_trans();dp=dp*(trans^(n-1));ll ans=0;for(int j=1;j>T;while(T--) solve();cerr<<"\n";cerr<<"[Info] End Of Running. \n";cerr<<"       Time Used   : "<<(clock()-stT)/CLOCKS_PER_SEC*1000.00<<" ms\n";cerr<<"       Memory Used : "<<(&stMem-&edMem)/1048576.00<<" Mib\n";return 0;
}

::::

生成函数优化

你注意到当 \(m\)（对应原题的 \(k\)）开到 \(10^4\) 之后这样写就不可能能过了。且我们目前的状态设计的复杂度是 \(\mathcal O(nm)\) 的。我们考虑将一位的全部合法状态集合起来为 \(f_i\)，那么有如下生成函数设计：

\[F(x)=\sum_{i=0}^\infty f_ix^i \]

我们使用这个生成函数刻画该 DP 的原理在于对最终函数进行降幂处理，以使得其能刻画最终答案。

\[F(x)=\frac{P(x)}{Q(x)} \]

::::info[题解中的 \(P(x)\) 与 \(Q(x)\) 设计]{open}

  • \(P(x)\) 的设计：

    \[P(x)=\sum_{i=0}^{m-1}i!x^i \]

  • \(Q(x)\) 的设计：

    \[Q(x)=1-\sum_{i=1}^{m-1}(m-i)(i-1)!x^i \]

::::

这非常难懂。我表述一下这个东西是怎么来的。首先，我们考虑刻画非法串如何保证不重不漏。首先，这个东西是方便抵消掉 \(x^i\)，最简单的式子貌似并不是这个，而是另一条比较简易的式子，然后给分子配上 \(P(x)\) 就能自然得出 \(Q(x)\)。

首先我们延续刚刚矩阵快速幂的思路来设计 \(f_j\)。根据我们上面的转移式子，新接的字符可以让 \(\forall 1\le k\le j,f_j\to f_k\)，以及 \(f_j\to f_{j+1}\)。当然这些状态都是基于 \(j\neq m-1\) 的，否则绝对非法。

那么我们就考虑 \(F_i(x)\) 为最终停在 \(j=i\) 的串的生成函数，则我们要求的 \(F(x)=\sum_{i=0}^{m-1} F_i(x)\)，且有 \(F(0)=1\)（矩阵快速幂初值）。此时根据转移矩阵中 \(i\ge j\) 的点，我们显然需要求一个后缀和状物的东西 \(T_i(x)=\sum_{j\ge i} F_j(x)\)，以处理接重复后缀的情况（这里是递推式的，和刷表式相反），则直接根据转移所获得的状态递推式就是：

\[F_i(x)=xT_i(x)+x(m-i+1)F_{i-1}(x),1\le i\le m-1 \]

那么根据这个 \(T_i(x)\) 的定义有 \(F_i(x)=T_{i}(x)-T_{i+1}(x)\)，直接暴力带入以下三式：\(F_i(x)=T_{i}-T_{i+1},T_0(x)=F(x),T_1(x)=F(x)-1\)（注意这里后两式带的是 \(F(x)\) 而不是 \(F_i(x)\)，两者有本质区别），得：

\[\begin{aligned} T_{i+1}(x)&=(1+(m-i)x)T_{i}(x)-x(m-i+1)T_{i-1}(x)\\ T_0(x)&=F(x)\\ T_{1}(x)&=F(x)-1\\ T_m(x)&=0 \end{aligned} \]

然后你就发现各 \(T_i(x)\) 都是 \(F(x)\)（注意不是 \(F_i(x)\)） 的一次函数，且 \(T_i(x)\) 构成与自己形成递推关系，也就是二阶线性递推。那么我们直接钦定（本质是暴力展开 /汗颜）：

\[T_i(x)=K_i(x)F(x)+B_i(x) \]

其中我们的初值就是 \(K_0(x)=1,B_0(x)=0,K_1(x)=1,B_1(x)=-1\)（也就是 \(T_0(x),T_1(x)\) 的函数关系）。那么我们直接根据 \(T_m(x)=0\) 就可以解出 \(F=\frac{-B_m(x)}{K_m(x)}\)，于是就有 \(P(x)=-B_m(x),Q(x)=A_m(x)\)。然后就有小朋友要问了那怎么求 \(A_i(x),B_i(x)\) 呢？考虑到我们的 Token 非常的多，暴力带入计算即可。

当然这里直接暴力去算是非常对的，但是太费草稿纸了。带出来就是上面提到过的 \(P(x),Q(x)\)。下面会介绍一个比较好懂的计算方法。

::::success[另一种计算方法]{open}

有点难懂，参考了 AI，但是内容都是我自己总结的。

我们注意到组合意义天地灭，但是代数推导不够纸。首先我们考虑写容斥的思路。钦定 \(S(w)\) 为所有非法区间左端点的点集。那么根据我们非常好用的二项式定理可以写出容斥：

\[[|S|=0]=\sum_{s\subseteq S(w)}(-1)^{|s|} \]

（注意区分这里大小写的区别）然后我们直接根据这个写出 \(F(x)\) 的式子：

\[F(x)=\sum_{w}x^{|w|}\sum_{s\subseteq S(w)}(-1)^{|s|}=\sum_{(w,S)}(-1)^{|S|}x^{|w|} \]

也就是全部都要让 \(|S|=0\)，保证不存在排列。那么我们转化这个 \(S(w)\)，把不同的 \([w,w+m-1]\) 拆成若干个连通块，那么我们使用这个块序列的思想来刻画 \(S\)：

  • 对于没标记的数，也就是合法拼接，有 \(mx\) 作为系数，代表其可以作为任意数字。
  • 否则构成一个连通块，其贡献到其所属的连通块内，设 \(C(x)=\sum_{s}(-1)^{|S|}x^{\text{len}}\)。

那么我们直接对这个东西跑一下生成函数，由任意标记情况都能刻画答案，有：

\[F(x)=\frac{1}{1-mx-C(x)} \]

然后我们拆这个 \(C(x)\)。目前来看这还是一个感性理，设 \(s\) 所代表的非法段集合内左端点为 \(s_1，两段为一个连续段当且仅当 \(1\le (d_t=s_{t+1}-s_t) 。刻画这个连续段情况：

  • 第一个非法段是排列，则 \(k!\) 种。
  • 两个非法段可以看做是一段把另一段吃进去了，中间还缺一些值没填，本质是一个排列，则 \(d_t!\) 种，也就是：

然后我们就能写出大份式子了：

\[C(x)=\sum_{i\ge 1}m!x^m\prod_{t=1}^{i-1}\left(\sum_{d=1}^{m-1} m!x^m\right) \]

这一大坨是什么意思呢？若干连续段所代表的生成函数直接卷积在一起。非常诡异，但是就是符合上述刻画的。然后开始大份推式子（下文钦定 \(G(x)=\sum_{i=1}^{m-1} i!x^i\)）：

\[\begin{aligned} C(x)&=\sum_{i\ge 1}(-1)^im!x^m\prod_{t=1}^{i-1}\left(\sum_{d=1}^{m-1} d!x^d\right)\\ &=m!x^m\sum_{i\ge 1}(-1)^i\prod_{t=1}^{i-1}\left(\sum_{d=1}^{m-1}d!x^d\right)\\ &=m!x^m\sum_{i\ge 1}(-1)^i\prod_{t=1}^{i-1} G(x)\\ &=m!x^m\sum_{i\ge 1}(-1)^i G(x)^{i-1}\\ &=m!x^m\left(-1+G(x)-G(x)^2+\cdots\right)\\ &=m!x^m\cdot-\left((-G(x))^0+(-G(x))^1+(-G(x))^2+\cdots\right)\\ &=m!x^m\left(-\frac{1}{1-(-G(x))}\right)\\ &=-\frac{m!x^m}{1+G(x)} \end{aligned} \]

这些转化都比较显然，这里不再赘述。然后带入最开头的大份式子，可得：

\[\begin{aligned} F(x)&=\frac{1}{1-mx-C(x)}\\ &=\frac{1}{1-mx+\frac{m!x^m}{1+G(x)}}\\ &=\frac{1+G(x)}{(1-mx)(1+G(x))+m!x^m}\\ \end{aligned} \]

即：

\[P(x)=1+G(x)=1+\sum_{i=1}^{m-1} i!x^i=\sum_{i=0}^{m-1}i!x^i \]

以及：

\[Q(x)=1-\sum_{i=1}^{m-1}(m-i)(i-1)!x^i \]

底下这个 \(Q(x)\) 是等价的，我有点懒了，不想打 LaTax，直接复制了 qwq。

::::

然后我们现在的目标就是求 \([x^n]F(x)=[x^n]\frac{P(x)}{Q(x)}\)。然后我们使用不知道是什么的 Bostan-Mori 算法。注意到其就是专门求解 \([x^n]\frac{P(x)}{Q(x)}\) 的算法（且可以满足 \(n\) 很大的约束），那么我们直接套用。

  • 首先分子分母同时乘 \(Q(-x)\)，得到 \(\frac{P(x)Q(-x)}{Q(x)Q(-x)}\)。
  • 对这个东西整理，得 \(\frac{U_0(x^2)+xU_1(x^2)}{V(x^2)}\)。
  • 然后注意到只与 \(x^2\) 相关，那么递归下去跑卷积即可。

复杂度 \(\mathcal O(n\log n\log m)\)，递归一只 \(\log\) 卷积一只 \(\log\)。

感觉这题神了吧。

::::success[参考代码（含 poly，极度卡顿警告）]

//#define FileIO
#include
using namespace std;
#define ll long long
#define fi first
#define se second
#define il inline
bool stMem; double stT=clock();
// ---------------------------------------------------- //
namespace exfish_poly{#define ll long long#define ld long double#define il inlineconst ll P=998244353, G=3, iG=332748118;const ll N=(1<<22)+5;const ld pi=acos(-1);// model: mod int (998244353)ll qp(ll a,ll b=P-2){ll res=1;for(;b;b>>=1,a=(a*a)%P) if(b&1) res=(res*a)%P;return res;}il ll add(ll a,ll b){ return a+b>=P?a+b-P:a+b; }il ll sub(ll a,ll b){ return a>1]>>1|((i&1)?fn>>1:0);}il void fft(Vector *f, int op){ // +fft -ifftfor(int i=0;i>1);for(int s=0;s>1]>>1|((i&1)?fn>>1:0);}il void ntt(ll *f, int op){ // +fntt -ifnttfor(int i=0;i>1);for(int s=0;s>1);ntt::ntt_init(n<<1);copy_to(f,f+n+1,a);clear(a+n+1,a+ntt::fn); clear(g+n+1,g+ntt::fn);ntt::ntt(a,1); ntt::ntt(g,1);for(int i=0;i>1;cdq(f,g,l,mid);multiple(f,r-l,g+l,mid-l,a); // [l,mid] -> [mid+1,r]for(int i=mid+1;i<=r;++i) g[i]=add(g[i],a[i-l]);cdq(f,g,mid+1,r);}// model: divisive (remainder)ll irg[N];il void reverse(ll *f,ll *rf,int n){ for(int i=0;i<=n/2;++i){ ll t=f[i]; rf[i]=f[n-i]; rf[n-i]=t; } }il void divisive(ll *f,int n,ll *g,int m,ll *q,ll *r){reverse(f,q,n); reverse(g,g,m);inverse(g,irg,n-m);multiple(q,n,irg,n-m,q); // irg(x)·rf(x) clear(q+n-m+1,q+ntt::fn);reverse(q,q,n-m); reverse(g,g,m);multiple(g,n,q,n-m,a); // g(x)q(x)for(int i=0;i>1);inverse(g,ig,n);multiple(f,n,ig,n,a);for(int i=0;i<=n;++i) g[i]=add(g[i],a[i])*qp(2)%P;clear(g+n+1,g+ntt::fn);}// model: ex-sqrt (bsgs)il ll bsgs(ll x,ll y){ll lmt=sqrtl(P)+1, mul=1;unordered_mapmp;for(int i=0;i>1);inverse(g,ig,n);multiple(f,n,ig,n,a);for(int i=0;i<=n;++i) g[i]=add(g[i],a[i])*qp(2)%P;clear(g+n+1,g+ntt::fn);}	// model: calculusil void derivative(ll *f,ll *g,int n){ for(int i=1;i<=n;++i) g[i-1]=f[i]*i%P; g[n]=0; }il void integral(ll *f,ll *g,int n){ for(int i=n;~i;--i) g[i+1]=f[i]*qp(i+1)%P; g[0]=0; }// model: lnll fi[N];il void ln(ll *f,ll *g,int n){inverse(f,fi,n);derivative(f,g,n);multiple(g,n,fi,n,a);clear(g+n+1,g+ntt::fn);integral(a,g,n-1);}// model: expll lg[N];void exp(ll *f,ll *g,int n){if(!n){ g[0]=1; return; }exp(f,g,n>>1);ln(g,lg,n);for(int i=0;i<=n;++i) lg[i]=sub((i==0)+f[i],lg[i]); // notice: only x^0multiple(g,n,lg,n,g);clear(g+n+1,g+ntt::fn);}// model: quick powerll lf[N];void qpow(ll *f,ll *g,int n,ll k){ln(f,lf,n);for(int i=0;i<=n;++i) lf[i]=lf[i]*k%P;exp(lf,g,n);}// model: sin / cosll eg[N],ieg[N],_i,r_i;il void sin(ll *f,ll *g,int n){_i=residue(P-1), r_i=P-_i;for(int i=0;i<=n;++i) g[i]=f[i]*_i%P;exp(g,eg,n);inverse(eg,ieg,n);for(int i=0;i<=n;++i) g[i]=(eg[i]-ieg[i]+P)%P*qp(2)%P*r_i%P;}il void cos(ll *f,ll *g,int n){_i=residue(P-1);for(int i=0;i<=n;++i) g[i]=f[i]*_i%P;exp(g,eg,n);inverse(eg,ieg,n);for(int i=0;i<=n;++i) g[i]=(eg[i]+ieg[i])%P*qp(2)%P;}// model: arcsin / arctanll sqg[N],pln[N],nln[N];il void arcsin(ll *f,ll *g,int n){_i=residue(P-1), r_i=P-_i;multiple(f,n,f,n,g); // x^2for(int i=0;i<=n;++i) g[i]=sub((i==0),g[i]); // notice: only x^0square_root(g,sqg,n);for(int i=0;i<=n;++i) sqg[i]=(sqg[i]+_i*f[i]%P)%P;ln(sqg,g,n);for(int i=0;i<=n;++i) g[i]=g[i]*r_i%P;}il void arctan(ll *f,ll *g,int n){_i=residue(P-1), r_i=P-_i;for(int i=0;i<=n;++i) g[i]=add((i==0),_i*f[i]%P); // notice: only x^0ln(g,pln,n); for(int i=0;i<=n;++i) g[i]=sub((i==0),_i*f[i]%P); // notice: only x^0ln(g,nln,n); for(int i=0;i<=n;++i) g[i]=sub(pln[i],nln[i])*qp(2)%P*r_i%P;}// model: linear recurrence (fiduccia)ll tq[N],tr[N],tm[N],md[N],pa[N];il ll recurrence_fiduccia(ll *f,ll *p,int k,ll n){for(int i=0;i<=k;++i) md[i]=(P-(f[k-i-1]%P+P)%P)%P, pa[i]=p[i]; md[k]=1;for(int i=0;i<=k;++i) f[i]=p[i]=0; f[1]=1; p[0]=1;for(;n;n>>=1){if(n&1){multiple(p,k,f,k,tm);divisive(tm,2*k,md,k,tq,tr);copy_to(tr,tr+k+1,p);}multiple(f,k,f,k,tm); 	divisive(tm,2*k,md,k,tq,tr);copy_to(tr,tr+k+1,f);}ll res=0;for(int i=0;i>1]=rp[i]; psz=(i>>1)-1;for(i=0  ;i<=qsz+qsz;i+=2) q[i>>1]=rq[i]; qsz=(i>>1)-1;n>>=1;}return p[0]*qp(q[0])%P;} il ll recurrence_bostan_mori(ll *f,ll *p,int k,ll n){for(int i=0;i<=k;++i) md[i]=(P-(f[i-1]%P+P)%P)%P, pa[i]=(p[i]%P+P)%P; md[0]=1;multiple(md,k,pa,k-1,tp);return bostan_mori(tp,k-1,md,k,n);}
}
using namespace exfish_poly;
ll n,m,p[N],q[N],fac[N];
// ---------------------------------------------------- //
il void solve(){cin>>n>>m;for(int i=0;i>T;while(T--) solve();cerr<<"\n";cerr<<"[Info] End Of Running. \n";cerr<<"       Time Used   : "<<(clock()-stT)/CLOCKS_PER_SEC*1000.00<<" ms\n";cerr<<"       Memory Used : "<<(&stMem-&edMem)/1048576.00<<" Mib\n";return 0;
}

::::

::::warning[AIGC]{open}

蓝色大肥鱼一点都不好用，骗我推的半天的式子跟我说她错完了。还是 Arena 好用 /hec /hec。

本文完全由人类生成，辅助理解的部分用了 AI，因为我实在是无法看懂你们的题解 /xk。 ::::

https://www.luogu.com.cn/article/uid99ij1