1. 程式人生 > >bzoj 3771: Triple【生成函式+FFT+容斥原理】

bzoj 3771: Triple【生成函式+FFT+容斥原理】

瞎搞居然1A,真是吃鯨
欸今天寫不完了orz

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=500005;
int n,m,lm,bt,ans[N],re[N],x[N],q[N];
struct cd
{
    double a,b;
    cd(double A=0,double B=0)
    {
        a=A,b=B;
    }
    cd operator + (const cd &x) const
    {
        return cd(a+x.a,b+x.b);
    }
    cd operator - (const cd &x) const
    {
        return cd(a-x.a,b-x.b);
    }
    cd operator * (const cd &x) const
    {
        return cd(a*x.a-b*x.b,a*x.b+b*x.a);
    }
}a[N],b[N],c[N],d[N],e[N];
int read()
{
    int r=0,f=1;
    char p=getchar();
    while(p>'9'||p<'0')
    {
        if(p=='-')
            f=-1;
        p=getchar();
    }
    while(p>='0'&&p<='9')
    {
        r=r*10+p-48;
        p=getchar();
    }
    return r*f;
}
void dft(cd a[],int f)
{
    for(int i=0;i<lm;i++)
        if(i<re[i])
            swap(a[i],a[re[i]]);
    for(int i=1;i<lm;i<<=1)
    {
        cd wi=cd(cos(M_PI/i),f*sin(M_PI/i));
        for(int k=0;k<lm;k+=(i<<1))
        {
            cd w=cd(1.0,0.0),x,y;
            for(int j=0;j<i;j++)
            {
                x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y,a[j+k+i]=x-y;
                w=w*wi;
            }
        }
    }
    if(f==-1)
        for(int i=0;i<lm;i++)
            a[i].a/=lm;
}
void fft(cd a[],cd b[])
{
    dft(a,1);
    dft(b,1);
    for(int i=0;i<lm;i++)
        a[i]=a[i]*b[i];
    dft(a,-1);
}
int main()
{
    n=read();
    for(int i=1;i<=n;i++)
    {
        x[i]=read(),m=max(m,3*x[i]);
        a[x[i]].a+=1,b[x[i]].a+=1,c[x[i]].a+=1,d[x[i]*2].a+=1,e[x[i]].a+=1;
        ans[x[i]]+=1;
    }
    for(bt=0;(1<<bt)<=2*m;bt++);
    lm=(1<<bt);
    for(int i=0;i<lm;i++)
        re[i]=(re[i>>1]>>1)|((i&1)<<(bt-1));
    fft(a,b);
    for(int i=0;i<lm;i++)
        q[i]=int(a[i].a+0.5);
    for(int i=1;i<=n;i++)
        q[2*x[i]]-=1;
    for(int i=0;i<lm;i++)
        q[i]/=2,ans[i]+=q[i];
    fft(c,a);
    fft(d,e);
    for(int i=0;i<lm;i++)
        q[i]=3*(int)(d[i].a+0.5);
    for(int i=1;i<=n;i++)
        q[3*x[i]]-=2;
    for(int i=0;i<lm;i++)
        ans[i]+=((int)(c[i].a+0.5)-q[i])/6;
    for(int i=0;i<lm;i++)
        if(ans[i])
            printf("%d %d\n",i,ans[i]);
    return 0;
}