BZOJ 3771 Triple

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

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

数学语言表达：
给出一个有$n$个数的数列$A$，求集合
{ A i , A i + A j , A i + A j + A k ∣ 1 ≤ i < j < k ≤ n } \{A_i,A_i+A_j,A_i+A_j+A_k\quad|\quad1\le i<j<k\le n\} {Ai​,Ai​+Aj​,Ai​+Aj​+Ak​∣1≤i<j<k≤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$如下表：

0	1	2	3	4	5	6	7	8	9
0	1	1	0	1	0	0	0	0	0

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

	0	1($A_i$)($A_j$)	2($A_i+A_j$)	3	4	5	6	7	8	9
$f$	0	1	1	0	1	0	0	0	0	0
$g$	0	0	1	0	0	0	0	0	0	0

	0	1($A_i$)	2($A_j$)	3($A_i+A_j$)	4	5	6	7	8	9
$f$	0	1	1	0	1	0	0	0	0	0
$g$	0	0	1	1	0	0	0	0	0	0

	0	1($A_i$)	2	3	4($A_j$)	5($A_i+A_j$)	6	7	8	9
$f$	0	1	1	0	1	0	0	0	0	0
$g$	0	0	1	1	0	1	0	0	0	0

	0	1($A_j$)	2($A_i$)	3($A_i+A_j$)	4	5	6	7	8	9
$f$	0	1	1	0	1	0	0	0	0	0
$g$	0	0	1	2	0	1	0	0	0	0

	0	1	2($A_i$)($A_j$)	3	4($A_i+A_j$)	5	6	7	8	9
$f$	0	1	1	0	1	0	0	0	0	0
$g$	0	0	1	2	1	1	0	0	0	0

	0	1	2($A_i$)	3	4($A_j$)	5	6($A_i+A_j$)	7	8	9
$f$	0	1	1	0	1	0	0	0	0	0
$g$	0	0	1	2	1	1	1	0	0	0

……直到最后， g = { 1 , 0 , 2 , 2 , 1 , 2 , 1 , 0 , 0 } g=\{1,0,2,2,1,2,1,0,0\} 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,j\},\{j,i,j\},\{j,j,i\} {i,j,j},{j,i,j},{j,j,i}三种情况，但是其中包括了$i=j$，即三个数相等的情况，于是加回两个$I$，再为了消除排列先后顺序带来的多余项除以6

大概流程如下：

1. $G_{A_i}++,H_{2A_i}++,I_{3A_i}++$
2. 对$G,H,I$进行DFT变换
3. $D=G_i+\frac{{G_i}^2-H_i}{2}+\frac{{G_i}^3-3GH+2I}{6}$
4. 对$D$进行IDFT变换
5. 检测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;
}