BZOJ 3771 Triple

想了我好久，然后发现我好像还是没有学到FFT的精髓…… 万物皆可FFT!

由于是权限题，我放一下大意： 给出一个数列，让你从中间按顺序挑出1~3个数，求它们的和可能是多少，并且从小到大输出每种可能的方案数。

数学语言表达： 给出一个有$n$个数的数列$A$，求集合 $\{A_i,A_i+A_j,A_i+A_j+A_k|1\le i<j<k\le n\}$ 并对于每个数，输出组成它的可能方案数

贴个样例就完事儿了

样例输入：

4 4 5 6 8

样例输出：

4 1 5 1 6 1 8 1 9 1 10 1 11 1 12 1 13 1 14 1 15 1 17 1 18 1 19 1

记得之前OZY师兄给我们讲生成函数，就讲了类似的一道题 一眼看下去是个背包嘛，但是数据范围大而且有一些奇奇怪怪的限制，还要输出方案数 不管怎么样先打个暴力再说嘛 首先抛开大数据和个数的限制范围，我们来看固定使用两个物品的解决方案，并且包括重复使用一个物品的情况 转移方程就是这样的啦： $F_k=\sum _{i=1}^n\sum _{j=i}^n(A_i+A_j=k)?1:0$ 我们再回顾一下卷积： $F_k=\sum _{i=1}^k{A}_i{B}_{k-i}$ 乱搞一下可以变成这样： $F_k=\sum _{i=1}^n\sum _{j=i}^n(i+j=k)?A_i{B}_j:0$ 强行转化成卷积啦

实际上我们可以按照这个思路去做啊： 建立一个多项式$G(x)$，使$\forall i\in [1,n]$，$G_{A_i}++$（相当于背包中$f[a[i]]++$） 将$G$自己跟自己卷一下（卷积讲起来好捞啊，其实就是求一个平方啦），求解过程$G_{A_i+A_j}+=G_{A_i}\ast G_{A_j}$与$g[a[i]+a[j]]+=f[a[i]]\ast f[a[j]]$对应 写一个表出来可能更好理解：

对于本问题中使用两个物品的部分 如果允许出现$i=j$并且$(i,j)$与$(j,i)$不等价 设数列$A$为$1,2,4$ 暴力中的$f$如下表：

01234567890110100000

求解答案表$g$的过程如下：

01($A_i$)($A_j$)2($A_i+A_j$)3456789$f$0110100000$g$0010000000 01($A_i$)2($A_j$)3($A_i+A_j$)456789$f$0110100000$g$0011000000 01($A_i$)234($A_j$)5($A_i+A_j$)6789$f$0110100000$g$0011010000 01($A_j$)2($A_i$)3($A_i+A_j$)456789$f$0110100000$g$0012010000 012($A_i$)($A_j$)34($A_i+A_j$)56789$f$0110100000$g$0012110000 012($A_i$)34($A_j$)56($A_i+A_j$)789$f$0110100000$g$0012111000

……直到最后，$g=\{1,0,2,2,1,2,1,0,0\}$ 恰是我们卷积时进行的一模一样的步骤，一模一样的结果

扯了这么多，终于可以用FFT对它进行加速了 对于3个物品的情况，我们只需要求$G$的三次方即可。

接下来就是解决重复的问题了 如果$D$是答案，在不解决重复问题的情况下$D=G^3+G^2+G$ 为了解决重复问题，我们引入$H(x)$和$I(x)$ 其中，$G_i=H_{2i}=I_{3i}$ $H,I$分别代表重复取两个相同数和三个相同数的情况。 当取一个数时，没有重复的情况。 当取两个数时，我们为消除$i=j$的情况减去一个$H$，为了消除$j<i$的情况除以2 当取三个数时，我们减去$3GH$，分别对应$\{i,j,j\},\{j,i,j\},\{j,j,i\}$三种情况，但是其中包括了$i=j$，即三个数相等的情况，于是加回两个$I$，再为了消除排列先后顺序带来的多余项除以6

大概流程如下：

$G_{A_i}++,H_{2A_i}++,I_{3A_i}++$对$G,H,I$进行DFT变换$D=G_i+\frac{{G_i}^2-H_i}{2}+\frac{{G_i}^3-3GH+2I}{6}$对$D$进行IDFT变换检测D的各项，若不为零则输出

代码如下：

#include<cstdio> #include<cmath> #include<algorithm> using namespace std; const double pi=acos(-1); struct cp { double x,i; cp(double _x=0,double _i=0):x(_x),i(_i) {} cp(const cp &b):x(b.x),i(b.i) {} cp operator +(cp b)const{return cp(x+b.x,i+b.i);} cp operator -(cp b)const{return cp(x-b.x,i-b.i);} cp operator *(cp b)const{return cp(x*b.x-i*b.i,x*b.i+i*b.x);} cp operator +(double b)const{return (*this)+cp(b);} cp operator -(double b)const{return (*this)-cp(b);} cp operator *(double b)const{return (*this)*cp(b);} cp& operator *=(cp b){return (*this)=(*this)*b;} };//复数 cp operator *(int a,cp b) {return b*a;} int R[1<<17]; int x[41000]; cp a[1<<17|1],b[1<<17|1],c[1<<17|1],ans[1<<17|1]; void fft(cp arr[],int len,int inv) { for(int n=2;n<=len;n<<=1) { const cp ur(cos(2*pi/n),sin(2*inv*pi/n)); const int mid=n>>1; for(cp *a=arr;a<arr+len-1;a+=n) { cp r(1),t; for(int i=0;i<mid;i++,r*=ur)//被迫从0开始————存在等于零的情况 { t=a[i+mid]*r; a[i+mid]=a[i]-t; a[i]=a[i]+t; } } } } long long int k[1<<17]; int main() { int n,m=0;scanf("%d",&n); for(int i=1;i<=n;i++) { scanf("%d",x+i); a[x[i]].x++; b[2*x[i]].x++; c[3*x[i]].x++; m=max(m,3*x[i]); } for(n=1;n<=m;n<<=1); for(int i=1;i<n;i++) R[i]=R[i>>1]>>1|(i&1?n>>1:0); for(int i=0;i<n;i++) if(i<R[i])swap(a[i],a[R[i]]); fft(a,n,1); for(int i=0;i<n;i++) if(i<R[i])swap(b[i],b[R[i]]); fft(b,n,1); for(int i=0;i<n;i++) if(i<R[i])swap(c[i],c[R[i]]); fft(c,n,1); for(int i=0;i<=n;i++) ans[i]=a[i]+(a[i]*a[i]-b[i])*(1.0/2)+(a[i]*a[i]*a[i]-3.0*a[i]*b[i]+2.0*c[i])*(1.0/6); for(int i=0;i<=n;i++) if(i<R[i])swap(ans[i],ans[R[i]]); fft(ans,n,-1); for(int i=0;i<=n;i++) { k[i]=ans[i].x/n+0.5; if(k[i])printf("%d %lld\n",i,k[i]); } return 0; }