0条评论
还没有人评论过~
「学习笔记」多项式
把一直学不懂的各种大常数 \(O(n\log n)\) 的神奇多项式算法总结一下……
证明什么的比较简略……
还有我今天才知道预处理一下单位根会快很多 \(qwq\)
1、\(FFT\)
最没用的一个
首先我们能写出一个 \(O(n^2)\) 的暴力 这个你都不会就可以退役了（某位dalao题解中的）
要确定一个多项式，我们发现只要代 \(f(1),f(2),f(3),...,f(n+m)\) 然后暴力插值就可以了。
还是不行。但是聪明的数学家们发现单位根有分治的特性，然后 \(FFT\) 就发明出来了。
\(FFT\) 分两步，一步 \(DFT\)，一步 \(IDFT\)。
\(Complex:\)
struct Complex{
double x,y;
Complex(double u=0,double v=0){x=u,y=v;}
};
Complex operator + (Complex a,Complex b){
return Complex(a.x+b.x,a.y+b.y);
}
Complex operator - (Complex a,Complex b){
return Complex(a.x-b.x,a.y-b.y);
}
Complex operator * (Complex a,Complex b){
return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
\(FFT:\)
inline void FFT(Complex *f,int op){
for(int i=0;i<n;i++)
if(i<r[i]) swap(f[i],f[r[i]]);
Complex tmp,buf,x,y;
for(int len=1;len<n;len<<=1){
tmp=Complex(cos(Pi/len),op*sin(Pi/len));
for(int R=len<<1,j=0;j<n;j+=R){
buf=Complex(1,0);
for(int k=0;k<len;k++){
x=f[j+k],y=buf*f[j+k+len];
f[j+k]=x+y;f[j+k+len]=x-y;
buf=buf*tmp;
}
}
}
if(op==1) return ;
for(int i=0;i<n;i++) f[i].x=fabs(f[i].x/n);
}
2、\(NTT\)
这个更有用。。。
把复数换成原根什么的就可以了。
预处理：
inline int add(int x,int y){
x+=y;x>=mod?x-=mod:0;
return x;
}
inline int sub(int x,int y){
x-=y;x<0?x+=mod:0;
return x;
}
inline int mul(int x,int y){
return 1ll*x*y%mod;
}
inline void calcrev(int lim){
for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)?(lim>>1):0);
}
inline int fpow(int a,int b){
int ret=1;
for(;b;b>>=1,a=mul(a,a))
if(b&1) ret=mul(ret,a);
return ret;
}
for(int len=1,l=0;len<=mod;len<<=1,l++){
G[l][0]=fpow(3,(mod-1)/len);
G[l][1]=fpow(G[l][0],mod-2);
}
\(NTT:\)
inline void NTT(int *f,int n,int op){
for(int i=0;i<n;i++)
if(i<r[i]) swap(f[i],f[r[i]]);
int buf,tmp,x,y;
for(int len=1,l=1;len<n;len<<=1,l++){
tmp=G[l][0];
if(op==-1) tmp=G[l][1];
for(int i=0;i<n;i+=len<<1){
buf=1;
for(int j=0;j<len;j++){
x=f[i+j];y=mul(buf,f[i+j+len]);
f[i+j]=add(x,y);f[i+j+len]=sub(x,y);
buf=mul(buf,tmp);
}
}
}
if(op==1) return ;
int inv=fpow(n,mod-2);
for(int i=0;i<n;i++) f[i]=mul(f[i],inv);
}
还有比较懒写了几个辅助函数：
inline void Mul(int *A,int *B,int *C,int n,int m){
int lim;
for(lim=1;lim<(n+m);lim<<=1);
calcrev(lim);
NTT(A,lim,1);NTT(B,lim,1);
for(int i=0;i<lim;i++) C[i]=mul(A[i],B[i]);
NTT(C,lim,-1);
}
inline void Clear(int *A,int *B,int *C,int *D,int n){
int lim;
for(lim=1;lim<(n<<1);lim<<=1);
for(int i=0;i<lim;i++) A[i]=B[i]=C[i]=0;
for(int i=n;i<lim;i++) D[i]=0;
}
inline void Clear(int *A,int *B,int *C,int n){
int lim;
for(lim=1;lim<n;lim<<=1);
for(int i=0;i<lim;i++) A[i]=B[i]=C[i]=0;
}
3、多项式求逆
倍增。
\(F(x)*G(x)\equiv 1(mod\ x^{\lceil \frac n2 \rceil})\)
\(F(x)*H(x)\equiv 1(mod\ x^n)\)
\(G(x)-H(x)\equiv 0(mod\ x^{\lceil \frac n2 \rceil})\)
\((G(x)-H(x))^2\equiv 0(mod\ x^n)\)
\(G^2(x)-2*G(x)*H(x)+H^2(x)\equiv 0(mod\ x^n)\)
\(F(x)*G^2(x)-2*G(x)*F(x)*H(x)+F(x)*H^2(x)\equiv 0(mod\ x^n)\)
\(F(x)*G^2(x)-2*G(x)+H(x)\equiv 0(mod\ x^n)\)
\(H(x)=2*G(x)-F(x)*G^2(x)(mod\ x^n)\)
\(Inv:\)
inline void Inv(int *a,int *b,int n){
b[0]=fpow(a[0],mod-2);
static int A[maxn],B[maxn],C[maxn],len,lim;
for(len=1;len<(n<<1);len<<=1){
lim=len<<1;
for(int i=0;i<len;i++) A[i]=a[i],B[i]=b[i];
calcrev(lim);
NTT(A,lim,1);NTT(B,lim,1);
for(int i=0;i<lim;i++) C[i]=mul(sub(2,mul(A[i],B[i])),B[i]);
NTT(C,lim,-1);
for(int i=0;i<len;i++) b[i]=C[i];
}
Clear(A,B,C,b,n);
}
4、多项式除法+取模
\(F(x)=Q(x)*G(x)+R(x)\)
\(F(\frac 1x)=Q(\frac 1x)*G(\frac 1x)+R(\frac 1x)\)
\(x^nF(\frac 1x)=x^{n-m}Q(\frac 1x)*x^mG(\frac 1x)+x^{n-m+1}*x^{m-1}*R(\frac 1x)\)
\(F_R(x)=Q_R(x)*G_R(x)+x^{n-m+1}*R_R(x)\)
\(F_R(x)=Q_R(x)*G_R(x)(mod\ x^{n-m+1})\)
然后到这里 \(Q_R(x)\) 就可以求出来了，求出来由 \(R(x)=F(x)-Q(x)*G(x)\) 求出 \(R(x)\)
\(Div:\)
inline void Div(int *a,int *b,int *c,int n,int m){
static int A[maxn],B[maxn],C[maxn];
reverse(a,a+n);reverse(b,b+m);
for(int i=0;i<n-m+1;i++) A[i]=a[i];
Inv(b,B,n-m+1);
Mul(A,B,C,n-m+1,n-m+1);
for(int i=0;i<n-m+1;i++) c[i]=C[i];
reverse(c,c+n-m+1);
reverse(a,a+n);reverse(b,b+m);
Clear(A,B,C,(n-m+1)<<1);
}
\(Mod:\)
inline void Mod(int *a,int *b,int *c,int n,int m){
static int A[maxn],B[maxn],C[maxn];
for(int i=0;i<m;i++) A[i]=b[i];
Div(a,b,B,n,m);
Mul(A,B,C,n-m+1,m);
for(int i=0;i<m-1;i++) c[i]=sub(a[i],C[i]);
Clear(A,B,C,n+1);
}
5、多项式开根
\(H^2(x)\equiv F(x)(mod\ x^{\lceil\frac n2 \rceil})\)
\(G(x)\equiv H(x)(mod\ x^{\lceil\frac n2 \rceil})\)
\(G(x)-H(x)\equiv 0(mod\ x^{\lceil\frac n2 \rceil})\)
\((G(x)-H(x))^2\equiv 0(mod\ x^{\lceil\frac n2 \rceil})\)
\(G^2(x)-2H(x)*G(x)+H^2(x)\equiv 0(mod\ x^n)\)
\(F(x)-2H(x)*G(x)+H^2(x)\equiv 0(mod\ x^n)\)
\(G(x)=\Large \frac {F(x)+H^2(x)}{2H(x)}\)
\(Sqrt:\)
inline void Sqrt(int *a,int *b,int n){
b[0]=1;
static int A[maxn],B[maxn],C[maxn],len;
for(len=1;len<(n<<1);len<<=1){
for(int i=0;i<len;i++) A[i]=a[i];
for(int i=0;i<len;i++) B[i]=0;//用static的话这里一定要清零
Inv(b,B,len);
Mul(A,B,C,len,len);
for(int i=0;i<len;i++) b[i]=mul(add(b[i],C[i]),inv2);
}
Clear(A,B,C,b,n);
}
6、分治 \(FFT\)
用求逆的话不是很傻吗。。。
我们考虑 \([l,mid]\) 对 \([mid+1,r]\) 的贡献，构造多项式卷积。
\(CDQFFT:\)
void CDQFFT(int *f,int *g,int l,int r){
if(l==r) return ;
int mid=(l+r)>>1,lim;
CDQFFT(f,g,l,mid);
for(lim=1;lim<=((r-l+1)<<1);lim<<=1);
for(int i=l;i<=mid;i++) A[i-l]=f[i];
for(int i=0;i<=r-l;i++) B[i]=g[i];
calcrev(lim);
NTT(A,lim,1);NTT(B,lim,1);
for(int i=0;i<lim;i++) A[i]=1ll*A[i]*B[i]%mod;
NTT(A,lim,-1);
for(int i=mid+1;i<=r;i++) f[i]=(f[i]+A[i-l])%mod;
for(int i=0;i<lim;i++) A[i]=B[i]=0;
CDQFFT(f,g,mid+1,r);
}
7、多项式多点求值
构造多项式 \(\prod_{i=l}^{r}(x-a_i)\)，类似分治 \(FFT\) 的思路，但是要用到多项式取模。
我们用类似线段树的方式存下所有多项式。
\(Build+Eval:\)
void Build(int *a,int l,int r,int x){
if(l==r){
P[x].push_back(mod-a[l]);
P[x].push_back(1);
return ;
}
int mid=(l+r)>>1;
Build(a,l,mid,x<<1);
Build(a,mid+1,r,x<<1|1);
int n=P[x<<1].size(),m=P[x<<1|1].size();
static int A[maxn],B[maxn],C[maxn];
for(int i=0;i<n;i++) A[i]=P[x<<1][i];
for(int i=0;i<m;i++) B[i]=P[x<<1|1][i];
Mul(A,B,C,n,m);
for(int i=0;i<n+m-1;i++) P[x].push_back(C[i]);
Clear(A,B,C,n+m);
}
void Eval(int *a,int *b,int dep,int n,int l,int r,int x){
int mid=(l+r)>>1,m=P[x].size();
static int A[maxn];int *B=Q[dep];
for(int i=0;i<m;i++) A[i]=P[x][i];
Mod(a,A,B,n,m);
if(l==r){b[l]=B[0];return;}
Eval(B,b,dep+1,m-1,l,mid,x<<1);
Eval(B,b,dep+1,m-1,mid+1,r,x<<1|1);
Clear(A,B,B,m-1);
}
8、多项式快速插值
咕咕咕。
\(calc+solve:\)
void calc(int *a,int *b,int dep,int l,int r,int x){
if(l==r){b[0]=a[l];return;}
int mid=(l+r)>>1,n=P[x<<1].size(),m=P[x<<1|1].size();
static int A[maxn],C[maxn];int *B=Q[dep];
calc(a,B,dep+1,l,mid,x<<1);
for(int i=0;i<m;i++) A[i]=P[x<<1|1][i];
Mul(A,B,C,m,n);
for(int i=0,j=P[x].size();i<j;i++) b[i]=add(b[i],C[i]);
Clear(A,B,C,n+m);
calc(a,B,dep+1,mid+1,r,x<<1|1);
for(int i=0;i<n;i++) A[i]=P[x<<1][i];
Mul(A,B,C,n,m);
for(int i=0,j=P[x].size();i<j;i++) b[i]=add(b[i],C[i]);
Clear(A,B,C,n+m);
}
inline void solve(int *a,int *b,int *c,int n){
Build(a,0,n-1,1);
int m=P[1].size();
static int A[maxn],B[maxn];
for(int i=0;i<m;i++) A[i]=P[1][i];
Direv(A,A,m);Eval(A,B,0,m-1,0,n-1,1);
for(int i=0;i<n;i++) B[i]=mul(b[i],fpow(B[i],mod-2));
calc(B,c,0,0,n-1,1);
for(int i=0;i<n;i++) A[i]=B[i]=0;
}
前五个是 \(O(n\log n)\) 的，后面三个是 \(O(n\log^2 n)\) 的。
\(Upd:\) 居然忘了多项式 \(Ln\) 和 \(Exp\)。。。
9、\(MTT\)
拆系数 \(FFT\)，实际使用比三模数要快。不过我的 \(MTT\) 必须要开 \(long\ double\) 是什么鬼啊。。。
\(MTT:\)
inline void Mul(int *a,int *b,int *c,int n,int m){
int lim;
for(lim=1;lim<(n+m);lim<<=1);
for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)?(lim>>1):0);
for(int i=0;i<lim;i++){
a[i]%=mod;b[i]%=mod;
A[i]=Complex(a[i]/base,0);
B[i]=Complex(a[i]%base,0);
C[i]=Complex(b[i]/base,0);
D[i]=Complex(b[i]%base,0);
}
FFT(A,lim,1);FFT(B,lim,1);FFT(C,lim,1);FFT(D,lim,1);
Complex tmp;
for(int i=0;i<lim;i++){
tmp=A[i]*C[i];C[i]=B[i]*C[i];
B[i]=B[i]*D[i];D[i]=D[i]*A[i];
A[i]=tmp;C[i]=C[i]+D[i];
}
FFT(A,lim,-1);FFT(B,lim,-1);FFT(C,lim,-1);
for(int i=0;i<lim;i++){
c[i]=0;
c[i]=(c[i]+(ll)A[i].x%mod*base%mod*base%mod)%mod;
c[i]=(c[i]+(ll)C[i].x%mod*base%mod)%mod;
c[i]=(c[i]+(ll)B[i].x%mod)%mod;
c[i]=(c[i]+mod)%mod;
}
}
来源：https://www.cnblogs.com/owencodeisking/p/10353700.html
标签：