FFT/NTT/FWT专题

50 阅读 0 评论 33 点赞

我是靠谱客的博主矮小白昼，最近开发中收集的这篇文章主要介绍FFT/NTT/FWT专题，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

FFT可快速计算卷积，卷积为两个生成函数的乘积，生成函数为 $g(n)=sum_{i=0}^{n}f(i)x^i$

FFT 卷积生成函数的关系

多项式之快速傅里叶变换(FFT)/数论变换(NTT)/例题与常用套路【入门】

HDU1402 A * B Problem Plus

大数乘法

生成函数: $g(n)=sum_{i=0}^{n}f(i)*10^i$ ,设两个大数位： $a(l1)$ 、 $b(l2)$ ，L1/L2为两个大数的长度，

$gsum(n)=g1(n)*g2(n)=sum_{i+j=n}g1(i)g2(j)=sum_{i=0}^{n}g1(i)g2(n-i)$

利用FFT求两个函数的卷积，将所得 $ans(n)$ 改为十进制整数即为所求答案

#include<bits/stdc++.h>
using namespace std;
const int MAX=2e5+5;
const double PI=acos(-1.0);
typedef complex<double> C;
C a[MAX],b[MAX];
char s1[MAX],s2[MAX];
int ans[MAX];
int L;
void rader(C *c_)
{
int i,j,k;
for(i=1,j=L/2;i<L-1;++i)
{
if(i<j) swap(c_[i],c_[j]);
k=L/2;
while(j>=k){j-=k;k/=2;}
if(j<k) j+=k;
}
}
void fft(C *c_,int v)
{
rader(c_);
for(int i=2;i<=L;i<<=1)
{
C wn(cos(-2.0*v*PI/i),sin(-2.0*v*PI/i));
for(int j=0;j<L;j+=i)
{
C w(1,0);
for(int k=j;k<j+i/2;++k)
{
C u=c_[k];
C t=w*c_[k+i/2];
c_[k]=u+t;
c_[k+i/2]=u-t;
w=w*wn;
}
}
}
if(v==-1) for(int i=0;i<L;++i) c_[i].real(c_[i].real()/L);
}
int main()
{
while(~scanf("%s%s",s1,s2))
{
int L1=strlen(s1),L2=strlen(s2);
L=1;
while(L<L1<<1||L<L2<<1) L<<=1;
for(int i=0;i<L1;++i) a[i]=C(s1[L1-i-1]-'0',0);
for(int i=L1;i<L;++i) a[i]=C(0,0);
for(int i=0;i<L2;++i) b[i]=C(s2[L2-i-1]-'0',0);
for(int i=L2;i<L;++i) b[i]=C(0,0);
fft(a,1);fft(b,1);
for(int i=0;i<L;++i) a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<L;++i) ans[i]=(int)(a[i].real()+0.5);
for(int i=0;i<L;++i)
{
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
L=L1+L2-1;
while(ans[L]<=0&&L>0) --L;
for(int i=L;i>=0;--i)
printf("%d",ans[i]);
printf("n");
}
return 0;
}

SPOJ - TSUMSP Triple Sums

N个数中取3个，求和的值和其方案数

将N个数复制为3堆，从每堆中取出1个数,FFT+组合计数+容斥：

设一堆固定的生成函数为 $a1(n)$ ,两堆固定的生成函数为 $a2(n)$ ,三堆固定的生成函数为 $a3(n)$ ,答案的生成函数为 $ans(n)$

$ans(n)=frac{a1(n)^3-a2(n)a1(n)*3+a3(n)*2}{6}$

#include<bits/stdc++.h>
using namespace std;
const int MAX=1<<17;
const double PI=acos(-1.0);
typedef complex<double> C;
map<int,int>mp;
map<int,int>::iterator p;
C a1[MAX],a2[MAX],a3[MAX],ans[MAX];
int n,L,cnt[MAX];
void rader(C *c_)
{
int i,j,k;
for(i=1,j=L/2;i<L-1;++i)
{
if(i<j) swap(c_[i],c_[j]);
k=L/2;
while(j>=k){j-=k;k/=2;}
if(j<k) j+=k;
}
}
void fft(C *c_,int v)
{
rader(c_);
for(int i=2;i<=L;i<<=1)
{
C wn(cos(-2.0*v*PI/i),sin(-2.0*v*PI/i));
for(int j=0;j<L;j+=i)
{
C w(1,0);
for(int k=j;k<j+i/2;++k)
{
C u=c_[k];
C t=w*c_[k+i/2];
c_[k]=u+t;
c_[k+i/2]=u-t;
w=w*wn;
}
}
}
if(v==-1) for(int i=0;i<L;++i) c_[i].real(c_[i].real()/L);
}
int main()
{
int x;
scanf("%d",&n);
mp.clear();
for(int i=0;i<n;++i)
{
scanf("%d",&x);x+=20000;
if(mp[x]) ++mp[x];
else mp[x]=1;
}
for(p=mp.begin();p!=mp.end();++p) a1[p->first]=a2[p->first*2]=a3[p->first*3]=C(p->second,0);
L=MAX;
fft(a1,1);fft(a2,1);fft(a3,1);
for(int i=0;i<L;++i) ans[i]=a1[i]*a1[i]*a1[i]-a1[i]*a2[i]*3.0+a3[i]*2.0;
fft(ans,-1);
for(int i=0;i<L;++i)
{
long long tmp=(long long)(ans[i].real()+0.5)/6;
if(tmp) printf("%d : %I64dn",i-60000,tmp);
}
return 0;
}

HDU5829 Rikka with Subset

一个集合的价值是前K大的数之和，有N个数，求其各集合的价值之和

将N个数从大到小排列，所有树作为第 $k$ 大的贡献之和为： $f(k)=sum_{i=k}^{n}c(i-1,k-1)2^{n-i}c[i]$

即从前 $i$ 个数中选取 $k-1$ 个，使 $c[i]$ 成为第 $k$ 大，因为其作为第 $j$ 大 $(j<k)$ 大仍有贡献，所以取这 $k-1$ 个数的所有子集 $2^{k-1}$ 个

将该式化为卷积形式：

$f[k]=ff[n-k]=frac{1}{2^k*(k-1)!}sum_{i=0}^{n-k}frac{2^{n-i}}{i!}(i+k-1)!c[i+k]$ , $m=n-k$

$f[k]=frac{1}{2^k*(k-1)!}ff[m]=sum_{i=0}^{m}frac{2^{n-i}}{i!}(i+k-1)!c[i+k]$

$a[i]=frac{2^{n-i}}{i!}$ , $b[i+k]=(n-i-k-1)!c[n-i-k]->b[i]=(n-i-1)!c[n-i]$

$f[k]=frac{1}{2^k*(k-1)!}ff[m<=n]=sum_{i=0}^{m<=n}a[i]b[i]$

求其前缀和

#include<bits/stdc++.h>
using namespace std;
const int MAX=2e5+5;
const long long MOD=998244353;
const int g=3;
long long fact[MAX],inv[MAX],bit[MAX],a[MAX<<1],b[MAX<<1],c[MAX],ni;
int t,n,L;
long long pow(long long a,int k,long long mod)//快速幂取模
{
long long b=1;
while(k)
{
if(k&1) b=b*a%mod;
a=a*a%mod;
k>>=1;
}
return b;
}
void init()
{
fact[0]=inv[0]=bit[0]=1;
for(int i=1;i<MAX;i++)
{
fact[i]=fact[i-1]*i%MOD;
inv[i]=inv[i-1]*pow(i,MOD-2,MOD)%MOD;//利用欧拉函数的积性
bit[i]=bit[i-1]*2%MOD;
}
}
void rader(long long *f_)
{
int i,j,k;
for(i=1,j=L/2;i<L-1;++i)
{
if(i<j) swap(f_[i],f_[j]);
k=L/2;
while(j>=k){j-=k;k>>=1;}
if(j<k) j+=k;
}
}
void ntt(long long *f_,int t)
{
rader(f_);
for(int i=2;i<=L;i<<=1)
{
long long wn=pow(g,(MOD-1)/i,MOD);
if(t==-1) wn=pow(wn,MOD-2,MOD);
for(int j=0;j<L;j+=i)
{
long long e=1;
for(int k=j;k<j+i/2;++k)
{
long long u=f_[k];
long long v=e*f_[k+i/2]%MOD;
f_[k]=(u+v)%MOD;
f_[k+i/2]=(u-v+MOD)%MOD;
e=e*wn%MOD;
}
}
}
if(t==-1) for(int i=0;i<L;++i) f_[i]=f_[i]*ni%MOD;
}
int main()
{
init();
scanf("%d",&t);
while(t--)
{
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
scanf("%d",&n);
L=1;
while(L<n<<1) L<<=1;
ni=pow(L,MOD-2,MOD);
for(int i=1;i<=n;++i) scanf("%I64d",&c[i]);
sort(c+1,c+n+1,greater<long long>());
for(int i=0;i<n;++i)
{
a[i]=bit[n-i]*inv[i]%MOD;
b[i]=c[n-i]*fact[n-i-1]%MOD;
}
ntt(a,1);ntt(b,1);
for(int i=0;i<L;++i) a[i]=a[i]*b[i]%MOD;
ntt(a,-1);
long long r=inv[2],ans=0;
for(int i=1;i<=n;++i)
{
ans=(ans+a[n-i]*inv[i-1]%MOD*r%MOD)%MOD;
r=r*inv[2]%MOD;
printf("%I64d ",ans);
}
printf("n");
}
return 0;
}

FWT

void FWT_or(int *a,int opt)
{
for(int i=1;i<N;i<<=1)
for(int p=i<<1,j=0;j<N;j+=p)
for(int k=0;k<i;++k)
if(opt==1)a[i+j+k]=(a[j+k]+a[i+j+k])%MOD;
else a[i+j+k]=(a[i+j+k]+MOD-a[j+k])%MOD;
}
void FWT_and(int *a,int opt)
{
for(int i=1;i<N;i<<=1)
for(int p=i<<1,j=0;j<N;j+=p)
for(int k=0;k<i;++k)
if(opt==1)a[j+k]=(a[j+k]+a[i+j+k])%MOD;
else a[j+k]=(a[j+k]+MOD-a[i+j+k])%MOD;
}
void FWT_xor(int *a,int opt)
{
for(int i=1;i<N;i<<=1)
for(int p=i<<1,j=0;j<N;j+=p)
for(int k=0;k<i;++k)
{
int X=a[j+k],Y=a[i+j+k];
a[j+k]=(X+Y)%MOD;a[i+j+k]=(X+MOD-Y)%MOD;
if(opt==-1)a[j+k]=1ll*a[j+k]*inv2%MOD,a[i+j+k]=1ll*a[i+j+k]*inv2%MOD;
}
}