$part ~ 1 ~ $多项式除法
01 问题描述
给定一个 \(n\) 次多项式 \(F(x)\) 和一个 \(m\) 次多项式 \(G(x)\) ，请求出多项式 \(Q(x)\), \(R(x)\)，满足以下条件：
- \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 次数小于 \(m\)
- \(F(x) = Q(x) * G(x) + R(x)\)
所有的运算在模 \(998244353\) 意义下进行。
输入格式
第一行两个整数 \(n\)，\(m\)，意义如上。
第二行 \(n+1\) 个整数，从低到高表示 \(F(x)\) 的各个系数。
第三行 \(m+1\) 个整数，从低到高表示 \(G(x)\) 的各个系数。
输出格式
第一行 \(n-m+1\) 个整数，从低到高表示 \(Q(x)\) 的各个系数。
第二行 \(m\) 个整数，从低到高表示 \(R(x)\) 的各个系数。
如果 \(R(x)\) 不足 \(m-1\) 次，多余的项系数补 \(0\)。
样例 #1
样例输入 #1
5 1
1 9 2 6 0 8
1 7
样例输出 #1
237340659 335104102 649004347 448191342 855638018
760903695
提示
对于所有数据，\(1 \le m < n \le 10^5\)，给出的系数均属于 \([0, 998244353) \cap \mathbb{Z}\)。
02 巧妙的转化
首先由题意，我们得到的是如下关系式：
\[F(x)=G(x)Q(x)+R(x) \]先求一发对应的几个多项式的次数范围：
\(F(x)\)次数\(\left[0,n\right]\)，共\(n+1\)项
\(G(x)\)次数\(\left[0,m\right]\)，共\(m+1\)项
\(Q(x)\)次数\(\left[0,n-m\right]\)，共\(n-m+1\)项
\(R(x)\)次数\(\left[0,m-1\right]\)，共\(m\)项
做完多项式求逆后会很自然地想到求逆就是为了除法服务的，于是我们会发现如果没有\(R(x)\)余项的话很显然我们的除法会直接变为一道求逆题. 由于整个式子的多项式范围都在\(x^n\)范围内，如果等式是简单的\(F(x)=G(x)Q(x)\)的话，我们可以化为在\(\pmod{x^{n+1}}\)的意义下求\(G^{-1}(x)\)，而对应的会有\(Q(x)\equiv F(x)Q^{-1}(x)\pmod{x^{n+1}}\)，十分简单
而对于目前我们出现的\(R(x)\)增加的复杂性，我们考虑的思路便是先将\(R(x)\)给消去，这里一个很巧妙的思路便是系数翻转，即
\[F(\frac1x)=G(\frac1x)Q(\frac1x)+R(\frac{1}x) \]此时我们再看看对应的多项式的系数范围：
\(F(\frac1x)\)次数\(\left[-n,0\right]\)，共\(n+1\)项
\(G(\frac1x)\)次数\(\left[-m,0\right]\)，共\(m+1\)项
\(Q(\frac1x)\)次数\(\left[m-n,0\right]\)，共\(n-m+1\)项
\(R(\frac1x)\)次数\(\left[1-m,0\right]\)，共\(m\)项
而接下来为了计算便利我们得把次数都化为自然数，那么两边同时乘以\(x^n\)得
\[x^nF(\frac1x)=x^mG(\frac1x)*x^{n-m}Q(\frac1x)+x^nR(\frac1x) \]看出来这样划分的深意和妙处了吗？接下来我们再看看一些多项式的系数范围：
\(x^nF(\frac1x)\)次数\(\left[0,n\right]\)，共\(n+1\)项
\(x^mG(\frac1x)\)次数\(\left[0,m\right]\)，共\(m+1\)项
\(x^{n-m}Q(\frac1x)\)次数\(\left[0,n-m\right]\)，共\(n-m+1\)项
\(x^nR(\frac1x)\)次数\(\left[n-m+1,n\right]\)，共\(m\)项
接下来我们不妨设\(\begin{cases}F^{'}(x)=x^nF(\frac1x)\\G^{'}(x)=x^mG(\frac1x)\\Q^{'}(x)=x^{n-m}Q(\frac1x)\\R^{'}(x)=x^nR(\frac1x) \end{cases}\)，那么很显然F'、G'、Q'分别对应着F、G、Q系数翻转之后的多项式，而R'则与我们所需要的Q'的系数完全错开，于是我们在\(\pmod{x^{n-m+1}}\)意义下可以刚好完全消去R'而完全不影响Q'的任何一项，即
\[F^{'}(x)\equiv G^{'}(x)Q^{'}(x)\pmod{x^{n-m+1}}\Leftrightarrow Q^{'}(x)\equiv F^{'}(x)G^{'-1}(x)\pmod{x^{n-m+1}} \]于是我们求出Q'之后再将系数翻转便可以得到我们所需要的Q，之后的R便可通过\(R(x)=F(x)-Q(x)G(x)\)求出，而求逆得到\(Q^{'-1}(x)\)的时候一定要注意是对\(\pmod{x^{n-m+1}}\)的逆元
03 代码实现
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int xx=(1<<21)+1,mod=998244353;
const int g1=3,g_1=332748118;
int rev[xx],ite[21],sup;
inline ll power(ll x,int y){ll ans=1;for(;y;(x*=x)%=mod,y>>=1) if(y&1) (ans*=x)%=mod;return ans;}
inline void init(int n)
{
memset(rev,0,sizeof(rev));
for(sup=1;sup<=n;sup<<=1);
for(n=sup>>1;n;n>>=1)
for(register int i=n;i<sup;i+=(n<<1))
for(register int j=0;j<n;++j)
rev[i+j]+=(sup/n/2);
}
inline void ntt(ll *a,const bool op)
{
for(register int i=0;i<sup;++i)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(register int k=1;k<sup;k<<=1)
{
const ll ge=power(op?g1:g_1,(mod-1)/(k<<1));
for(register int i=0;i<sup;i+=(k<<1))
{
ll gj=1;
for(register int j=0;j<k;++j,(gj*=ge)%=mod)
{
const ll l=a[i+j],r=(gj*a[i+j+k])%mod;
a[i+j]=(l+r)%mod;
a[i+j+k]=(l-r+mod)%mod;
}
}
}
}
ll tmp1[xx],tmp2[xx];
inline void inv(ll *f,ll *g,int n)
{
g[0]=power(f[0],mod-2);
for(;n>1;(n+=1)>>=1)
ite[++ite[0]]=n;
ite[ite[0]+1]=1;
for(register int i=ite[0];i;--i)
{
init(ite[i]+ite[1]);
for(register int j=0;j<sup;++j)
{
tmp1[j]=f[j];
if(j<ite[i+1])
tmp2[j]=g[j];
else
tmp2[j]=0;
}
ntt(tmp1,1),ntt(tmp2,1);
for(register int j=0;j<sup;++j)
tmp1[j]=tmp2[j]*((2-tmp1[j]*tmp2[j]%mod+mod)%mod)%mod;
ntt(tmp1,0);
ll inv=power(sup,mod-2);
for(register int j=0;j<ite[1];++j)
g[j]=tmp1[j]*inv%mod;
}
}
ll F[xx],G[xx],Q[xx],R[xx];
ll Fr[xx],Gr[xx],invGr[xx];
int main()
{
int n,m;
scanf("%d%d",&n,&m);
++n,++m;
for(register int i=0;i<n;++i)
scanf("%lld",&F[i]),Q[n-i-1]=F[i];
for(register int i=0;i<m;++i)
scanf("%lld",&G[i]),Gr[m-i-1]=G[i];
for(register int i=n-m+1;i<m;++i)
Gr[i]=0;
inv(Gr,invGr,n-m+1);
init(n<<2);
ntt(Q,1),ntt(invGr,1);
for(register int i=0;i<sup;++i)
Q[i]=Q[i]*invGr[i]%mod;
ntt(Q,0);
ll inv=power(sup,mod-2);
for(register int i=0;i<sup;++i)
Q[i]=Q[i]*inv%mod;
for(register int i=0;i<n-m-i;++i)
swap(Q[i],Q[n-m-i]);
for(register int i=0;i<n-m+1;++i)
printf("%lld ",Q[i]);
puts("");
for(register int i=n-m+1;i<sup;++i)
Q[i]=0;
ntt(F,1),ntt(G,1),ntt(Q,1);
for(register int i=0;i<sup;++i)
R[i]=(F[i]-G[i]*Q[i]%mod+mod)%mod;
ntt(R,0);
for(register int i=0;i<m-1;++i)
printf("%lld ",R[i]*inv%mod);
puts("");
system("pause");
return 0;
}
\(part~ 2 ~\)常系数齐次线性递推
01 问题描述
求一个满足 \(k\) 阶齐次线性递推数列 \({a_i}\) 的第 \(n\) 项，即：
\[a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i} \]输入格式
第一行两个数 \(n\),\(k\)，如题面所述。
第二行 \(k\) 个数，表示 \(f_1 \ f_2 \ \cdots \ f_k\)
第三行 \(k\) 个数，表示 \(a_0 \ a_1 \ \cdots \ a_{k-1}\)
输入格式
一个数，表示 \(a_n \bmod 998244353\) 的值
样例 #1
样例输入 #1
6 4
3 -1 0 4
-2 3 1 5
样例输出 #1
73
提示
对于100%的数据我们有$N \leqslant 10^{9} , K \leqslant 32000 $
保证读入的数字均为 \([-10^9,10^9]\) 内的整数。
02 问题转化
首先我们可以发现这个数列显然是可以通过矩阵构造转化的
首先我们假设如下矩阵：
\[\alpha_i=\begin{pmatrix}a_i\\a_{i+1}\\\vdots\\a_{i+k-1}\end{pmatrix}^T，\beta=\begin{pmatrix}0&0&0&\cdots&0&f_k\\1&0&0&\cdots&0&f_{k-1}\\0&1&0&\cdots&0&f_{k-2}\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&0&f_2 \\0&0&0&\cdots&1&f_{1} \end{pmatrix} \]那么根据矩阵乘法的定义以及我们的递推关系式，我们可以得到
\[\alpha_i\beta=\alpha_{i+1} \]于是
\[\alpha_n=\alpha_0\beta^n \]对于\(\alpha_0\)则是输入数据\(a_{0-(k-1)}\)组成的列向量，而\(\alpha_n\)的第一项则是我们所需要求的\(a_n\)，因此我们的任务可以说直接转化成了求解\(\alpha_0\beta^n\)，其中最核心的一点就是求解\(\beta^n\)
然而对于求解一个矩阵，或者说方阵的n次幂，我们的办法最直接的一项便是计算机直接利用矩阵乘法的定义进行求解，这样求解的复杂度是\(O(k^3n)\). 而由于矩阵乘法的结合律我们可以利用快速幂的思想将这个复杂度化为\(O(k^3\log_2n)\)，然而对这个问题来说这里最大的问题就在于\(k\)这一维也十分难办，不是很好降维，于是我们考虑其他办法
03 特征多项式
我们对于求解一个方阵的n次幂，由线性代数中学到的知识会首先想到将这个方阵给相似对角化然后对应的我们就可以求解其对角矩阵的n次幂之后便可直接求出我们所需要的数值
然而对于这道题而言我们的这种想法也受到了限制，我们假设我们将所需要的转移矩阵给相似对角化如下
\[\beta=P^{-1}\Lambda P \]然而对于我们所目标的矩阵\(\beta^n\)我们可以由\(O(k\log_2n)\)的复杂度求解出\(\Lambda^n\)，然而对于后续的矩阵乘法操作无论是\(P^{-1}\Lambda^nP\)还是\(\alpha\beta^n\)，尽管再怎么优化我们也还是会至少有\(O(k^2)\)的复杂度，对于这道题的范围来说十分不友好，更别说对于一个特殊的\(\beta\)我们能不能在实数范围内顺利将其相似对角化，这是一个十分有难度的问题
然而这个思路也给我们提供了一个十分巧妙的转化想法，无论怎样我们首先的方向是尽量将这个矩阵的特征多项式给求出，只有这样我们才能观察出来相似对角化的可行性
而通过数学归纳法我们最终可以得到我们的特征多项式如下
\[g(\lambda)=(-1)^{k+1}\sum\limits_{i=1}^kf_i\lambda^{k-i}+(-1)^k\lambda^k \]即
\[g(\lambda)=(-1)^{k+1}(-\lambda^k+\sum\limits_{i=1}^kf_i\lambda^{k-i}) \]而由我们线性代数中学到的特征值的相关性质我们可以知道，由于\(g(\lambda)=0\)，那么
\[g(\beta)=(-1)^{k+1}(-\beta^k+\sum\limits_{i=1}^kf_i\beta^{k-i})=0 \]于是通过我们之前得到的多项式除法的相关知识，我们一定能得到
\[\beta^n=Q(\beta)g(\beta)+R(\beta)=R(\beta)\,\,\cdots\cdots ① \]其中\(R(\beta)\)是\(\beta^n\pmod{g(\beta)}\)的余项，可以表示为\(\beta^{0-(k-1)}\)的线性组合
之后最优秀的一步转化，就是在①式两边同时左乘一个\(\alpha_0\)，那么
\[\alpha_0\beta^n=\alpha_0 R(\beta) \]而由我们之前所说的\(R(\beta)=\left<\beta^0,\beta^1,\beta^2,\cdots,\beta^{k-1}\right>\)，可以得到
\[\alpha_n=\left<\alpha\beta^0,\alpha\beta^1,\alpha\beta^2,\cdots,\alpha\beta^{k-1}\right>=\left<\alpha_0,\alpha_1,\alpha_2,\cdots,\alpha_{k-1}\right> \]而由我们向量相等的充要条件，我们可以直接得到
\[a_n=\left<a_0,a_1,a_2,\cdots,a_{k-1}\right> \]且线性表示的系数与我们最开始得到的\(\left<\beta^0,\beta^1,\beta^2,\cdots,\beta^{k-1}\right>\)的系数相同
04 快速幂的思想依旧
而我们使用多项式除法的话复杂度会在\(O((n+k)\log_2(n+k))\)，而这也成为了我们最后一个障碍，也就是过大的n应该怎么处理
而对于这个我们会用到我们的老思路也就是快速幂，因为我们被除的整式的结构十分简单也就是\(x^n\)，那么我们可以轻松地利用快速幂的思路将其转化为\(x^1,x^2,x^4,\cdots,x^{2^m}\)的乘积，最后面我们可以大大缩短我们所需要进行多项式除法的被除式的长度，可以保持每次的多项式除法中的被除式的长度不会超过2k，那么一次多项式除法的复杂度不会超过\(\Theta(3k\log_2 3k)\)，而快速幂的乘法次数不会超过\(\Theta(2\log_2n)\)，那么最终的时间复杂度化为\(O(k\log_2k\log_2n)\)，相较于\(O(k^3\log_n)\)大大简化
当然k的大小还是制约我们算法的最大限制，但能将其化为线性复杂度已经能够解决极大一部分的问题了
最后面还有一个小的优化细节，那就是我们每次多项式除法的除式都是相同的，那么我们除式对应的翻转逆也可以预处理求出，这样可以大大节省我们用在多项式除法中的费用
05 程序实现
最后的程序实现一定要注意由于我们需要进行多次多项式除法，有些数组需要及时清理，不然无法保证多项式除法多次执行的正确性
#include<bits/stdc++.h>
using namespace std;
#define cl(a) memset(a,0,sizeof(a))
typedef long long ll;
const int xx=(1<<17)+10,mod=998244353;
const int g1=3,_g=332748118;
int rev[xx],ite[21],sup;
inline int read()
{
char ch=getchar();
int x=0,f=1;
for(;!isalnum(ch);ch=getchar()) if(ch=='-') f=-1;
for(;isalnum(ch);ch=getchar()) x=(x<<3)+(x<<1)+ch-'0';
return f*x;
}
inline ll power(ll x,int y)
{
ll ans=1;
for(;y;(x*=x)%=mod,y>>=1)
if(y&1)
(ans*=x)%=mod;
return ans;
}
inline void init(int n)
{
memset(rev,0,sizeof(rev));
for(sup=1;sup<=n;sup<<=1);
for(register int k=sup>>1;k;k>>=1)
{
n=(sup/k)>>1;
for(register int i=k;i<sup;i+=(k<<1))
for(register int j=0;j<k;++j)
rev[i+j]+=n;
}
}
inline void ntt(ll *a,const bool op)
{
for(register int i=0;i<sup;++i)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(register int k=1;k<sup;k<<=1)
{
const ll ge=power(op?g1:_g,(mod-1)/(k<<1));
for(register ll i=0;i<sup;i+=(k<<1))
{
ll gj=1;
for(register int j=0;j<k;++j,(gj*=ge)%=mod)
{
const ll l=a[i+j],r=(gj*a[i+j+k])%mod;
a[i+j]=(l+r)%mod;
a[i+j+k]=(l-r+mod)%mod;
}
}
}
}
inline void finish(ll *a)
{
ntt(a,0);
ll invsup=power(sup,mod-2);
for(register int i=0;i<sup;++i)
(a[i]*=invsup)%=mod;
}
ll tmp1[xx],tmp2[xx];
inline void inv(ll *a,ll *b,int n)
{
cl(ite);
b[0]=power(a[0],mod-2);
for(;n>1;n=(n+1)>>1)
ite[++ite[0]]=n;
ite[ite[0]+1]=1;
for(register int i=ite[0];i;--i)
{
init(ite[i]<<1);
for(register int j=0;j<sup;++j)
{
tmp1[j]=j<ite[i]?a[j]:0;
tmp2[j]=j<ite[i+1]?b[j]:0;
}
ntt(tmp1,1),ntt(tmp2,1);
for(register int j=0;j<sup;++j)
tmp1[j]=tmp2[j]*((2-tmp1[j]*tmp2[j]%mod+mod)%mod)%mod;
finish(tmp1);
for(register int j=0;j<ite[i];++j)
b[j]=tmp1[j];
}
}
ll tmp3[xx],tmp4[xx],tmp5[xx],tmp6[xx];
inline void getmod(ll *a,int n,ll *b,int m)
{
if(n<m) return;
cl(tmp3);cl(tmp4);cl(tmp5);
for(register int i=0;i<n;++i)
tmp3[i]=a[n-i-1];
for(register int i=0;i<n-m+1;++i)
tmp5[i]=tmp6[i];
init(2*n-m+1);
ntt(tmp3,1),ntt(tmp5,1);
for(register int i=0;i<sup;++i)
tmp4[i]=(tmp3[i]*tmp5[i])%mod;
finish(tmp4);
for(register int i=0;i<sup;++i)
{
if(i<m) tmp5[i]=b[i];
else tmp5[i]=0;
if(i<n-m+1) tmp3[i]=tmp4[n-m-i];
else tmp3[i]=0;
}
ntt(tmp3,1),ntt(a,1),ntt(tmp5,1);
for(register int i=0;i<sup;++i)
a[i]=(a[i]-tmp5[i]*tmp3[i]%mod+mod)%mod;
finish(a);
for(register int i=m-1;i<sup;++i)
a[i]=0;
}
ll Ans[xx],kid[xx];
inline void PoW(ll *a,int n,ll *Mod,int m)
{
Ans[0]=1;
for(register int x=2,y=1;n;n>>=1)
{
if(n&1)
{
init(y+=x);
for(register int i=0;i<sup;++i)
kid[i]=a[i];
ntt(Ans,1),ntt(kid,1);
for(register int i=0;i<sup;++i)
(Ans[i]*=kid[i])%=mod;
finish(Ans);
getmod(Ans,y,Mod,m);
if(y>m) y=m-1;
}
init(x+=x);
ntt(a,1);
for(register int i=0;i<sup;++i)
(a[i]*=a[i])%=mod;
finish(a);
getmod(a,x,Mod,m);
if(x>m) x=m-1;
}
}
ll f[xx>>1],a[xx>>1],tmp[xx];
int main()
{
int n,k;
n=read();f[k=read()]=mod-1;
for(register int i=k-1;~i;--i)
{
f[i]=read();
if(f[i]<0)
f[i]+=mod;
}
for(register int i=0;i<k+1;++i)
tmp4[i]=f[k-i];
inv(tmp4,tmp6,k);
for(register int i=0;i<k;++i)
{
a[i]=read();
if(a[i]<0)
a[i]+=mod;
}
tmp[1]=1;
PoW(tmp,n,f,k+1);
ll ans=0;
for(register int i=0;i<k;++i)
(ans+=(Ans[i]*a[i]%mod))%=mod;
printf("%lld\n",ans);
return 0;
}