多项式细节梳理&模板（多项式）
基础
现在的码风又变了。。。
FFT和NTT的板子
typedef complex<double> C;
const double PI=acos(-1);
void FFT(C*a,R op){
for(R i=0;i<N;++i)
if(i<r[i])swap(a[i],a[r[i]]);
for(R i=1;i<N;i<<=1){
C wn=C(cos(PI/i),sin(PI/i)*op),w=1,t;
for(R j=0;j<N;j+=i<<1,w=1)
for(R k=j;k<j+i;++k,w*=wn)
t=a[k+i]*w,a[k+i]=a[k]-t,a[k]+=t;
}
}
const int YL=998244353;
LL Pow(LL b,R k=YL-2,LL a=1){
for(;k;k>>=1,b=b*b%YL)
if(k&1)a=a*b%YL;
return a;
}
void NTT(R*a,R n,R op){
for(R i=0;i<n;++i)
if(i<r[i])swap(a[i],a[r[i]]);
for(R i=1;i<n;i<<=1){
LL wn=Pow(3,(YL-1)/(i<<1)),w=1,x;
if(op==-1)wn=Pow(wn);
for(R j=0;j<n;j+=i<<1,w=1)
for(R k=j;k<j+i;++k,w=w*wn%YL)
x=w*a[k+i]%YL,a[k+i]=Mod(a[k]-x+YL),a[k]=Mod(a[k]+x);
}
if(op==-1){
LL x=Pow(n);
for(R i=0;i<n;++i)a[i]=a[i]*x%YL;
}
}
任意模数NTT
怎么看都觉得MTT给人感觉不好，用double巨佬说怕掉精度，用longdouble那常数想象一下～
然后就去学了三模NTT
参考Memory of Winter巨佬的题解
细节不多吧，首先要把\(998244353,1004535809,469762049\)背下来，它们都有一个原根为\(3\)。
然后手动两两合并同余方程，不要直接用CRT公式。
蒟蒻为了减少码量直接开二维数组了，实际上像Memory of Winter巨佬一样封装结构体会跑得更快。
upd:重写了一遍，二维数组又长又丑又慢qwq
洛谷P4245 【模板】任意模数NTT
#include<bits/stdc++.h>
#define LL long long
#define I inline
#define R register int
#define G if(++ip==ie)if(fread(ip=buf,1,SZ,stdin))
#define Wn(A) Pow(3,Mod((A-1)/(i<<1)*op,A-1),A)
using namespace std;
const int SZ=1<<19,N=1<<18,YL=1e9+7,A=998244353,B=1004535809,C=469762049;
char buf[SZ],*ie=buf+SZ,*ip=ie-1;
inline int in(){
G;while(*ip<'-')G;
R x=*ip&15;G;
while(*ip>'-'){x*=10;x+=*ip&15;G;}
return x;
}
int L,r[N];
I int Mod(R x,R YL){
return x+(x>>31&YL);
}
I int Pow(LL b,R k,LL YL,LL a=1){
for(;k;k>>=1,b=b*b%YL)
if(k&1)a=a*b%YL;
return a;
}
I int Inv(LL b,LL YL){
return Pow(b%YL,YL-2,YL);
}
struct Z{
int a,b,c;
Z(){}
Z(LL a):a(a%A),b(a%B),c(a%C){}
Z(R a,R b,R c):a(a),b(b),c(c){}
I Z operator!(){a+=a>>31&A,b+=b>>31&B,c+=c>>31&C;return*this;}
I Z operator+(const Z&x){return!Z(a+x.a-A,b+x.b-B,c+x.c-C);}
I Z operator-(const Z&x){return!Z(a-x.a ,b-x.b ,c-x.c );}
I Z operator*(const Z&x){return Z((LL)a*x.a%A,(LL)b*x.b%B,(LL)c*x.c%C);}
I int CRT(LL YL){
static LL AB=(LL)A*B%YL,I0=Inv(A,B),I1=Inv((LL)A*B,C),x;
x=Mod(b-a,B)*I0%B*A+a;
return(Mod(c-x%C,C)*I1%C*AB+x)%YL;
}
}f[N],g[N];
void NTT(Z*a,R op){
for(R i=0;i<L;++i)
if(i<r[i])swap(a[i],a[r[i]]);
for(R i=1;i<L;i<<=1){
Z wn(Wn(A),Wn(B),Wn(C)),w(1),x;
for(R j=0;j<L;j+=i<<1,w=1)
for(R k=j;k<j+i;++k,w=w*wn)
x=w*a[k+i],a[k+i]=a[k]-x,a[k]=a[k]+x;
}
if(op==-1){
Z w(Inv(L,A),Inv(L,B),Inv(L,C));
for(R i=0;i<L;++i)a[i]=a[i]*w;
}
}
int main(){
R n=in(),m=in(),p=in();
for(R i=0;i<=n;++i)f[i]=in();
for(R i=0;i<=m;++i)g[i]=in();
for(L=1;L<=n+m;L<<=1);
for(R i=1;i<L;++i)r[i]=(r[i>>1]|L*(1&i))>>1;
NTT(f,1);NTT(g,1);
for(R i=0;i<L;++i)f[i]=f[i]*g[i];
NTT(f,-1);
for(R i=0;i<=n+m;++i)printf("%d ",f[i].CRT(p));
return 0;
}
另外，三模NTT不能将三个及以上的多项式一次都乘起来，因为实际值域太大，不能保证模意义下的唯一性。
正确的做法是两个两个乘。
比如，任意模数多项式求逆的核心代码
洛谷P4239 【模板】多项式求逆（加强版）
void PolyInv(Z*a,Z*b,Z*a1,R m){
if(m==1){b[0]=Inv(a[0].b,YL);return;}
PolyInv(a,b,a1,(m+1)>>1);memcpy(a1,a,12*m);
Init(2*m);NTT(b,1);NTT(a1,1);
for(R i=0;i<L;++i)a1[i]=b[i]*a1[i];
NTT(a1,-1);
for(R i=0;i<L;++i)a1[i]=Mod(-a1[i].CRT()+2*(i==0),YL);
NTT(a1,1);
for(R i=0;i<L;++i)b[i]=b[i]*a1[i];
NTT(b,-1);memset(a1,0,12*L);memset(b+m,0,12*(L-m));
for(R i=0;i<m;++i)b[i]=b[i].CRT();
}
多项式运算
分治乘法
若干个总项数不超过\(n\)的多项式的卷积，根据合并果子的方法，可以在不超过\(O(n\log^2n)\)的时间内完成。
有时候是若干个二项式相乘，还可以写成循环的形式，好处是减少Rader排序r数组的预处理次数。
第一次尝试这种写法：洛谷CF981H K Paths
void Solve(R x){
R n=4*e[x].size();
for(R i=0;i<n;i+=4)
a[i]=1,a[i+1]=s[e[x][i>>2]];
for(R i=4;i<n;i<<=1){
for(R j=1;j<i;++j)
r[j]=(r[j>>1]>>1)|(1&j?i>>1:0);
for(R j=0;j+i<n;j+=i<<1){
R*p=a+j,*q=p+i;
NTT(p,i,1);NTT(q,i,1);
for(R k=0;k<i;++k)p[k]=(LL)p[k]*q[k]%YL;
NTT(p,i,-1);
memset(q,0,4*i);
}
}
}
分治FFT
某些卷积式有一些特点，直接算计算量庞大，却可以通过CDQ分治思想，考虑左边对右边的贡献。
洛谷P4721 【模板】分治 FFT
void Solve(R l,R r){
if(l+1==r)return;
R m=(l+r)>>1;
Solve(l,m);R n=Init(m-l+r-l);
memcpy(a,f+l,4*(m-l));memset(a+m-l,0,4*(n-m+l));
memcpy(b,g ,4*(r-l));memset(b+r-l,0,4*(n-r+l));
NTT(a,n,1);NTT(b,n,1);
for(R i=0;i<n;++i)a[i]=(LL)a[i]*b[i]%YL;
NTT(a,n,-1);
for(R i=m;i<r;++i)f[i]=Mod(f[i]+a[i-l]);
Solve(m,r);
}
泰勒展开
将函数\(f(x)\)在\(x_0\)处展开得
\(f(x)=f(x_0)+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+...+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+\xi\)
牛顿迭代推导倍增求解式
已知多项式函数\(F(x)\)，用倍增法求解\(B\)使得\(F(B)\equiv0(\mod x^{2^k})\)。
假设已经求出\(B,F(B)\equiv0(\mod x^{2^k})\)，现在求\(B_1,F(B_1)\equiv0(\mod x^{2^{k+1}})\)。
\(F(B_1)\)在\(x=B\)处只展开一项得到
\(F(B_1)=F(B)+F'(B)(B_1-B)\equiv0(\mod x^{2^{k+1}})\)
\(B_1=B-\frac{F(B)}{F'(B)}\)，带入即可。
接下来的\(A\)为已知多项式，B为待求多项式。
求逆
\(F(B)=AB-1=0\)
\(B_1=B-\frac{AB-1}{A}=B-B(AB-1)=2B-AB^2\)
开方
\(F(B)=B^2-A=0\)
\(B_1=B-\frac{B^2-A}{2B}=\frac{B^2+A}{2B}\)
对数
没必要迭代。
\(B=\ln A\)
\(B=\int\frac{A'}{A}\)
指数
\(B=e^A\)
\(\ln B-A=0\)
\(B_1=B-\frac{\ln B-A}{\frac1B}=B(1+A-\ln B)\)
快速幂
\(B=A^k=e^{k\ln A}\)
这常数看着就吓人。。。
除法&取余
https://www.luogu.org/problemnew/solution/P4512
代码模板
写的时候保证了一定的稳定性（数组清空问题），因此牺牲了一丁点效率。
参数中，a为已知，b为待求，a1、b1为辅助数组。
传入前需自己确保b、a1、b1为空。
返回后保证a维持原状，b存好答案，a1、b1为空。
还有，线性递推，多点求值，插值待填坑
const int YL=998244353;
int Inv[N],r[N];//在外面初始化逆元
inline int Mod(const int x){
return x>=YL?x-YL:x;
}
LL Pow(LL b,R k=YL-2,LL a=1){
for(;k;k>>=1,b=b*b%YL)
if(k&1)a=a*b%YL;
return a;
}
int Init(R m){
R n=1;while(n<m<<1)n<<=1;
for(R i=0;i<n;++i)r[i]=(r[i>>1]|(1&i)*n)>>1;
return n;
}
void NTT(R*a,R n,R op){
for(R i=0;i<n;++i)
if(i<r[i])swap(a[i],a[r[i]]);
for(R i=1;i<n;i<<=1){
LL wn=Pow(3,(YL-1)/(i<<1)),w=1,x;
if(op==-1)wn=Pow(wn);
for(R j=0;j<n;j+=i<<1,w=1)
for(R k=j;k<j+i;++k,w=w*wn%YL)
x=w*a[k+i]%YL,a[k+i]=Mod(a[k]-x+YL),a[k]=Mod(a[k]+x);
}
if(op==-1){
LL x=Pow(n);
for(R i=0;i<n;++i)a[i]=a[i]*x%YL;
}
}
void PolyRev(R*a,R*b,R n){
for(R i=0;i<n;++i)b[i]=a[n-i-1];
}
void PolyDer(R*a,R*b,R n){
for(R i=1;i<n;++i)b[i-1]=(LL)a[i]*i%YL;
if(a==b)a[n-1]=0;
}
void PolyInt(R*a,R*b,R n){
for(R i=n;i;--i)b[i]=(LL)a[i-1]*Inv[i]%YL;
if(a==b)a[0]=0;
}
void PolyInv(R*a,R*b,R*a1,R m){
if(m==1){b[0]=Pow(a[0]);return;}
PolyInv(a,b,a1,(m+1)>>1);memcpy(a1,a,4*m);
R n=Init(m);NTT(b,n,1);NTT(a1,n,1);
for(R i=0;i<n;++i)b[i]=(YL+2-(LL)b[i]*a1[i]%YL)*b[i]%YL;
NTT(b,n,-1);memset(a1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolySqrt(R*a,R*b,R*a1,R*b1,R m){
if(m==1){b[0]=sqrt(a[0]);return;}
PolySqrt(a,b,a1,b1,(m+1)>>1);PolyInv(b,b1,a1,m);memcpy(a1,a,4*m);
R n=Init(m);NTT(a1,n,1);NTT(b1,n,1);
for(R i=0;i<n;++i)a1[i]=(LL)a1[i]*b1[i]%YL;
NTT(a1,n,-1);
for(R i=0;i<m;++i)b[i]=(LL)(b[i]+a1[i])*((YL+1)>>1)%YL;
memset(a1,0,4*n);memset(b1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolyLn(R*a,R*b,R*a1,R m){
PolyInv(a,b,a1,m);PolyDer(a,a1,m);
R n=Init(m);NTT(b,n,1);NTT(a1,n,1);
for(R i=0;i<n;++i)b[i]=(LL)b[i]*a1[i]%YL;
NTT(b,n,-1);memset(a1,0,4*n);PolyInt(b,b,m);
}
void PolyExp(R*a,R*b,R*a1,R*b1,R m){
if(m==1){b[0]=1;return;}
PolyExp(a,b,a1,b1,(m+1)>>1);PolyLn(b,b1,a1,m);memcpy(a1,a,4*m);
R n=Init(m);NTT(b,n,1);NTT(a1,n,1);NTT(b1,n,1);
for(R i=0;i<n;++i)b[i]=(LL)b[i]*(YL+1+a1[i]-b1[i])%YL;
NTT(b,n,-1);memset(a1,0,4*n);memset(b1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolyDiv(R*A,R*a,R*b,R*A1,R*a1,R M,R m){
PolyRev(a,a1,m);PolyInv(a1,b,A1,M-m+1);PolyRev(A,A1,M);
R n=Init(M);NTT(b,n,1);NTT(A1,n,1);
for(R i=0;i<n;++i)b[i]=(LL)b[i]*A1[i]%YL;
NTT(b,n,-1);memset(A1,0,4*n);memset(a1,0,4*n);
reverse(b,b+M-m+1);memset(b+M-m+1,0,4*(n-M+m-1));
}
void PolyMod(R*A,R*a,R*b,R*A1,R*a1,R M,R m){
PolyDiv(A,a,b,A1,a1,M,m);memcpy(A1,A,4*M);memcpy(a1,a,4*m);
R n=Init(M);NTT(b,n,1);NTT(A1,n,1);NTT(a1,n,1);
for(R i=0;i<n;++i)b[i]=Mod(A1[i]-(LL)b[i]*a1[i]%YL+YL);
NTT(b,n,-1);memset(A1,0,4*n);memset(a1,0,4*n);
}
