学习笔记:多项式变换

多项式学习笔记

卷积形式:

一般卷积形式

Fi=jk=iajbk

+是为常见的多项式加法卷积,×有时也能转为加法,位运算是FWT,整除是狄利克雷卷积

加法卷积常用形式:

Fi=j=0iaj×aij

卷积满足交换律,结合律,分配律

FFT

讲解

利用单位根和复数的性质分治的去把多项式系数表达式转成点值表达式, O(n) 把对应点值乘起来后再转成系数表达式

复杂度O(nlogn)

递归形式
/*
FFT递归实现
*/

#include <bits/stdc++.h>
#define endl putchar('\n')
using namespace std;
const int M = 7e6+10;
const int inf = 2147483647;
const int Mod = 1e9+7;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
struct I{
    double x,y;
    I(){x=y=0;}
    I(double a,double b){x=a;y=b;}
    friend I operator +(I a,I b){return I(a.x+b.x,a.y+b.y);}
    friend I operator -(I a,I b){return I(a.x-b.x,a.y-b.y);}
    friend I operator *(I a,I b){return I(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}a[M],b[M];
void FFT(int n,I *a,int type){
    if(n==1) return;
    I a1[n>>1],a2[n>>1];
    for(int i=0;i<n;i+=2){
        a1[i>>1]=a[i];a2[i>>1]=a[i+1];
    }
    FFT(n>>1,a1,type);
    FFT(n>>1,a2,type);
    I wn=I(cos(2.0*pi/n),type*sin(2.0*pi/n)),w=I(1,0);
    for(int i=0;i<(n>>1);i++,w=w*wn){
        a[i]=a1[i]+w*a2[i];
        a[i+(n>>1)]=a1[i]-w*a2[i];
    }
}
signed main(){
    int n=read(),m=read();
    for(int i=0;i<=n;i++) a[i].x=read();
    for(int i=0;i<=m;i++) b[i].x=read();
    int lim=1;
    while(lim<=n+m) lim<<=1;
    FFT(lim,a,1);
    FFT(lim,b,1);
    for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
    FFT(lim,a,-1);
    for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/lim+0.5));endl;
    return 0;
}

迭代形式
/*
FFT迭代实现
*/

#include <bits/stdc++.h>
#define endl putchar('\n')
using namespace std;
const int M = 5e6+10;
const int inf = 2147483647;
const int Mod = 1e9+7;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
struct I{
    double x,y;
    I(){x=y=0;}
    I(double a,double b){x=a;y=b;}
    friend I operator +(I a,I b){return I(a.x+b.x,a.y+b.y);}
    friend I operator -(I a,I b){return I(a.x-b.x,a.y-b.y);}
    friend I operator *(I a,I b){return I(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}a[M],b[M];
int lim=1,r[M];
void FFT(I *a,int type){
    for(int i=0;i<lim;i++){
        if(r[i]>i) swap(a[i],a[r[i]]);
    }
    for(int mid=1;mid<lim;mid<<=1){
        I wn(cos(pi/mid),type*sin(pi/mid));
        for(int len=mid<<1,j=0;j<lim;j+=len){
            I w=I(1,0);
            for(int k=0;k<mid;k++,w=w*wn){
                I x=a[j+k],y=w*a[j+k+mid];
                a[j+k]=x+y;
                a[j+k+mid]=x-y;
            }
        }
    }
}
int main(){
    int n=read(),m=read();
    for(int i=0;i<=n;i++) a[i].x=read();
    for(int i=0;i<=m;i++) b[i].x=read();
    int l=0;
    while(lim<=n+m) lim<<=1,l++;
    for(int i=0;i<lim;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    }
    FFT(a,1);FFT(b,1);
    for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/lim+0.5));
    return 0;
}

例题

[ZJOI2014]力


给出 n 个数 q1,q2,qn,定义

Fj = i=1j1qi×qj(ij)2  i=j+1nqi×qj(ij)2

Ei = Fiqi

1in,求 Ei 的值。


把式子换成卷积形式,常用技巧有翻转序列,定义初值以扩大枚举范围等

Fj=i=1j1qi×qj(ij)2i=j+1nqi×qj(ij)2=i=1jqi(ij)2i=jnqi(ij)2

gi=1i2,g0=0,q0=0qi=qni,t=nj

Fj=i=0jqi×gjii=0njqi+j×gi=i=0jqi×gjii=0tqti×gi

然后FFT

代码
#include <bits/stdc++.h>
#define endl putchar('\n')
using namespace std;
const int M = 5e5+10;
const int inf = 2147483647;
const int Mod = 1e9+7;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
struct I{
    double x,y;
    I(){x=y=0;}
    I(double a,double b){x=a;y=b;}
    friend I operator +(I a,I b){return I(a.x+b.x,a.y+b.y);}
    friend I operator -(I a,I b){return I(a.x-b.x,a.y-b.y);}
    friend I operator *(I a,I b){return I(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}q[M],qp[M],g[M],tmp[M];
int lim,l,r[M];
void FFT(I *a,int type){
    for(int i=0;i<lim;i++) if(i>r[i]) swap(a[i],a[r[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        I wn(cos(pi/mid),type*sin(pi/mid));
        for(int len=mid<<1,j=0;j<lim;j+=len){
            I w(1,0);
            for(int k=0;k<mid;k++,w=w*wn){
                I x=a[k+j],y=w*a[k+j+mid];
                a[k+j]=x+y;a[k+j+mid]=x-y;
            }
        }
    }
}
void RollChicken(I *a,I *b,int n,int m){
    lim=1;l=0;
    while(lim<=n+m) lim<<=1,l++;
    for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1);FFT(b,1);
    for(int i=0;i<lim;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=n+m;i++) a[i].x=a[i].x/lim; 
}
int main(){
    int n=read();
    for(int i=1;i<=n;i++){
        tmp[i].x=g[i].x=(double)(1.0/i/i);
        scanf("%lf",&q[i].x);
    }
    for(int i=0;i<=n;i++) qp[i]=q[n-i];
    RollChicken(q,g,n,n);
    RollChicken(qp,tmp,n,n);
    for(int i=1;i<=n;i++){
        printf("%.3lf\n",q[i].x-qp[n-i].x);
    }
    return 0;
}

NTT

原根:

阶:
m>1,gcd(a,m)=1,使得ar1(modm)的最小正整数 r 称为a对模m的阶,记作δm(a)

性质:δa(m)|ϕ(m) (欧拉定理证)

原根:若有 m>1δp(m)=ϕ(p),则称 mp 的原根

性质:

1.一个正整数 m 有原根的充要条件是m=2,4,pe,2pe

2.有原根的正整数 mϕ(ϕ(m)) 个原根

3.若 gm 的一个原根,则

g,g2gϕ(m)

互不相同,且值域为1到ϕ(m)

用处:

1.原根同样满足单位根的一些性质,可以代替单位根以实现NTT,用于保证精度和取模运算

2.用性质3可以把数 x 转成 loggx(modm),把乘法运算变成加法运算,注意如果卷积是取模意义下的根据欧拉定理要模ϕ(m)

即把

Fx=ijmodp=xaibj

变成

FX=I+Jmodϕ(p)aIbJ

其中X=loggxmodp,I=loggimodp,J=loggjmodp

求法:

ϕ(m)分解质因数,枚举g,若恒满足

gϕ(m)pi1(modm)

gm 的原根

NTT实现

把单位根换成 g(p1)/n即可,记得取模

封装带modint
#include <bits/stdc++.h>
#define endl putchar('\n')
using namespace std;
const int M = 3e6+10;
const int inf = 2147483647;
const int Mod = 1004535809,G = 3,Gi = 334845270;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
struct sg{
    int x;
    sg(){x=0;}
    sg(int a){x=a;}
    sg& operator =(int a){return x=a,*this;}
    friend sg operator +(sg a,sg b){return sg(a.x+b.x>=Mod?a.x+b.x-Mod:a.x+b.x);}
    friend sg operator -(sg a,sg b){return sg(a.x-b.x<0?a.x-b.x+Mod:a.x-b.x);}
    friend sg operator *(sg a,sg b){return sg(1ll*a.x*b.x>=Mod?1ll*a.x*b.x%Mod:a.x*b.x);}
    sg& operator +=(sg a){return *this=*this+a,*this;}
    sg& operator *=(sg a){return *this=*this*a,*this;}
};
sg qpow(int xx,int k){
    sg res=1,x=xx;
    while(k){
        if(k&1) res=res*x;
        x=x*x;
        k>>=1;
    }
    return res;
}
struct Polygen{
    sg a[M],b[M];
    int r[M],lim,l;
    void NTT(sg *a,int type){
        for(int i=0;i<lim;i++) if(r[i]>i) swap(a[i],a[r[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            sg wn=qpow(type==1?G:Gi,(Mod-1)/(mid<<1));
            for(int len=mid<<1,j=0;j<lim;j+=len){
                sg w=1;
                for(int k=0;k<mid;k++,w=w*wn){
                    sg x=a[j+k],y=w*a[j+k+mid];
                    a[j+k]=x+y;a[j+k+mid]=x-y;
                }
            }
        }
    }
    void Mul(int *A,int *B,int *C,int n,int m){
        l=0;lim=1;
        while(lim<=n+m) lim<<=1,l++;
        for(int i=0;i<lim;i++) a[i]=b[i]=0,r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        for(int i=0;i<=n;i++) a[i]=A[i];
        for(int i=0;i<=m;i++) b[i]=B[i];
        NTT(a,1);NTT(b,1);
        for(int i=0;i<lim;i++) a[i]=a[i]*b[i];
        NTT(a,-1);
        sg inv=qpow(lim,Mod-2);
        for(int i=0;i<lim;i++) C[i]=(a[i]*inv).x;
    }
}Poly;
int n,m,a[M],b[M],ans[M];
signed main(){
    n=read();m=read();
    for(int i=0;i<=n;i++) a[i]=read();
    for(int i=0;i<=m;i++) b[i]=read();
    Poly.Mul(a,b,ans,n,m);
    for(int i=0;i<=n+m;i++) printf("%d ",ans[i]);
    return 0;
}

例题

[SDOI2015]序列统计

ht(x) 表示从集合中选 t 个数加和为 i 的方案数,有

ht(x)=i=1nht1(i)×h1(xi)

这个过程可以用倍增加NTT的方法快速求解

然后考虑把乘法变成加法,应用上面提到的的技巧,把 x 变成 loggx即可

代码
#include <bits/stdc++.h>
#define int long long
#define endl putchar('\n')
using namespace std;
const int M = 1e6+10;
const int inf = 2147483647;
const int Mod = 1004535809,G = 3,Gi = 334845270;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
ll qpow(ll x,ll k,ll mod){
    ll res=1;
    while(k){
        if(k&1) res=res*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return res;
}
struct Polygen{
    int lim,l,rev[M];
    void NTT(ll *a,int type){
        for(int i=0;i<lim;i++) if(i>rev[i]) swap(a[i],a[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            ll wn=qpow(type==1?G:Gi,(Mod-1)/(mid<<1),Mod);
            for(int len=mid<<1,j=0;j<lim;j+=len){
                ll w=1;
                for(int k=0;k<mid;k++,w=w*wn%Mod){
                    ll x=a[j+k],y=w*a[j+k+mid]%Mod;
                    a[j+k]=(x+y)%Mod;
                    a[j+k+mid]=(x-y+Mod)%Mod;
                }
            }
        }
    }
    ll A[M],B[M];
    void Mul(ll *aa,ll *bb,ll *c,ll n,ll m,ll modx){
        lim=1;l=0;
        while(lim<=n+m) lim<<=1,l++;
        for(int i=0;i<lim;i++) A[i]=B[i]=0;
        for(int i=0;i<=n;i++) A[i]=aa[i];
        for(int i=0;i<=m;i++) B[i]=bb[i];
        for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
        NTT(A,1);NTT(B,1);
        for(int i=0;i<lim;i++) A[i]=A[i]*B[i]%Mod;
        NTT(A,-1);
        ll inv=qpow(lim,Mod-2,Mod);
        for(int i=0;i<=n+m;i++) A[i]=A[i]*inv%Mod;
        for(int i=0;i<modx;i++) c[i]=A[i];
        for(int i=modx;i<lim;i++) c[i%modx]=(c[i%modx]+A[i])%Mod;
        for(int i=modx;i<lim;i++) c[i]=0;
    }
}Ploy;
ll n,m,x,s,tot,f[34][M],fac[M],ans[M],trans[M];
void calc(int k){
    ans[0]=1;
    for(int i=32;i>=0;i--){
        if((k>>i)&1) Ploy.Mul(f[i],ans,ans,m-2,m-2,m-1);
    }
}
int getg(int x){
    int tmp=x-1;
    for(int i=2;i*i<=x-1;i++){
        if(tmp%i==0){
            fac[++tot]=i;
            while(tmp%i==0) tmp/=i;
        }
    }
    if(tmp>1) fac[++tot]=tmp;
    for(int i=2;i<x;i++){
        bool flag=1;
        for(int j=1;j<=tot;j++){
            if(qpow(i,(x-1)/fac[j],x)==1){
                flag=0;break;
            }
        }
        if(flag) return i;
    } 
    return -1;
}
signed main(){
    n=read();m=read();x=read();s=read();
    int g=getg(m);
    for(int i=0,tmp=1;i<m-1;i++,tmp=tmp*g%m){
        trans[tmp]=i;
    }
    for(int i=1;i<=s;i++){
        int x=read();
        if(!x) continue;
        f[0][trans[x]]++;
    }
    for(int i=1;i<=32;i++) Ploy.Mul(f[i-1],f[i-1],f[i],m-2,m-2,m-1);
    calc(n);
    printf("%lld\n",ans[trans[x]]);
    return 0;
}

分治FFT

用于求解:


给定序列 g1n1,求序列 f0n1

其中 fi=j=1ifijgj,边界为 f0=1


考虑分治,先求出 [l,mid]fx 值,把左边对右边的贡献加上,再算右边的值

考虑怎么计算中间贡献,设 flfmidfx 的贡献为 wx

wx=i=lmidfi×gxi

fmid+1r=0,有

wx=i=lx1fi×gxiwx=i=0xl1fi+l×gxil

ai=fi+l,bi=gi+1,有

wx=i=0xl1ai×bxl1i

也就是把 flmidg1rl做卷积,FFT/NTT即可

【模板】分治 FFT

代码
#include <bits/stdc++.h>
#define endl putchar('\n')
using namespace std;
const int M = 3e6+10;
const int inf = 2147483647;
const int Mod = 998244353,G = 3,Gi = 332748118;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
ll qpow(ll x,ll k,ll mod){
    ll res=1;
    while(k){
        if(k&1) res=res*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return res;
}
inline int add(int x,int y){
    return (x+y)>=Mod?x+y-Mod:(x+y);
}
inline int minu(int x,int y){
    return (x-y)<0?x-y+Mod:(x-y);
}
struct Polygen{
    int r[M],a[M],b[M],lim,l;
    void NTT(int *a,int type){
        for(int i=0;i<lim;i++) if(r[i]>i) swap(a[i],a[r[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            ll wn=qpow(type==1?G:Gi,(Mod-1)/(mid<<1),Mod);
            for(int len=mid<<1,j=0;j<lim;j+=len){
                ll w=1;
                for(int k=0;k<mid;k++,w=w*wn%Mod){
                    ll x=a[j+k],y=1ll*w*a[j+k+mid]%Mod;
                    a[j+k]=add(x,y);a[j+k+mid]=minu(x,y);
                }
            }
        }
    }
    void Mul(int *A,int *B,int *C,int n,int m){
        l=0;lim=1;
        while(lim<=n+m) lim<<=1,l++;
        for(int i=0;i<lim;i++) a[i]=b[i]=0,r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        for(int i=0;i<=n;i++) a[i]=A[i];
        for(int i=0;i<=m;i++) b[i]=B[i];
        NTT(a,1);NTT(b,1);
        for(int i=0;i<lim;i++) a[i]=1ll*a[i]*b[i]%Mod;
        NTT(a,-1);
        ll inv=qpow(lim,Mod-2,Mod);
        for(int i=0;i<lim;i++) C[i]=a[i]*inv%Mod;
    }
}Poly;
int n,f[M],g[M],a[M],b[M],c[M];
void divide(int l,int r){
    if(l==r) return;
    int mid=l+r>>1;
    divide(l,mid);
    for(int i=l;i<=mid;i++) a[i-l]=f[i];
    for(int i=1;i<=r-l;i++) b[i-1]=g[i];
    Poly.Mul(a,b,c,mid-l,r-l-1);
    for(int i=mid+1;i<=r;i++) f[i]=add(f[i],c[i-l-1]);
    divide(mid+1,r);
}
signed main(){
    n=read()-1;
    f[0]=1;
    for(int i=1;i<=n;i++) g[i]=read();
    divide(0,n);
    for(int i=0;i<=n;i++) printf("%d ",f[i]);endl;
    return 0;
}

FWT

用于解决位运算卷积问题

讲解

几个结论:

或卷积:

fwt[a]=merge(fwt[a0],fwt[a0]+fwt[a1])

a=merge(a0,a1a0)

与卷积:

fwt[a]=merge(fwt[a0]+fwt[a1],fwt[a1])

a=merge(a0a1,a1)

异或卷积:

fwt[a]=merge(fwt[a0]+fwt[a1],fwt[a0]fwt[a1])

a=merge(a0+a12,a0a12)

关于异或卷积的证明:

代码
#include <bits/stdc++.h>
#define endl putchar('\n')
using namespace std;
const int M = 1e6+10;
const int inf = 2147483647;
const int Mod = 998244353,G = 3,Gi = 332748118;
const double pi = acos(-1.0);
typedef long long ll;
inline int read(){
    int x=0,f=0;char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;c=getchar();
    }
    do{
        x=(x<<1)+(x<<3)+(c^48);
    }while(isdigit(c=getchar()));
    return f?-x:x;
}
ll qpow(ll x,ll k,ll mod){
    ll res=1;
    while(k){
        if(k&1) res=res*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return res;
}
struct sg{
    int x;
    sg(){x=0;}
    sg(int a){x=a;}
    sg& operator =(int a){return x=a,*this;}
    friend sg operator +(sg a,sg b){return sg(a.x+b.x>=Mod?a.x+b.x-Mod:a.x+b.x);}
    friend sg operator -(sg a,sg b){return sg(a.x-b.x<0?a.x-b.x+Mod:a.x-b.x);}
    friend sg operator *(sg a,sg b){return sg(1ll*a.x*b.x>=Mod?1ll*a.x*b.x%Mod:a.x*b.x);}
    sg& operator +=(sg a){return *this=*this+a,*this;}
    sg& operator *=(sg a){return *this=*this*a,*this;}
}a[M],b[M],A[M],B[M];
int n;
inline void OR(sg *f,sg t=1){
    for(int mid=1,w=2;w<=n;mid=w,w<<=1){
        for(int i=0;i<n;i+=w){
            for(int j=0;j<mid;j++){
                f[i+j+mid]+=(f[i+j]*t);
            }
        }
    }
}
inline void AND(sg *f,sg t=1){
    for(int mid=1,w=2;w<=n;mid=w,w<<=1){
        for(int i=0;i<n;i+=w){
            for(int j=0;j<mid;j++){
                f[i+j]+=(f[i+j+mid]*t);
            }
        }
    }
}
inline void XOR(sg *f,sg t=1){
    for(int mid=1,w=2;w<=n;mid=w,w<<=1){
        for(int i=0;i<n;i+=w){
            for(int j=0;j<mid;j++){
                f[i+j]+=f[i+j+mid];
                f[i+j+mid]=f[i+j]-f[i+j+mid]-f[i+j+mid];
                f[i+j]*=t;f[i+j+mid]*=t;
            }
        }
    }
}
void ctrlv(){
    for(int i=0;i<n;i++) A[i]=a[i],B[i]=b[i];
}
void print(){
    for(int i=0;i<n;i++) printf("%d ",A[i].x);endl;
}
void roll(){
    for(int i=0;i<n;i++) A[i]*=B[i];
}
signed main(){
    n=read();n=1<<n;
    for(int i=0;i<n;i++) a[i]=read();
    for(int i=0;i<n;i++) b[i]=read();
    ctrlv();OR(A);OR(B);roll();OR(A,sg(Mod-1));print();
    ctrlv();AND(A);AND(B);roll();AND(A,sg(Mod-1));print();
    ctrlv();XOR(A);XOR(B);roll();XOR(A,sg(qpow(2,Mod-2,Mod)));print();
    return 0;
}

posted @   blue_tsg  阅读(186)  评论(1编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 上周热点回顾(2.17-2.23)
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
点击右上角即可分享
微信分享提示