学习笔记:多项式变换
多项式学习笔记
卷积形式:
一般卷积形式
当\(\circ\)为 \(+\)是为常见的多项式加法卷积,\(\times\)有时也能转为加法,位运算是FWT,整除是狄利克雷卷积
加法卷积常用形式:
卷积满足交换律,结合律,分配律
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;
}
例题
给出 \(n\) 个数 \(q_1,q_2, \dots q_n\),定义
对 \(1 \leq i \leq n\),求 \(E_i\) 的值。
把式子换成卷积形式,常用技巧有翻转序列,定义初值以扩大枚举范围等
令 \(g_i=\frac{1}{i^2},g_0=0,q_0=0\),\(q'_i=q_{n-i}\),\(t=n-j\)
然后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\),使得$a^r \equiv 1 \pmod m $的最小正整数 \(r\) 称为\(a\)对模\(m\)的阶,记作\(\delta_m(a)\)
性质:\(\delta_a(m) \; | \; \phi(m)\) (欧拉定理证)
原根:若有 \(m>1\)且\(\delta_p(m)=\phi(p)\),则称 \(m\) 是 \(p\) 的原根
性质:
1.一个正整数 \(m\) 有原根的充要条件是\(m=2,4,p^e,2p^e\)
2.有原根的正整数 \(m\) 有 \(\phi(\phi(m))\) 个原根
3.若 \(g\) 是 \(m\) 的一个原根,则
互不相同,且值域为1到\(\phi(m)\)
用处:
1.原根同样满足单位根的一些性质,可以代替单位根以实现NTT,用于保证精度和取模运算
2.用性质3可以把数 \(x\) 转成 \(log_gx (\; \text{mod} \; m)\),把乘法运算变成加法运算,注意如果卷积是取模意义下的根据欧拉定理要模\(\phi(m)\)
即把
变成
其中\(X=log_gx \; \text{mod} \; p,I=log_gi \; \text{mod} \; p,J=log_gj\; \text{mod} \; p\)
求法:
对\(\phi(m)\)分解质因数,枚举\(g\),若恒满足
则 \(g\) 是 \(m\) 的原根
NTT实现
把单位根换成 \(g^{(p-1)/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;
}
例题
设 \(h_t(x)\) 表示从集合中选 \(t\) 个数加和为 \(i\) 的方案数,有
这个过程可以用倍增加NTT的方法快速求解
然后考虑把乘法变成加法,应用上面提到的的技巧,把 \(x\) 变成 \(log_gx\)即可
代码
#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
用于求解:
给定序列 \(g_{1\dots n - 1}\),求序列 \(f_{0\dots n - 1}\)。
其中 \(f_i=\sum_{j=1}^if_{i-j}g_j\),边界为 \(f_0=1\)。
考虑分治,先求出 \([l,mid]\) 的 \(f_x\) 值,把左边对右边的贡献加上,再算右边的值
考虑怎么计算中间贡献,设 $f_l \sim f_{mid} $ 对 \(f_x\) 的贡献为 \(w_x\) 有
令 \(f_{mid+1 \sim r} = 0\),有
设 \(a_i=f_{i+l},b_i=g_{i+1}\),有
也就是把 \(f_{l \sim mid}\)与 \(g_{1 \sim r-l}\)做卷积,FFT/NTT即可
代码
#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
用于解决位运算卷积问题
几个结论:
或卷积:
与卷积:
异或卷积:
关于异或卷积的证明:
代码
#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;
}