[BZOJ3771]Triple（快速傅里叶变换FFT）

【题意】
n（n<=500000）个物品，可以用1/2/3个不同的物品组成不同的价值，求每种价值有多少种方案(顺序不同算一种)

【输入】
第一行是整数N，表示有N个物品
接下来n行升序输入N个数字Ai，表示每个物品的价值。

【输出】
若干行，按升序对于所有可能的总价值输出一行x y，x为价值，y为方案数。

【样例输入】
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

引用ATP酱的话：

要求所有组合的方案数，并且还要求具体组合出来的方案，可以想到利用生成函数

先简单说一下生成函数，生成函数其实就是一个多项式，对于这一题我们可以把多项式中x的指数看成斧头的价格，x的系数看成斧头的数量，那么这样有什么用呢，我们举一个例子来说明，比如说我有2个价值为5的斧头，1个价值为7的斧头
那么他们的多项式分别表示为$2x^5$和$x^7$那么把两个多项式乘起来，得到$2x^{12}$，因为乘法的时候指数相加，系数相乘，那么得到的这个多项式的系数就是选他们两个他们的方案数，指数就是选他们两个的损失。

我们设拿一个斧头的生成函数为$F(x)$
根据生成函数的定义，那么拿两个，三个的生成函数就为$F^2(x),F^3(x)$
那么我们是不是可以理所当然的认为答案就是
$F(x)+F^2(x)+F^3(x)$呢
当然不是！因为这里面是有重复的（重复的拿一个方案和一个斧头的情况都有）！
比如说：
在一个$F^2(x)$生成函数中：
$(x+x^2)^2$
展开得：
$x^2+2x^3+x^4$
我们可以发现在这样展开的多项式里面有一个$x^2$，一个$x^4$和两个$x^3$
他们分别是选两个斧头1，选两个斧头2，和一个斧头1一个斧头2得出来的，
但是根据题意，我们不能选两个一样的斧头，而且选斧头的顺序不同只能算一种，也就是说在我们展开的多项式里面只有一个$x^3$是有用的！
那么我们怎么去掉其他的无用状态呢？那么我们就另$B(x)$、$C(x)$为一种斧头选两次和一种斧头选三次的生成函数，另$F^2(x)-B(x)$就可以去掉重复选一个的状态，至于选斧头的顺序呢，我们可以除一个2!就可以啦。
还有一个问题，怎么实现$B(x)$、$C(x)$呢，hin简单，我们在读入斧头的权值q的时候把q*2,q*3项（也就是指数为q*2,q*3的项）的系数+1就可以啦（可以理解为一次拿了双倍和三倍）
对于三个的也是差不多的，不过要注意容斥噢。

好啦还是讲一下容斥吧

三个叠起来的圈就是$F(x)^3$，重复加了三个绿色的部分，那每一个绿色的部分（也就是一片绿色+那块红色）就是$3F(x)B(x)$，重复减了两个红色那就要加上$2C(x)$

那么说了这么多答案就出来啦：

$Ans=\frac{F(x)}{1!}+\frac{F^2(x)-B(x)}{2!}+\frac{F(x)^3-3F(x)B(x)+2C(x)}{3!}$

再次提醒精度的问题噢

code:

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<iostream>
using namespace std;
const int maxn=100010;
const double pi=acos(-1.0);
struct complex
{
    double r,i;
    complex(){}
    complex(double _r,double _i){r=_r,i=_i;}
    friend complex operator + (const complex &x,const complex &y){return complex(x.r+y.r,x.i+y.i);}
    friend complex operator - (const complex &x,const complex &y){return complex(x.r-y.r,x.i-y.i);}
    friend complex operator * (const complex &x,const complex &y){return complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
    friend complex operator * (const complex &x,const double &y){return complex(x.r*y,x.i*y);}
    friend complex operator / (const complex &x,const double &y){return complex(x.r/y,x.i/y);}
}a[maxn*4],b[maxn*4],c[maxn*4],ans[maxn*4];
double h1[maxn*4],h2[maxn*4],q[maxn*4];
int R[maxn*4];
void fft(complex *y,int len,int on)
{
    for(int i=0;i<len;i++)if(i<R[i])swap(y[i],y[R[i]]);
    for(int i=1;i<len;i<<=1) 
    {
        complex wn(cos(pi/i),sin(on*pi/i)); 
        for(int j=0;j<len;j+=(i<<1)) 
        {
            complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn) 
            {
                complex u=y[j+k]; 
                complex v=w*y[j+k+i];
                y[j+k]=u+v;
                y[j+k+i]=u-v;
            }
        }
    }
    if(on==-1)for(int i=0;i<len;i++)y[i].r/=len; 
}
int main()
{
    int n;scanf("%d",&n);
    int maxx=-9999999;
    int N=n;
    for(int i=1;i<=N;i++)
    {
        int q;scanf("%d",&q);
        a[q].r+=1.0;
        b[q*2].r+=1.0;
        c[q*3].r+=1.0;
        n=max(n,q*3);
    }
    int m=n*2;
    int L=0;
    for(n=1;n<=m;n*=2) L++;
    for(int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|(i&1)<<(L-1);
    fft(a,n,1);fft(b,n,1);fft(c,n,1);
    for(int i=0;i<n;i++)
    {
        ans[i]=a[i]+(a[i]*a[i]-b[i])/2.0+(a[i]*a[i]*a[i]-a[i]*b[i]*3.0+c[i]*2.0)/6.0;
    }
    fft(ans,n,-1);
    for(int i=0;i<n;i++)
    {
        if(int(ans[i].r+0.5))
        {
            printf("%d %d\n",i,int(ans[i].r+0.5));
        }
    }
    return 0;
}