[学习笔记]多项式

把一直学不懂的各种大常数 \(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)\) 的。

[学习笔记]多项式

原文：https://www.cnblogs.com/owencodeisking/p/10353700.html