生成函数+FFT
之前做了几道感觉推完式子都是FFT的纯板子题,然后做这道题想到了容斥,但是不知道应该怎么上FFT,因为这个并不是纯卷积,或者说纯卷积十分累赘复杂。
FFT其实只是一个多项式点值表达式和系数表达式之间转换的一个工具,我们要做的就是理清楚各个多项式之间准确的关系,或者把卷积转化成原始的for循环理解
然而生成函数的原理是不能变的,生成函数i次方的系数就是代表对应权值为i的方案数,否则权值不对应,乘起来就乱套了。
所以我们设三个生成函数,A表示选一个权值为x的斧子,B表示选两个权值为x/2的斧子,C表示选三个权值为x/3的斧子。
那么选一个的方案数就是A,选两个的就是(A*A-B)/2,这里我们考虑原始的for循环,A*A是两层从1到n的循环,然而我们要去掉i=j的情况,就是B,然后还应该除去组合顺序,所以就是上式,选三个的也差不多,首先有三层循环,然后我们应该去掉i=j||j=k||i=k的情况,就是另外两层循环:A*B*3,但是这里i=j=k的情况多减了两次,所以我们要加上C,最后在除去组合顺序,就是(A*A*A-A*B*3+C*2)/6。
说的太多了,还是因为自己不明白。233
1 #include <cstdio> 2 #include <cstring> 3 #include <iostream> 4 #include <algorithm> 5 #include <cmath> 6 #define MAXN 133333 7 using namespace std; 8 const double pi=acos(-1.0); 9 struct cmx{ 10 double x,y; 11 cmx(){} 12 cmx(double a,double b){x=a,y=b;} 13 cmx operator + (cmx a){return cmx(x+a.x,y+a.y);} 14 cmx operator - (cmx a){return cmx(x-a.x,y-a.y);} 15 cmx operator * (cmx a){return cmx(x*a.x-y*a.y,x*a.y+y*a.x);} 16 cmx operator / (double a){return cmx(x/a,y/a);} 17 cmx operator * (double a){return cmx(x*a,y*a);} 18 }A[MAXN],B[MAXN],C[MAXN],ans[MAXN]; 19 int n,N,rev[MAXN],maxn; 20 void fft(cmx *a,int o){ 21 cmx dan,now,t; 22 register int i,j,k; 23 for(i=0;i<N;i++)if(i<rev[i])swap(a[i],a[rev[i]]); 24 for(k=2;k<=N;k<<=1){ 25 dan=cmx(cos(2*pi/k),o*sin(2*pi/k)); 26 for(j=0;j<N;j+=k){ 27 now=cmx(1,0); 28 for(i=0;i<(k>>1);i++,now=now*dan){ 29 t=now*a[i+j+(k>>1)]; 30 a[i+j+(k>>1)]=a[i+j]-t; 31 a[i+j]=a[i+j]+t; 32 } 33 } 34 } 35 if(o==-1)for(i=0;i<N;i++)a[i].x/=N; 36 } 37 int main(){ 38 scanf("%d",&n); 39 register int i; 40 register long long now; 41 for(int i=1,x;i<=n;i++){ 42 scanf("%d",&x); 43 A[x].x+=1; 44 B[x*2].x+=1; 45 C[x*3].x+=1; 46 maxn=max(maxn,3*x); 47 } 48 for(N=1;N<maxn;N<<=1); 49 for(i=0;i<N;i++){ 50 if(i&1)rev[i]=(rev[i>>1]>>1)|(N>>1); 51 else rev[i]=rev[i>>1]>>1; 52 } 53 fft(A,1);fft(B,1);fft(C,1); 54 for(i=0;i<N;i++) 55 ans[i]=A[i]+(A[i]*A[i]-B[i])/2+(A[i]*A[i]*A[i]-A[i]*B[i]*3+C[i]*2)/6; 56 fft(ans,-1); 57 for(i=0;i<maxn;i++){ 58 now=1ll*round(ans[i].x); 59 if(now>0)printf("%d %lld ",i,now); 60 } 61 return 0; 62 }