学习傅里叶的基本性质及其代码,可以参考大神理解
还有 ACdream 的博客
贴一下NTT的模板:
using namespace std; typedef long long ll; int n; const ll MOD=998244353; const int N = 400005; const int g=3; int len; ll A[N]; long long a[N],b[N],wn[30]; long long q_pow(long long x,long long y,long long P) { long long ans=1; while(y>0) { if(y&1)ans=ans*x%P; x=x*x%P; y>>=1; } return ans; } void init() { for(int i=0;i<21;i++) { int t=1<<i; wn[i]=q_pow(g,(MOD-1)/t,MOD); } } ///雷德算法,2^M=len,将第i位的数与“i的二进制反转之后的位”的数交换 void rader(long long F[],int len) { int j=len/2;///模拟二进制反转进位的的位置 for(int i=1;i<len-1;i++) { if(i<j)swap(F[i],F[j]);///该出手时就出手 int k=len/2; while(j>=k) { j-=k; k>>=1; } if(j<k)j+=k; } } void NTT(long long F[],int len,int t) { int id=0; rader(F,len); for(int h=2;h<=len;h<<=1) { id++; for(int j=0;j<len;j+=h) { long long E=1; for(int k=j;k<j+h/2;k++) { long long u=F[k]; long long v=(E*F[k+h/2])%MOD; F[k]=(u+v)%MOD; F[k+h/2]=((u-v)%MOD+MOD)%MOD; E=(E*wn[id])%MOD; } } } if(t==-1) { for(int i=1;i<len/2;i++)swap(F[i],F[len-i]); long long ni=q_pow(len,MOD-2,MOD); for(int i=0;i<len;i++)F[i]=(F[i]%MOD*ni)%MOD; } } void work()///卷积,点乘,插值 { NTT(a,len,1); NTT(b,len,1); for(int i=0;i<len;i++) a[i]=(a[i]*b[i])%MOD; NTT(a,len,-1); } int main() { #ifndef ONLINE_JUDGE freopen("in.txt","r",stdin); #endif init(); int T_T; scanf("%d",&T_T); for(int kase=1;kase<=T_T;kase++) { sc(n); len=1; while(len<=2*n) len<<=1;
/**
这部分就是你对公式的变形并把他转化为卷积形式,分别存在a[],b[]里
*/
work();
}
return 0;
}
对于傅里叶运用的要点,要认真对待对下面卷积公式的理解
训练:
传送门:hdu 5829 Rikka with Subset
题意:给你一个数组A[],然后计算F[k],F[k]指A[]所有子集中,先对每个子集前k大的数的和,然后再对把所以子集所求值求和
思路:因为我也是初学,所以我找的一个大神的博客学习的http://blog.csdn.net/cqu_hyx/article/details/52194696
对于其中的几点,我做一点补充:
- f[k]计算的是第k大的数对答案的贡献,而比他大的数的贡献没有计算;不过只需要累加分f[1]..f[k-1]就是前k大的数的和了
- 对比以上卷积公式,我们得到的是 a[i]=2n-i/i! , b[i]=A[i]*(i-1)! ;如果我们b[0]=A[1]*0!,那么我们reverse一下,b[i]=b[n-i]了;在我们推导的f[k]公式里,f[k]是与a[i]*b[i+k]相关的,且i最大为(n-k),经过我们对b数组reverse,f[k]是与a[i]*b[(n-k)-i]相关的,这刚好和上面的卷积公式一致,所以我们可以放心ntt了。
- 由卷积公式可知,y[n-k]是a[i]*b[(n-k)-i]的卷积结果,而这个结果存在a[]里;并且f[k]也等于a[i]*b[(n-k)-i]的卷积,所以f[k]与y[n-k]即a[n-k]相关
/************************************************************** Problem:hdu 5829 Rikka with Subset User: youmi Language: C++ Result: Accepted Time:2605MS Memory:11804K ****************************************************************/ //#pragma comment(linker, "/STACK:1024000000,1024000000") //#include<bits/stdc++.h> #include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <map> #include <stack> #include <set> #include <sstream> #include <cmath> #include <queue> #include <deque> #include <string> #include <vector> #define zeros(a) memset(a,0,sizeof(a)) #define ones(a) memset(a,-1,sizeof(a)) #define sc(a) scanf("%d",&a) #define sc2(a,b) scanf("%d%d",&a,&b) #define sc3(a,b,c) scanf("%d%d%d",&a,&b,&c) #define scs(a) scanf("%s",a) #define sclld(a) scanf("%I64d",&a) #define pt(a) printf("%d ",a) #define ptlld(a) printf("%I64d ",a) #define rep(i,from,to) for(int i=from;i<=to;i++) #define irep(i,to,from) for(int i=to;i>=from;i--) #define Max(a,b) ((a)>(b)?(a):(b)) #define Min(a,b) ((a)<(b)?(a):(b)) #define lson (step<<1) #define rson (lson+1) #define eps 1e-6 #define oo 0x3fffffff #define TEST cout<<"*************************"<<endl const double pi=4*atan(1.0); using namespace std; typedef long long ll; int n; const ll MOD=998244353; const int N = 400005; const int g=3; int len; ll inv2; ll A[N]; ll inv[N],fac[N],f[N],bit[N]; long long a[N],b[N],wn[30]; long long q_pow(long long x,long long y,long long P) { long long ans=1; while(y>0) { if(y&1)ans=ans*x%P; x=x*x%P; y>>=1; } return ans; } void init() { for(int i=0;i<21;i++) { int t=1<<i; wn[i]=q_pow(g,(MOD-1)/t,MOD); } inv[0]=inv[1]=1; rep(i,2,1e5) inv[i]=(inv[i-1]*q_pow(i,MOD-2,MOD))%MOD; fac[0]=fac[1]=1; rep(i,2,1e5) fac[i]=(fac[i-1]*i)%MOD; bit[0]=1; rep(i,1,1e5) bit[i]=bit[i-1]*2%MOD; inv2=q_pow(2,MOD-2,MOD); } ///雷德算法,2^M=len,将第i位的数与“i的二进制反转之后的位”的数交换 void rader(long long F[],int len) { int j=len/2;///模拟二进制反转进位的的位置 for(int i=1;i<len-1;i++) { if(i<j)swap(F[i],F[j]);///该出手时就出手 int k=len/2; while(j>=k) { j-=k; k>>=1; } if(j<k)j+=k; } } void NTT(long long F[],int len,int t) { int id=0; rader(F,len); for(int h=2;h<=len;h<<=1) { id++; for(int j=0;j<len;j+=h) { long long E=1; for(int k=j;k<j+h/2;k++) { long long u=F[k]; long long v=(E*F[k+h/2])%MOD; F[k]=(u+v)%MOD; F[k+h/2]=((u-v)%MOD+MOD)%MOD; E=(E*wn[id])%MOD; } } } if(t==-1) { for(int i=1;i<len/2;i++)swap(F[i],F[len-i]); long long ni=q_pow(len,MOD-2,MOD); for(int i=0;i<len;i++)F[i]=(F[i]%MOD*ni)%MOD; } } void work()///卷积,点乘,插值 { NTT(a,len,1); NTT(b,len,1); for(int i=0;i<len;i++) a[i]=(a[i]*b[i])%MOD; NTT(a,len,-1); } int main() { #ifndef ONLINE_JUDGE freopen("in.txt","r",stdin); #endif init(); int T_T; scanf("%d",&T_T); for(int kase=1;kase<=T_T;kase++) { sc(n); len=1; while(len<=2*n) len<<=1; rep(i,1,n) sclld(A[i]); sort(A+1,A+1+n,greater<ll>()); zeros(a);zeros(b); rep(i,0,n-1) a[i]=(bit[n-i]*inv[i])%MOD; rep(i,0,n-1) b[i]=(A[i+1]*fac[i])%MOD; reverse(b,b+n); work(); ll r=inv2; f[0]=0; rep(i,1,n) { f[i]=a[n-i]*inv[i-1]%MOD*r%MOD; r=r*inv2%MOD; f[i]=(f[i]+f[i-1])%MOD; printf("%I64d ",f[i]); } puts(""); } }