数据范围告诉我们要写两档的分
第一档:$M\leq200,N\leq10^9$,可以枚举$m$计算答案
直接矩阵快速幂:$O\left(M^4\log_2N\right)$,会超时,所以我们需要某些“技巧”来加速这个过程:矩阵特征多项式
矩阵$A$的特征多项式为$f(\lambda)=\left|\lambda I-A\right|$
Cayley–Hamilton定理指出,如果将$A$作为自变量代入特征多项式,那么$f(A)=0$
证明过于复杂,请自行wikipedia或找书(估计这个坑不会填)
2018.4.6 来填坑
果然我一点都不懂线性代数,还是滚回去重新学吧==
1.把$n$阶行列式$a$划去第$i$行第$j$列,形成$n-1$阶行列式,称其为$a_{ij}$的余子式$A_{ij}$
2.如果行列式有两行相等,那么它的行列式为$0$
用数学归纳法,结论对二阶行列式显然成立,假设结论对$n-1$阶行列式成立,且$n$阶行列式$a$的第$i$行和第$j$行相等,选取一个$k$使得$k\ne i$且$k\ne j$,把$n$阶行列式按第$k$行展开,结果为$\begin{align*}\sum\limits_{l=1}^na_{kl}A_{kl}\end{align*}$,其中每个$A_{kl}$都是$n-1$阶行列式,且有两行相等,所以$A_{kl}=0$,即$a$的行列式为$0$
对列的证明类似
3.行列式某一行的元素乘另一行对应元素的余子式之和为$0$
对于原行列式$a=\left|\begin{matrix}a_{11}&\cdots&a_{1n}\\\vdots&&\vdots\\a_{i1}&\cdots&a_{in}\\\vdots&&\vdots\\a_{j1}&\cdots&a_{jn}\\\vdots&&\vdots\\a_{n1}&\cdots&a_{nn}\end{matrix}\right|$构造行列式$b=\left|\begin{matrix}a_{11}&\cdots&a_{1n}\\\vdots&&\vdots\\a_{i1}&\cdots&a_{in}\\\vdots&&\vdots\\a_{i1}&\cdots&a_{in}\\\vdots&&\vdots\\a_{n1}&\cdots&a_{nn}\end{matrix}\right|$,显然这个行列式值为$0$
因为两个行列式只有第$j$行不同,所以$A_{j1}\cdots A_{jn}$和$B_{j1}\cdots B_{jn}$全相同,对$b$按第$j$行展开,有$\begin{align*}\sum\limits_{k=1}^na_{ik}B_{jk}=0\end{align*}$,即$\begin{align*}\sum\limits_{k=1}^na_{ik}A_{jk}=0\end{align*}$
对列的证明类似
4.对于矩阵$A$,设其行列式$|A|$中$a_{ij}$的余子式为$A_{ij}$,定义$A$的伴随矩阵$A^*=\left[\begin{matrix}A_{11}&A_{21}&\cdots&A_{n1}\\A_{12}&A_{22}&\cdots&A_{n2}\\\vdots&\vdots&&\vdots\\A_{1n}&A_{2n}&\cdots&A_{nn}\end{matrix}\right]$,我们可以得到$A^*A=|A|I$(直接乘起来即可,主对角线的每一个位置是$|A|$的对应行元素乘对应行余子式,值为$|A|$,非对角线上的值由3可得是$0$)
5.于是我们可以证Cayley-Hamilton定理了
对$n$阶矩阵$A$,设$B(\lambda)$为$(\lambda I-A)$的伴随矩阵,$f(\lambda)=|\lambda I-A|$为$A$的特征多项式,则$B(\lambda)(\lambda I-A)=|\lambda I-A|I=f(\lambda)I$
先看左边,注意到$B$的每一个元素都是$|A|$的一个余子式,也就是$n-1$阶行列式,所以$B$的每一个元素都是关于$\lambda$的次数不超过$n-1$的一元多项式,我们把它写成$\begin{align*}B=\sum\limits_{i=0}^{n-1}\lambda^iB_i\end{align*}$,其中$B_i$是$n$阶矩阵,那么$\begin{align*}B(\lambda)(\lambda I-A)=\sum\limits_{i=0}^{n-1}\lambda^{i+1}B_i-\lambda^iB_iA\end{align*}$
再看右边,因为$\begin{align*}f(\lambda)=\sum\limits_{i=0}^na_i\lambda^i\end{align*}$,所以$\begin{align*}f(\lambda)I=\sum\limits_{i=0}^na_i\lambda^iI\end{align*}$
比较两个式子,把$\lambda^i$约掉,我们得到$\begin{cases}B_{n-1}=a_nI\\B_{n-2}-B_{n-1}A=a_{n-1}I\\\cdots\\B_0-B_1A=a_1I\\-B_0A=a_0I\end{cases}$
用$A^n\cdots I$依次乘上面式子的第$1\cdots n+1$项,我们得到$\begin{cases}B_{n-1}A^n=a_nA^n\\B_{n-2}A^{n-1}-B_{n-1}A^n=a_{n-1}A^{n-1}\\\cdots\\B_0A-B_1A^2=a_1A\\-B_0A=a_0I\end{cases}$
全部加起来,我们得到$f(A)=0$,就证完了
容易发现如果$A$是$m$阶矩阵,那么$f(\lambda)$是一个$m$次多项式
假如$g(x)=q(x)p(x)+r(x)$且$p(x)$是$m$阶矩阵$A$的特征多项式,并且$r(x)$的次数小于$p(x)$的次数,那么将$A$代入上式,得到$g(A)=r(A)$
也就是说,如果要计算$g(A)$,把$g(x)$对$A$的特征多项式取模得到的$r(x)$满足$r(A)=g(A)$
有什么用?假如$g(x)$的次数很大,使得我们难以计算$g(A)$,而$m$比较小,那么我们可以算出次数更小的$r(x)$,把问题转化为计算$r(A)$
对于常系数线性递推,我们要求数列$f$的第$n$项,如果转移矩阵是$A$,那么相当于求$A^{n-1}$,直接做是$O\left(m^3\log_2n\right)$的,如果我们先预处理出$f_{1\cdots m}$,用快速幂计算$r(\lambda)=\lambda^{n}\%f(\lambda)$再根据它的系数计算答案(不用把$A$代进去是因为$r(\lambda)$的次数$\lt m$,只需使用$f_{1\cdots m}$即可算出答案)
对于这题,$f_n=(a-1)\sum\limits_{k=1}^mf_{n-k}$,转移矩阵如下
$$A=\left[\begin{matrix}0&0&0&\cdots&0&a-1\\1&0&0&\cdots&0&a-1\\0&1&0&\cdots&0&a-1\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&0&a-1\\0&0&0&\cdots&1&a-1\end{matrix}\right]$$
$$\begin{align*}f(\lambda)=\left|\lambda I-A\right|&=\left|\begin{matrix}\lambda&0&0&\cdots&0&1-a\\-1&\lambda&0&\cdots&0&1-a\\0&-1&\lambda&\cdots&0&1-a\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&\lambda&1-a\\0&0&0&\cdots&-1&\lambda-a+1\end{matrix}\right|\\&=\left|\begin{matrix}\lambda&0&0&\cdots&0&1-a\\0&\lambda&0&\cdots&0&\frac{\lambda+1}\lambda(1-a)\\0&0&\lambda&\cdots&0&\frac{\lambda^2+\lambda+1}{\lambda^2}(1-a)\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&\lambda&\frac{\lambda^{m-2}+\cdots+1}{\lambda^{m-2}}(1-a)\\0&0&0&\cdots&0&\lambda-a+1+\frac{\lambda^{m-2}+\cdots+1}{\lambda^{m-1}}(1-a)\end{matrix}\right|\\&=\lambda^m+(\lambda^{m-1}+\cdots+1)(1-a)\end{align*}$$
有些题目手算是算不出来的,需要插值
还有一些题目更加丧病,多项式取模和多项式乘法要求用NTT做,这个就看着办吧2333
枚举$m$,对于每个$m\lt M$如此计算,对于$m\geq M$,答案都是一样的,对$19260817$等比数列求和即可,我们在$O\left(M^3\log_2N\right)$的时间内解决此题
第二档:$M\leq10^9,N\leq3\times10^6$
这个就直接多了,直接推生成函数算出通项公式就好了
首先令$n\leq0$的$f_n=0$,因为$f_n-f_{n-1}=(a-1)(f_{n-1}-f_{n-m-1})$,所以$f_n=af_{n-1}-(a-1)f_{n-m-1}$
记$F_m(z)$为这个$m$下的$f$的生成函数(注意:以上变形的递推公式在$n=1$和$n=m+1$时不适用,应该补齐)
$$\begin{align*}F_m(z)&=\sum\limits_{i=0}f_iz^i\\F_m(z)&=az-az^{m+1}+\sum\limits_{i=0}\left[af_{i-1}-(a-1)f_{i-m-1}\right]z^i\\F_m(z)&=az-az^{m+1}+azF_m(z)-(a-1)z^{m+1}F_m(z)\\F_m(z)&=\dfrac{az-az^{m+1}}{1+(a-1)z^{m+1}-az}\end{align*}$$
把生成函数的分子提出来,分母先幂级数展开再用牛顿二项式定理展开
$$\begin{align*}F_m(z)&=\dfrac{az-az^{m+1}}{1+(a-1)z^{m+1}-az}\\&=\left(az-az^{m+1}\right)\sum\limits_{k=0}\left[az-(a-1)z^{m+1}\right]^k\\&=\left(az-az^{m+1}\right)\sum\limits_{k=0}\sum\limits_{j=0}^k\binom kj(1-a)^jz^{(m+1)j}a^{k-j}z^{k-j}\end{align*}$$
记$G_m(z)=\sum\limits_{k=0}\sum\limits_{j=0}^k\binom kj(1-a)^ja^{k-j}z^{k+jm}$,则$\left[z^n\right]F_m(z)=a\left[z^{n-1}\right]G_m(z)-a\left[z^{n-m-1}\right]G_m(z)$
现在来求$\left[z^n\right]G_m(z)$,令$k+jm=n$并枚举$j$,同时保证$k\geq j$即$n-jm\geq j$,得到$\left[z^n\right]G_m(z)=\sum\limits_{n-jm\geq j}\binom{n-jm}j(1-a)^ja^{n-jm-j}z^n$
对于所有$m,j$,因为$j(m+1)\leq n$,所以满足条件的$(j,m)$数量是$O\left(n\log_2n\right)$级别的,可以暴力算,时限7s很宽裕
于是就做完了,又学了一个神奇的东西==
#include<stdio.h>
#include<string.h>
typedef long long ll;
const int mod=998244353,ha=19260817;
int max(int a,int b){return a>b?a:b;}
int mul(int a,int b){return a*(ll)b%mod;}
int ad(int a,int b){return(a+b)%mod;}
int de(int a,int b){return(a-b)%mod;}
int pow(int a,int b){
int s=1;
while(b){
if(b&1)s=mul(s,a);
a=mul(a,a);
b>>=1;
}
return s;
}
int M,a,N;
namespace sol1{
struct poly{
int x[410],n;
};
int m;
poly operator+(poly a,poly b){
poly c;
c.n=max(a.n,b.n);
for(int i=0;i<=c.n;i++)c.x[i]=ad(a.x[i],b.x[i]);
return c;
}
poly operator*(poly a,poly b){
poly c;
int i,j;
c.n=a.n+b.n;
memset(c.x,0,sizeof(c.x));
for(i=0;i<=a.n;i++){
if(a.x[i]){
for(j=0;j<=b.n;j++){
if(b.x[j])c.x[i+j]=ad(c.x[i+j],mul(a.x[i],b.x[j]));
}
}
}
return c;
}
poly operator%(poly a,poly b){
int i,j,t,inv;
inv=pow(b.x[b.n],mod-2);
for(i=a.n-b.n;i>=0;i--){
t=mul(a.x[i+b.n],inv);
for(j=0;j<=b.n;j++)a.x[j+i]=de(a.x[j+i],mul(b.x[j],t));
}
while(a.n>0&&a.x[a.n]==0)a.n--;
return a;
}
poly id,one;
poly pow(poly p,int k){
poly res;
res.n=0;
res.x[0]=1;
while(k){
if(k&1)res=res*p%id;
p=p*p%id;
k>>=1;
}
return res;
}
int f[210];
int F(){
int i,s;
id.n=m;
for(i=0;i<m;i++)id.x[i]=1-a;
id.x[m]=1;
one.n=1;
memset(one.x,0,sizeof(one.x));
one.x[1]=1;
one=pow(one,N-1);
for(i=s=0;i<=one.n;i++)s=ad(s,mul(one.x[i],f[i+1]));
return s;
}
void gao(){
int i,s,bas;
f[0]=1;
for(i=1;i<=M;i++)f[i]=mul(f[i-1],a);
s=0;
for(m=bas=1;m<=M;m++){
bas=mul(bas,ha);
s=ad(s,mul(F(),bas));
}
printf("%d",s);
}
}
namespace sol2{
#define maxn 3000010
int an[maxn],a1n[maxn],fac[maxn],rfac[maxn];
int C(int n,int k){return mul(fac[n],mul(rfac[k],rfac[n-k]));}
int G(int m,int n){
int s,j;
for(j=s=0;n-j*m>=j;j++)s=ad(s,mul(mul(C(n-j*m,j),a1n[j]),an[n-j*m-j]));
return s;
}
int F(int m){
return mul(a,de(G(m,N-1),G(m,N-m-1)));
}
void gao(){
int i,res,bas;
an[0]=a1n[0]=fac[0]=1;
for(i=1;i<=N;i++){
an[i]=mul(an[i-1],a);
a1n[i]=mul(a1n[i-1],1-a);
fac[i]=mul(fac[i-1],i);
}
rfac[N]=pow(fac[N],mod-2);
for(i=N;i>0;i--)rfac[i-1]=mul(rfac[i],i);
res=0;
for(i=bas=1;i<N;i++){
bas=mul(bas,ha);
res=ad(res,mul(F(i),bas));
}
res=ad(res,mul(an[N],mul(de(pow(ha,M+1),pow(ha,N)),pow(ha-1,mod-2))));
if(res<0)res+=mod;
printf("%d",res);
}
}
int main(){
scanf("%d%d%d",&M,&a,&N);
if(M>200)
sol2::gao();
else
sol1::gao();
}