「学习笔记」多项式

「学习笔记」多项式

把一直学不懂的各种大常数 (O(nlog 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 ceil}))

(F(x)*H(x)equiv 1(mod x^n))

(G(x)-H(x)equiv 0(mod x^{lceil frac n2 ceil}))

((G(x)-H(x))^2equiv 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^{lceilfrac n2 ceil}))

(G(x)equiv H(x)(mod x^{lceilfrac n2 ceil}))

(G(x)-H(x)equiv 0(mod x^{lceilfrac n2 ceil}))

((G(x)-H(x))^2equiv 0(mod x^{lceilfrac n2 ceil}))

(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(nlog n)) 的，后面三个是 (O(nlog^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;
    }
}