学习笔记:多项式变换
多项式学习笔记
卷积形式:
一般卷积形式
当为 是为常见的多项式加法卷积,有时也能转为加法,位运算是FWT,整除是狄利克雷卷积
加法卷积常用形式:
卷积满足交换律,结合律,分配律
FFT
利用单位根和复数的性质分治的去把多项式系数表达式转成点值表达式, 把对应点值乘起来后再转成系数表达式
复杂度
递归形式
/*
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;
}
例题
给出 个数 ,定义
对 ,求 的值。
把式子换成卷积形式,常用技巧有翻转序列,定义初值以扩大枚举范围等
令 ,,
然后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
原根:
阶:
设,使得的最小正整数 称为对模的阶,记作
性质: (欧拉定理证)
原根:若有 且,则称 是 的原根
性质:
1.一个正整数 有原根的充要条件是
2.有原根的正整数 有 个原根
3.若 是 的一个原根,则
互不相同,且值域为1到
用处:
1.原根同样满足单位根的一些性质,可以代替单位根以实现NTT,用于保证精度和取模运算
2.用性质3可以把数 转成 ,把乘法运算变成加法运算,注意如果卷积是取模意义下的根据欧拉定理要模
即把
变成
其中
求法:
对分解质因数,枚举,若恒满足
则 是 的原根
NTT实现
把单位根换成 即可,记得取模
封装带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;
}
例题
设 表示从集合中选 个数加和为 的方案数,有
这个过程可以用倍增加NTT的方法快速求解
然后考虑把乘法变成加法,应用上面提到的的技巧,把 变成 即可
代码
#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
用于求解:
给定序列 ,求序列 。
其中 ,边界为 。
考虑分治,先求出 的 值,把左边对右边的贡献加上,再算右边的值
考虑怎么计算中间贡献,设 对 的贡献为 有
令 ,有
设 ,有
也就是把 与 做卷积,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;
}
This is your dream.Anything you can dream,you can do it now.
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 上周热点回顾(2.17-2.23)
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)