BZOJ 3771: Triple [快速傅里叶变换 生成函数 容斥原理]

#include <iostream>
#include <cstdio>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;
const int N=131072+5;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1;c=getchar();}
    while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘;c=getchar();}
    return x*f;
}
const double PI=acos(-1);
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
};
typedef Vector CD;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,Vector b){return Vector(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
Vector operator *(Vector a,double b){return Vector(a.x*b,a.y*b);}
Vector operator /(Vector a,double b){return Vector(a.x/b,a.y/b);}
Vector conj(Vector a){return Vector(a.x,-a.y);}

struct FastFourierTransform{
    int n,rev[N];
    CD omega[N],omegaInv[N];
    void ini(int m){
        n=1;
        while(n<m) n<<=1; 
        for(int k=0;k<n;k++) 
            omega[k]=CD(cos(2*PI/n*k),sin(2*PI/n*k)),
            omegaInv[k]=conj(omega[k]);
        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }
    }
    void transform(CD *a,CD *omega){
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            for(CD *p=a;p!=a+n;p+=l)
                for(int k=0;k<m;k++){
                    CD t=omega[n/l*k]*p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                }
        }
    }
    void DFT(CD *a,int flag){
        if(flag==1) transform(a,omega);
        else{
            transform(a,omegaInv);
            for(int i=0;i<n;i++) a[i].x/=(double)n;
        }
    }
}fft;
int n,m,w;
CD A[N],B[N],C[N],ans[N];
int main(){
    freopen("in","r",stdin);
    n=read();
    for(int i=1;i<=n;i++){
        w=read(); m=max(m,w);
        A[w].x=1;
        B[w*2].x=1;
        C[w*3].x=1;
    }
    m=m*3;
    fft.ini(m);
    fft.DFT(A,1);fft.DFT(B,1);fft.DFT(C,1);
    
    for(int i=0;i<fft.n;i++)
        ans[i]=ans[i]+A[i]+(A[i]*A[i]-B[i])/2.0+(A[i]*A[i]*A[i]-3.0*A[i]*B[i]+2.0*C[i])/6.0;
    fft.DFT(ans,-1);
    for(int i=0;i<m;i++) if(int(ans[i].x+0.5)) printf("%d %d\n",i,int(ans[i].x+0.5));
}