多项式小记
多项式小记
基本介绍
多项式，听起来非常抽象的东西，的确，这个东西是一个省选的知识点，所以它确实不简单。但是随着算法学习的深入，它也显得非常重要。其主要功能，或者说其中包含算法的主要功能主要是更高效率的进行一些多项式的计算，其中的一些思想可以说是极为优美。
FFT
前置芝士
- \(F(x)\) 表示多项式 \(F\)，例如 \(A(x)\) 可以为 \(x^2+3x+2\)。
- 系数表示法：一般来说，我们表示多项式用的都是系数表示法。比如对于上面的 \(A(x)\)，我们可以用数组 \([1,3,2]\) 来表示它每一项的系数，简单明了。
- 系数表示法计算乘积：在系数表示发的情况下，如何计算多项式的乘积？非常简单，若 \(C(x)=A(x)\times B(x)\) 则 \(C[k]=\sum_{i=0}^{k}A[i]\times B[k-i]\)，假设 \(A\) 有 \(n\) 项，\(B\) 有 \(m\) 项，则计算的时间复杂度就是 \(\mathcal{O}(nm)\) 的
DFT(IDFT)思想
为了加速计算过程，我们首先从表示方法考虑。有没有其他可以表示多项式的方法？
是有的。考虑一个一次函数，我们知道其图像一定是一个直线，我们知道两点确定一个直线，因此我们可以用两个点来表示一个一次函数。
把这个结论推广（可以去思考更详细的证明），我们可以知道 \((d+1)\) 个不同点可以确定一个 \(d\) 次的多项式。
所以对于 \(F(x)\) 我们可以用 \((x_0,F(x_0)),...,(x_n,F(x_n))\) 来表示
这样子我们会发现，我们再次计算 \(A(x)\times B(x)\) 的时候，我们知道 \(C\) 至多有 \(mn\) 项，因此我们在 \(A,B\) 中用 \(mn+1\) 个不同的值计算出其对应的函数值，相乘，就可以得到 \(C\) 里面所有这些值对应的函数值，之后可以反推出来 \(C\) 的表达式。
但是问题是，如何把系数表达式，化为点表达式？
如果随机取点，我们需要取 \(\ge n+1\) 个点，会发现复杂度是\(O(nd)\)回到原点。
系数 \(\to\) 点，DFT
思考\(F(x)=x^2\)，我们发现其一个性质 \(F(-x)=F(x)\)，这样我们我们是不是只需要取 \(\dfrac{n}{2}\) 个点？
思考 \(F(x)=x^3\)，有 \(F(-x)=-F(x)\)，是不是也是只需要 \(\dfrac{n}{2}\) 个点？
浅浅推广一些，我们假设\(F_e(x),F_o(x)\)分别表示一个函数当中所以偶数和奇数项之和，他们的系数依赖于\(F(x)\)，比如说 \(F_e(x)=F[0]+F[2]x+...+F[n-2]x^{n/2-1}\)。
我们会有以下两个式子：\(F(x)=F_e(x^2)+xF_o(x^2)\) 和 \(F(-x)=F_e(x^2)-xF_o(x^2)\)。
这样做的目的是为了分治，即我们想用尽可求少的点具体的值，来得到尽可能多的点——少算、多点。上面式子即如果知道了 \(n\) 个 \(F_e(x),F_o(x)\) 的值那么就可以得到 \(2n\) 个 \(F(x)\) 的值。
但是我们分成一半之后，会发现 \([x_1^2,...,x_{n/2}^{2}]\) 不是正负匹配的了，即我们在分治的过程中计算 \(F_e(x)\) 的值时，被强制使用了平方，即正数，用不了刚才这个套路。那么如何做到可以继续正负匹配？
这里我们就考虑到虚数了！
令 \(\omega_{i}^j\) 表示 \(x^i=1\) 的第 \(j\) 个解。即在复平面上的单位圆。
对于它，有几个点：
- \(\omega_{n}^{k}\) 可以理解为把 \((1,0)\) 绕原点旋转 \(\dfrac{k}{n}\) 个圆周得到的点
- \(\omega_{n}^{k}=\omega_{pn}^{pk}\)，也就是把这个圆分成一定份数去某一份，将这个过程翻个倍，不会影响
- \(\omega_{n}^{j}\times \omega_{n}^{k}=\omega_{n}^{j+k}\)
- 若 \(n\) 是 \(2\) 的倍数，\(\omega_{n}^{k+n/2}=-\omega_{n}^{k}\)
我们代入一下分析：
Case 1:
\(F(x)=F_e(x^2)+xF_o(x^2)\)
经过一些化简可以得到：
\(F(\omega_{n}^{k})=F_e(\omega_{n}^{2k})+\omega_{n}^{k}F_o(\omega_{n}^{2k})\)
Case 2:
\(F(-x)=F_e(x^2)-xF_o(x^2)\)
经过一些化简可以得到：
\(F(-\omega_{n}^{k})=F_e(\omega_{n}^{2k})+\omega_{n}^{k}F_o({\omega_{n}^{2k}})\)
注：这一部分化简技巧性不大，暴力带入 \(\omega_{i}^j\) 即可，然后通过上面的几个性质带入化简。
更好的操作
我们把上面的式子转换一下，变成这样子会好做一点（因为上面的 \(2k\) 可能会超过 \(n\)，代码做不到极简，当上面的写法也是完全没有问题的）
\(F(ω^{k+n/2}_n)=F_e(ω^{k}_{n/2})-ω^k_nF_o(ω^{k}_{n/2})\)
\(F(ω^k_n)=F_e(ω^k_{n/2})+ω^k_nF_o(ω^k_{n/2})\)
这意味着什么？这意味着知道后面两个值的情况我们可以 \(O(n)\) 求出来 \(F\) 的值了！
最后，我们再考虑 \(\omega\) 怎么求。这里运用一点非常非常基础的三角函数的值就知道 \(\omega_n^1=(\cos \dfrac{2\pi}{n},\sin \dfrac{2\pi}{n})\)（单位圆上转一个特定角度）。于是乎，经过我们这样简单的操作，我们就做到了能够 \(O(nlogn)\) 做到系数\(\to\)点
系数 \(\gets\) 点 IDFT
最后反推的过程，这里有一个结论就是把上面的 \(\sin \dfrac{2\pi}{n}\) 换成 \(-\sin \dfrac{2\pi}{n}\)，最后结果除 \(n\) 就可以啦！这里也可以通过推理\(\omega\)得出来，这里不太多赘述，只要理解上面的部分，下面也就基本可以理解了。
补，考虑证明：
设：
需证明（即类似反演的东西）：
右边带入 \(G_i\) 得：
考虑 \(\omega_n^{i(j-k)}\) 的一些性质。令 \(v=j-k\)，有对于单个 \(j\) 有：
注意到第一个 \(\sum\) 可以是等比求和，公比是 \(\omega_n^v\)，根据公式得 \(\sum\omega_n^{iv}=\dfrac{\omega_n^0(1-(\omega_n^v)^n)}{1-\omega_n^v}=\dfrac{1-(\omega_n^n)^v}{1-\omega_n^v}=0,\omega_n^v\neq 1\)。当 \(\omega_n^v=0\) 时，常数列和为 \(n\omega_n^0=n\)。综合起来，前面贡献为 \(0\)，因此答案等于后半部分，得证。
代码实现
跟一些模板一样，FFT 的实现和理论还是有不小的距离的。结果上来说，只记住上面的结论式子就能够应对绝大部分需要 FFT 的多项式题了，但是模板需要全部记忆或者理解。
难点在于各种优化，由简到难（没看懂参考代码）：
DFT \(\rightleftarrows\) IDFT
可以封装成一个函数，如上所述，差距只有一个负号。
“分”
按照推理的逻辑，我们会按照 \(\bmod~2=0/1\) 来分成两个部分，开两个数组，然后递归计算。这不禁让我们想，能否开始直接算出一个对于下标 \(i\) 最后会变化到的位置 \(t_i\)，这样分治的过程中不用实际的“分”？那么需要找到 \(i\rightleftarrows t_i\) 的关系了。
直接上结论，这个过程为“蝴蝶变换”，即满足 \(n\) 长度是 \(2^k\) 的时候 \(i,t_i\) 互为反转二进制之后的结果（举例：\((3)_{10}=(011)_2\leftrightarrows (6)_{10}=(110)_2\)）。 考虑快速进行这个过程，设 \((x)_2=(V)_2+(a)_2\) 其中 \(a\) 表示 \(x\) 的第一位。那么 \(t_x=a\times 2^{len}+t_V\)。代码：
to[x]=((i&1)?(n>>1):0)|(to[(x>>1)]>>1);
“治”
回顾：
注意到直接在原数组上修改即可，不需要单独拷贝数组。
for(int k=j;k<j+len;k++){
Cp tmp=now*f[len+k];
f[len+k]=f[k]-tmp;
f[k]=f[k]+tmp;
}
三角函数 + 递推
三角函数常数较大，改为递推的形式。因为确定了 \(n\)，就知道了 \(\omega_n^1=(\cos \dfrac{2\pi}{n},\sin \dfrac{2\pi}{n})\)，只需要求 \(\log\) 次。
这些优化合起来就是代码实现的 FFT：
#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
using namespace std;
const double pi=acos(-1);
const int maxn=5e6+10;
struct Cp{//自定义复数类
Cp (double xx=0,double yy=0){x=xx,y=yy;}
double x,y;
Cp operator + (Cp const &B) const{return Cp(x+B.x,y+B.y);}
Cp operator - (Cp const &B) const{return Cp(x-B.x,y-B.y);}
Cp operator * (Cp const &B) const{return Cp(x*B.x-y*B.y,x*B.y+y*B.x);}
}F[maxn],G[maxn];
int n,m;
int to[maxn];
void FFT(Cp *f, bool flag){//flag=1:DFT
for(int i=0;i<n;i++)if(i<to[i])swap(f[i],f[to[i]]);
for(int i=2;i<=n;i<<=1){
int len=i>>1;
//当前的长度，len表示长度的一半
Cp bas(cos(2*pi/i),sin(2*pi/i));
//固定了长度之后单位自然也固定了
if(!flag)bas.y*=-1;
for(int j=0;j<n;j+=i){//枚举起点
Cp now(1,0);
for(int k=j;k<j+len;k++){//枚举具体位置
//直接在数组上操作，动手写出式子写出来然后画一下可以立刻明白
Cp tmp=now*f[len+k];
f[len+k]=f[k]-tmp;
f[k]=f[k]+tmp;
now=now*bas;
//omega_n^i->omega_n^(i+1)
}
}
}
}
//蝴蝶变化
void initto(int lim){for(int i=0;i<lim;i++)to[i]=((i&1)?(n>>1):0)|(to[i>>1]>>1);}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++)scanf("%lf",&F[i].x);
for(int i=0;i<=m;i++)scanf("%lf",&G[i].x);
for(m+=n,n=1;n<=m;n<<=1);//补齐
initto(n);
FFT(F,1),FFT(G,1);
for(int i=0;i<n;i++)F[i]=F[i]*G[i];
FFT(F,0);
for(int i=0;i<=m;i++){
printf("%d ",(int)(F[i].x/n+0.49));//取到最近的整数
}
}
感觉前面的代码理解铺垫应该已经够详细了。
终于把 FFT 更新了，也不知道拖了多久。
NTT
FFT 虽然可以快速计算多项式，但是有一些天然劣势。因为其用的是 double，因此当数字很大的时候精度会不够。而在 OI 当中数字过大的时候一般采用取模的方式解决。然而数学家已经证明了在 \(\mathbb{C}\) 中只有单位根满足要求的一类数。那么就需要考虑，有没有可以在模意义下单位根的替代品？
原根
原根涉及到很多的数学证明，由于本人水平不足加上时间问题，这里省略大量证明，毕竟作为 OIer 其实有一些证明也不需要太明白。
定义：\(a\) 在模 \(p\) 意义下的阶为最小的 \(k\) 满足 \(a^k\equiv 1\bmod p\)，记为 \(\delta_p(a)\)。即最小的循环节。
定义：若有 \(g\) 满足 \(\delta_p(g)=\varphi(p)\) 那么 \(g\) 就是 \(p\) 的一个原根。这里的 \(\varphi(p)\) 就是 \(p\) 意义下阶的上界，根据欧拉定理，这个循环最长就是 \(\varphi(p)\)。
定理：只有 \(2,4,p^k,2p^k\) （其中 \(p\) 为奇素数）这一类数有原根。
定理：若 \(n\) 有原根，那么其是 \(\mathcal{O}(n^\frac{1}{4})\) 级别的。
简述一下如何找 \(n\) 的所有原根。有定理，对于 \(\gcd(a,p)=1,a^k\equiv 1\bmod p\) 有 \(k\) 为 \(\varphi(p)\) 的因数。又有，如果 \(a^{ck}\equiv 1 \bmod p\) 那么 \(a^k\equiv 1\bmod p\)。因此只需要枚举 \(p\)，之后检查是否存在 \(\dfrac{\varphi(p)}{p_i}\) 等于 \(1\) 即可，其中 \(p_i\) 是 \(\varphi(p)\) 的因数（请忽略变量重名）。
结论：满足 \(n|(p-1)\) 的前提下，\(g^{\frac{p-1}{n}}\) 是一个可以作为单位根的数。证明带入单位根的各种要求即可。
基本的，\(998244353\) 的原根为 \(g=3\)，\(10^9+7\) 的原根为 \(g=5\)。
分治 FFT / NTT
好像也并没有真正的算法知识，是一个技巧性的东西。看例题：分治 FFT
给定多项式 \(G\)，求多项式 \(F\) 满足：
考虑进行类似 cdq 分治一样的东西，对于一个区间分治之后，计算左边对于右边的贡献。令当前计算的区间为 \([l,r]\) 那么对于右边的 \(w_x\) 有 \(w_x=\sum_{i=l}^{mid}F_{i}G_{x-i}\)，由于左边的答案已经求出来了，通过卷积计算这一部分即可。这样可以做到 \(\mathcal{O}(n\log^2 n)\)。
这一题还有另解，见下。
求逆
引入：界，即一个多项式 \(\bmod~x^k\) 的结果，即我们只想提取前 \(k\) 项。不过这个东西好像也只是一个形式化的表示罢了。
对于多项式求逆，即对于 \(F(x)\) 求 \(G(x)\) 满足 \(F(x)G(x)=1\)
考虑倍增，已知 \(R(x)=F^{-1}(x)\bmod x^{\frac{n}{2}}\)，欲求 \(T(x)=F^{-1}(x)\bmod x^n\)
考虑 \(T(x)-R(x)\equiv 0 \bmod x^{\frac{n}{2}}\)。因为余数是 \(0\)，那么当模数平方的时候前面跟着平方等式仍成立：
清醒点，这是同余意义下，用什么求根公式呢。利用逆的性质，同乘 \(F(x)\)：
注意因为界变成了 \(n\)，所以抵消的只能是 \(T(x)\)。最终有
Poly inv(Poly A){
int sizA=A.siz(),n=1;
for(;n<sizA;n<<=1);
Poly T,R,sav;
T.resiz(n<<1),R.resiz(n<<1),sav.resiz(n<<1),A.resiz(n<<1);
T.a[0]=qpow(A.a[0],mod-2);
for(int len=2;len<=n;len<<=1){
for(int i=0;i<(len>>1);i++)R.a[i]=(T.a[i]<<1)%mod;
sav.cpy(A,len);
NTT(T.a,len<<1,1),px(T,T,len<<1);
NTT(sav.a,len<<1,n),px(T,sav,len<<1);
NTT(T.a,len<<1,0),T.cl(len,len<<1);
for(int i=0;i<len;i++)T.a[i]=(R.a[i]-T.a[i]+mod)%mod;
}
T.resiz(sizA),A.resiz(sizA);
return T;
}
附：对于上面那题，考虑两个多项式的关系可以表示为
化简就有：
求逆直接秒了，而且比分治少一个 \(\log\)。
求导 + 积分
和正常求导和积分一模一样：
都是可以 \(\mathcal{O}(n)\) 直接写出的。
Poly dr(Poly A){
Poly res;
res.resiz(A.siz()-1);
for(int i=0;i<res.siz();i++)res.a[i]=A.a[i+1]*(i+1)%mod;
return res;
}
Poly itg(Poly A){
vector<int>pinv;
Poly res;
pinv.resize(A.siz()+1),res.resiz(A.siz()+1);
pinv[1]=1,pinv[0]=0;
for(int i=2;i<pinv.size();i++)pinv[i]=pinv[mod%i]*(mod-mod/i)%mod;
for(int i=1;i<res.siz();i++)res.a[i]=A.a[i-1]*pinv[i]%mod;
res.a[0]=0;
return res;
}
牛顿迭代（复合）
定义复合多项式：
重要定理
若 \(G(F(x))=0\) 有
这个结论就是牛迭的关键，证明大概是在 \(F_*(x)\) 处进行泰勒展开，没有仔细研究。
补：尝试证明。
有泰勒展开：
即 \(F(x)\) 在 \(x=a\) 位置的泰勒展开。
考虑 \(G(F(x))\) 在 \(G(F_*(x))\) 位置的展开，有：
考虑 \(F(x)-F_*(x)\) 是什么东西，因为 \(F_*(x)\equiv F(x)\bmod x^\frac{n}{2}\)，也就是说他们 \([0,x^\frac{n}{2})\) 部分的系数是一样的，最低项至少是 \(x^\frac{n}{2}\) 级别的。又因为要求的是 \(\bmod x^n\) 级别下的答案，因此 \((F(x)-F_*(x))^k,k\ge 2\) 都是 \(=0\) 的，得到：
然后移项加上 \(G(F(x))=0\) 的条件就得到了：
这个公式非常厉害，如上面的多项式求逆的推导，对于 \(A(x)\) 求 \(B(x)\) 满足 \(A(x)B(x)=1\)。
改一下，变成 \(A(x)B(x)-1=0\)，令 \(G(B(x))=A(x)B(x)-1\)，则 \(G'(B(x))=A(x)\) 注意是对着 \(B(x)\) 导的。下面令 \(F_*(x)\) 为 \(F(x)\) 在 \(\bmod~2^\frac{n}{2}\) 的结果。
\(B(x)=B_*(x)-\frac{1}{A(x)}(A(x)B_*(x)-1)\)，注意到因为导的对象是 \(B_*(x)\) 那么前面是 \(\frac{1}{A_*(x)}\)，可以用 \(B_*(x)\) 替换：
比较之前的 \(T,R\) 式子会发现结果是吻合的。
开根
我们尝试用一下刚才说的牛迭。
对于 \(A(x)\) 求 \(B^2(x)\equiv A(x) \bmod p\)。
有 \(B^2(x)-A(x)\equiv 0\)，令 \(G(B(x))=B^2(x)-A(x)\)，则 \(G'(B(x))=2B(x)\)，带入公式：
对于 \(B_0\) 结果是二次剩余，出门左转到……哦我好像没有把取模方程放到博客上，谔谔。
又要定义复数域，感觉回到了 FFT。不过注意 Cipolla 的复数运算为 \(i^2=a^2-n\)，因为那不是真正的复数。
Poly sq(Poly A){
int sizA=A.siz(),n=1;
for(;n<sizA;n<<=1);
Poly T,R;
T.resiz(n<<1),R.resiz(n<<1),A.resiz(n<<1);
T.a[0]=Cipolla(A.a[0]);
if(mod-T.a[0]<T.a[0])T.a[0]=mod-T.a[0];
for(int len=2;len<=n;len<<=1){
for(int i=0;i<(len>>1);i++)R.a[i]=(T.a[i]<<1)%mod;
R.cl(len>>1,len);
R=inv(R);
NTT(T.a,len,1),px(T,T,len),NTT(T.a,len,0);
for(int i=0;i<len;i++)T.a[i]=(T.a[i]+A.a[i])%mod;
T.resiz(len<<1),R.resiz(len<<1);
T=mut(T,R);
T.cl(len,len<<1);
}
T.resiz(sizA),A.resiz(sizA);
return T;
}
带余除法
多项式求逆好比计算 \(7\div 3=\dfrac{7}{3}\) 那么带余除法就是计算 \(7\div 3=2\dots\dots 1\)。具体地，给定 \(n\) 次多项式 \(A(x)\) 和 \(m\) 次多项式 \(B(x)\) 求 \(A(x)=Q(x)*B(x)+R(x)\) 其中满足 \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 的次数小于 \(m\)。
这一部分好像和牛顿迭代没有太大的关系，但是这个推导也是比较巧妙的。
先明确一下思路，我们需要想办法通过一些操作把 \(R(x)\) 消掉，然后用多项式求逆来做。
对于一个 \(n\) 次多项式，定义操作 \(\mathcal{T}\)：
发现其实相当于把系数全部反过来了。
对于要求的式子换元 + 同乘系数：
写成反转的形式有：\(A^{\mathcal{T}}(x)=Q^{\mathcal{T}}(x)B^{\mathcal{T}}(x)+x^{n-m+1}R^{\mathcal{T}}(x)\)
看起来没有实质上的改变，但是注意到 \(Q^{\mathcal{T}}(x)\) 是 \(n-m\) 次的，那么直接对于这个式子 \(\bmod x^{n-m+1}\) 就有：
通过多项式求逆得到 \(Q^{\mathcal{T}}(x)\) 之后得到 \(R(x)\) 就是平凡的了。
pair<Poly,Poly>mdiv(Poly A,Poly B){
Poly AT,BT,Q,tmp,R;
R.resiz(B.siz()-1);
int L=A.siz()-B.siz()+1;
AT.resiz(A.siz()),BT.resiz(B.siz());
AT.cpy(A,AT.siz()),BT.cpy(B,BT.siz());
AT.rev(),BT.rev(),AT.resiz(L),BT.resiz(L);
BT=inv(BT),Q=mut(AT,BT);
Q.resiz(L),Q.rev();
tmp=mut(B,Q);
for(int i=0;i<B.siz()-1;i++){
R.a[i]=A.a[i]-tmp.a[i];
if(R.a[i]<0)R.a[i]+=mod;
}
return make_pair(Q,R);
}
求 \(\ln\)
给定多项式 \(A(x)\)，求 \(B(x)=\ln A(x)\)。
首先先明确一些东西，有：
证明直接对两边求导得到导函数相等，加上 \(x=0\) 时两边相等即可证明。
这个对于多项式同样有用，再加上换元就有：
当然，这个式子对于直接做 \(\ln\) 是没什么帮助的，因为复杂度显然会太大。考虑利用一些性质，有 \(\ln'(x)=\dfrac{1}{x}\)，那么可以对于 \((\ln(A(x)))\) 求导，就有：
这样就可以算出来 \(B'(x)\bmod x^{n-1}\)，那么对其积分就可以得到 \(B(x)\bmod n\)。
注意到这样做的 \(\ln\) 会有一个前置条件，就是常数项为 \(1\)。否则需要提取前面的项将常熟化为 \(1\)。
Poly ln(Poly A){
Poly res,B=A;
B=inv(A),A=dr(A);
res=mut(A,B);
res=itg(res);
res.resiz(A.siz()+1);
return res;
}
求 \(\exp\)
给定 \(A(x)\) 求 \(B(x)=\exp A(x)\)。
同样我们只考虑首项为 \(0\) 的情况。
继续使用牛顿迭代，第一步先转化一下，有 \(\ln B(x) = A(x)\)，令 \(G(B(x))=\ln B(x)-A(x)\)，有 \(G'(B(x))=B(x)^{-1}\)，牛顿迭代：
Poly exp(Poly A){
Poly R,T;
int n=1,siz=A.siz();
for(;n<A.siz();n<<=1);
A.resiz(n<<1),T.resiz(n<<1),R.resiz(n<<1),T.a[0]=1;
for(int len=2;len<=n;len<<=1){
R.cpy(T,len>>1),R=ln(R);
for(int i=0;i<len;i++)R.a[i]=(A.a[i]-R.a[i]+mod)%mod;
R.a[0]=(R.a[0]+1)%mod;
T.resiz(len<<1),R.resiz(len<<1);
T=mut(T,R);
T.cl(len,len<<1);
}
T.resiz(siz);
return T;
}
快速幂
给定 \(A(x)\) 求 \(A^k(x)\)。
比较直观的是就像数字快速幂一样做，但是乘法的代价是 \(\mathcal{O}(n\log n)\) 的，因此最后的复杂度是 \(\mathcal{O}(n\log^2 n)\) 的。
注意到有 \(A^k(x)=e^{\ln A(x)\times k}\)。直接做即可。但是这样做其实常数也是挺大的。
这三个部分本人都还没有写常数项是一般情况下的做法，暂时放下，放一个完整的板子。常数不是很优秀，简介程度也说不上好，但是能够应对大部分情况了。
任意模数 NTT
做 NTT 的时候大部分情况遇到的模数就是 \(998244353\)，这是一个非常和谐的数字，因为 \(998244353=2^{23}\times 7\times 17+1\)，而 \(2^{23}=8\times 10^6\)，因此几乎所有情况下我们构造出来的 \(n\) 都能有 \(n|(p-1)\)。但是，如果模数换成同样常见但是其实并不和谐的 \(10^9+7=2\times 500000003+1\)，就会比较麻烦。
解决这种问题比较好些的方法就是三模数 NTT。需要知道的是中国剩余定理 CRT 不知道 Ex 版没关系，普通版本结论知道即可。
一般题目中模数会 \(\le 10^9\)，\(n\le 10^5\)，也就是说最大的系数应该是 \(10^{23}\) 的范围，根据中国剩余定理，如果用三个模数 \(m_1,m_2,m_3\) 满足他们 \(\prod\ge 10^{23}\)，分别求出来在他们意义下的解，那么最终结果在 \(\bmod \prod m\) 意义下有唯一解。\(m_1,m_2,m_3\) 我们一般用 \(998244353,1004535809,469762049\)，因为这三个数都包含大量 \(2\) 的因子。
但是要直接做 CRT 会爆 longlong 范围，考虑怎么解决这个问题。令三个解分别是 \(a_1,a_2,a_3\)，先算出来前两个情况，能够得到：
其中逆元都是 \(\bmod m_1m_2\) 意义下的。
可以建立如下关系：
在 \(\bmod m_3\) 意义下有 \(k\equiv(a_3-r)(m_1m_2)^{-1}\bmod m_3\)
全家桶
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod=998244353,_G=3,maxn=5e5+10;
int qpow(int x,int y){
int res=1;
for(;y;y>>=1){
if(y&1)res=res*x%mod;
x=x*x%mod;
}
return res;
}
const int invG=qpow(_G,mod-2);
struct cmx{
int x,y,w;
cmx operator * (cmx const & a)const{
cmx res;
res.x=((x*a.x)%mod+(y*a.y)%mod*w%mod)%mod;
res.y=((x*a.y)%mod+(y*a.x)%mod+mod)%mod;
res.w=w;
return res;
}
};
int chk_rem(int x){
int res=qpow(x,(mod-1)/2);
if(res==mod-1)return -1;
return 1;
}
cmx cmxqpow(cmx x,int y){
cmx res=(cmx){1,0,x.w};
for(;y;y>>=1){
if(y&1)res=res*x;
x=x*x;
}
return res;
}
int Cipolla(int X){
if(X==0)return 0;
if(chk_rem(X)==-1)assert(0);//no answer
int a;
cmx ans;
while(1){
a=rand()%mod;
ans.w=(a*a%mod-X+mod)%mod;
if(chk_rem(ans.w)==-1)break;
}
ans.x=a,ans.y=1;
return cmxqpow(ans,(mod+1)/2).x;
}
void NTT(vector<int> &g,int n,int flag){//flag=0:inverse
vector<int>w,f,to;
w.resize(n),f.resize(n),to.resize(n);
for(int i=0;i<n;i++)to[i]=((i&1)?(n>>1):0)|(to[i>>1]>>1);
for(int i=0;i<n;i++)f[i]=g[to[i]];
for(int T=2;T<=n;T<<=1){
int bas=flag?qpow(_G,(mod-1)/T):qpow(invG,(mod-1)/T),len=T>>1;
w[0]=1;
for(int i=1;i<T;i++)w[i]=w[i-1]*bas%mod;
for(int j=0;j<n;j+=T){
for(int k=j;k<j+len;k++){
int tmp=w[k-j]*f[len+k]%mod;
f[len+k]=f[k]-tmp;
f[k]=f[k]+tmp;
if(f[len+k]<0)f[len+k]+=mod;
if(f[k]>=mod)f[k]-=mod;
}
}
}
if(!flag){
int ivn=qpow(n,mod-2);
for(int i=0;i<n;i++)g[i]=f[i]*ivn%mod;
} else for(int i=0;i<n;i++)g[i]=f[i];
}
struct Poly{
vector<int>a;
int siz(){return a.size();}
void resiz(int x){a.resize(x);}
void cpy(Poly V,int len){for(int i=0;i<len;i++)a[i]=V.a[i];}
void cl(int l,int r){for(int i=l;i<r;i++)a[i]=0;}
void rev(){reverse(a.begin(),a.end());}
void pr(){
for(int i=0;i<a.size();i++)cout<<a[i]<<" ";
cout<<"\n";
}
};
void px(Poly &X,Poly Y,int len){
for(int i=0;i<len;i++)X.a[i]=X.a[i]*Y.a[i]%mod;
}
Poly mut(Poly X,Poly Y){
Poly res;
int sizX=X.siz(),sizY=Y.siz(),n=1,tsiz=sizX+sizY;
for(;n<sizX+sizY;n<<=1);
X.resiz(n),Y.resiz(n),res.resiz(n),res.cpy(X,n);
NTT(res.a,n,1),NTT(Y.a,n,1);
px(res,Y,n),NTT(res.a,n,0);
res.resiz(tsiz);
return res;
}
Poly inv(Poly A){
int sizA=A.siz(),n=1;
for(;n<sizA;n<<=1);
Poly T,R,sav;
T.resiz(n<<1),R.resiz(n<<1),sav.resiz(n<<1),A.resiz(n<<1);
T.a[0]=qpow(A.a[0],mod-2);
for(int len=2;len<=n;len<<=1){
for(int i=0;i<(len>>1);i++)R.a[i]=(T.a[i]<<1)%mod;
sav.cpy(A,len);
NTT(T.a,len<<1,1),px(T,T,len<<1);
NTT(sav.a,len<<1,n),px(T,sav,len<<1);
NTT(T.a,len<<1,0),T.cl(len,len<<1);
for(int i=0;i<len;i++)T.a[i]=(R.a[i]-T.a[i]+mod)%mod;
}
T.resiz(sizA),A.resiz(sizA);
return T;
}
Poly sq(Poly A){
int sizA=A.siz(),n=1;
for(;n<sizA;n<<=1);
Poly T,R;
T.resiz(n<<1),R.resiz(n<<1),A.resiz(n<<1);
T.a[0]=Cipolla(A.a[0]);
if(mod-T.a[0]<T.a[0])T.a[0]=mod-T.a[0];
for(int len=2;len<=n;len<<=1){
for(int i=0;i<(len>>1);i++)R.a[i]=(T.a[i]<<1)%mod;
R.cl(len>>1,len);
R=inv(R);
NTT(T.a,len,1),px(T,T,len),NTT(T.a,len,0);
for(int i=0;i<len;i++)T.a[i]=(T.a[i]+A.a[i])%mod;
T.resiz(len<<1),R.resiz(len<<1);
T=mut(T,R);
T.cl(len,len<<1);
}
T.resiz(sizA),A.resiz(sizA);
return T;
}
pair<Poly,Poly>mdiv(Poly A,Poly B){
Poly AT,BT,Q,tmp,R;
R.resiz(B.siz()-1);
int L=A.siz()-B.siz()+1;
AT.resiz(A.siz()),BT.resiz(B.siz());
AT.cpy(A,AT.siz()),BT.cpy(B,BT.siz());
AT.rev(),BT.rev(),AT.resiz(L),BT.resiz(L);
BT=inv(BT),Q=mut(AT,BT);
Q.resiz(L),Q.rev();
tmp=mut(B,Q);
for(int i=0;i<B.siz()-1;i++){
R.a[i]=A.a[i]-tmp.a[i];
if(R.a[i]<0)R.a[i]+=mod;
}
return make_pair(Q,R);
}
Poly dr(Poly A){
Poly res;
res.resiz(A.siz()-1);
for(int i=0;i<res.siz();i++)res.a[i]=A.a[i+1]*(i+1)%mod;
return res;
}
Poly itg(Poly A){
vector<int>pinv;
Poly res;
pinv.resize(A.siz()+1),res.resiz(A.siz()+1);
pinv[1]=1,pinv[0]=0;
for(int i=2;i<pinv.size();i++)pinv[i]=pinv[mod%i]*(mod-mod/i)%mod;
for(int i=1;i<res.siz();i++)res.a[i]=A.a[i-1]*pinv[i]%mod;
res.a[0]=0;
return res;
}
Poly ln(Poly A){
Poly res,B=A;
B=inv(A),A=dr(A);
res=mut(A,B);
res=itg(res);
res.resiz(A.siz()+1);
return res;
}
Poly exp(Poly A){
Poly R,T;
int n=1,siz=A.siz();
for(;n<A.siz();n<<=1);
A.resiz(n<<1),T.resiz(n<<1),R.resiz(n<<1),T.a[0]=1;
for(int len=2;len<=n;len<<=1){
R.cpy(T,len>>1),R=ln(R);
for(int i=0;i<len;i++)R.a[i]=(A.a[i]-R.a[i]+mod)%mod;
R.a[0]=(R.a[0]+1)%mod;
T.resiz(len<<1),R.resiz(len<<1);
T=mut(T,R);
T.cl(len,len<<1);
}
T.resiz(siz);
return T;
}
Poly pqpow(Poly A,int k){
A=ln(A);
for(int i=0;i<A.siz();i++)A.a[i]=A.a[i]*k%mod;
A=exp(A);
return A;
}
int n,m;
Poly f,g;
signed main(){
srand(time(0));
cin.tie(0),cout.tie(0),ios::sync_with_stdio(false);
return 0;
}
跨越大半个 OI 生涯的笔记。
浙公网安备 33010602011771号