BZOJ3771 Triple（FFT+容斥原理）

　　思路比较直观。设A(x)=Σx^a_i。先把只选一种的统计进去。然后考虑选两种，这个直接A(x)自己卷起来就好了，要去掉选同一种的情况然后除以2。现在得到了选两种的每种权值的方案数，再把这个卷上A(x)。得到这个后考虑去重，其中重复的就是选了两个相同的和另外一个，那么再把选两个相同的生成函数搞出来卷上A，减掉选三个相同的。把这个东西减掉之后再除以3。说了半天也不知道在说啥，总之是容斥原理很基础的应用。

　　有些卡精度，用long double才过，可能是我写丑了。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
#define N 270000
#define double long double
const double PI=3.14159265358979324;
struct complex
{
    double x,y;
    complex operator +(const complex&a) const
    {
        return (complex){x+a.x,y+a.y};
    }
    complex operator -(const complex&a) const
    {
        return (complex){x-a.x,y-a.y};
    }
    complex operator *(const complex&a) const
    {
        return (complex){x*a.x-y*a.y,x*a.y+y*a.x};
    }
}w[N],v[N],u[N];
int n,m,t,a[N],r[N];
long long f[N];
void DFT(int n,complex *a,int p)
{
    for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
    for (int i=2;i<=n;i<<=1)
    {
        complex wn=(complex){cos(2*PI/i),p*sin(2*PI/i)};
        for (int j=0;j<n;j+=i)
        {
            complex w=(complex){1,0};
            for (int k=j;k<j+(i>>1);k++,w=w*wn)
            {
                complex x=a[k],y=w*a[k+(i>>1)];
                a[k]=x+y,a[k+(i>>1)]=x-y;
            }
        }
    }
}
void mul(int n,complex *a,complex *b)
{
    for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1)*(n>>1);
    for (int i=0;i<n;i++) a[i].y=a[i].x-b[i].x,a[i].x=a[i].x+b[i].x;
    DFT(n,a,1);
    for (int i=0;i<n;i++) a[i]=a[i]*a[i];
    DFT(n,a,-1);
    for (int i=0;i<n;i++) a[i].x=a[i].x/n/4;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj3771.in","r",stdin);
    freopen("bzoj3771.out","w",stdout);
    const char LL[]="%d %I64d
";
#else
    const char LL[]="%d %lld
";
#endif
    n=read();
    for (int i=1;i<=n;i++) 
    {
        int x=read();
        m=max(m,x);
        w[x].x=v[x].x=f[x]=a[x]=1;
    }
    t=1;while (t<=(m<<1)) t<<=1;
    mul(t,w,v);
    for (int i=0;i<=m;i++) if (a[i]) w[i<<1].x--;
    for (int i=0;i<=m*2;i++) f[i]+=w[i].x=(int)(w[i].x/2+0.5);
    for (int i=m*2+1;i<t;i++) w[i].x=w[i].y=0;
    for (int i=0;i<=m;i++) v[i].x=a[i],v[i].y=0;
    for (int i=m+1;i<t;i++) v[i].x=v[i].y=0;
    t=1;while (t<=m*3) t<<=1;
    mul(t,w,v);
    for (int i=0;i<t;i++) u[i].x=(i&1)?0:a[i>>1];
    for (int i=0;i<=m;i++) v[i].x=a[i],v[i].y=0;
    for (int i=m+1;i<t;i++) v[i].x=0,v[i].y=0;
    mul(t,u,v);
    for (int i=0;i<=m;i++) if (a[i]) u[i*3].x--;
    for (int i=0;i<=m*3;i++) f[i]+=(long long)((w[i].x-u[i].x)/3+0.5);
    for (int i=0;i<=m*3;i++)
    if (f[i]) printf(LL,i,f[i]);
    return 0;
}