【bzoj3771】Triple 【FFT 生成函数 容斥原理】

Triple

Description

我们讲一个悲伤的故事。
从前有一个贫穷的樵夫在河边砍柴。
这时候河里出现了一个水神，夺过了他的斧头，说：
“这把斧头，是不是你的？”
樵夫一看：“是啊是啊！”
水神把斧头扔在一边，又拿起一个东西问：
“这把斧头，是不是你的？”
樵夫看不清楚，但又怕真的是自己的斧头，只好又答：“是啊是啊！”
水神又把手上的东西扔在一边，拿起第三个东西问：
“这把斧头，是不是你的？”
樵夫还是看不清楚，但是他觉得再这样下去他就没法砍柴了。
于是他又一次答：“是啊是啊！真的是！”
水神看着他，哈哈大笑道：
“你看看你现在的样子，真是丑陋！”
之后就消失了。
樵夫觉得很坑爹，他今天不仅没有砍到柴，还丢了一把斧头给那个水神。
于是他准备回家换一把斧头。
回家之后他才发现真正坑爹的事情才刚开始。
水神拿着的的确是他的斧头。
但是不一定是他拿出去的那把，还有可能是水神不知道怎么偷偷从他家里拿走的。
换句话说，水神可能拿走了他的一把，两把或者三把斧头。
樵夫觉得今天真是倒霉透了，但不管怎么样日子还得过。
他想统计他的损失。
樵夫的每一把斧头都有一个价值，不同斧头的价值不同。总损失就是丢掉的斧头价值和。
他想对于每个可能的总损失，计算有几种可能的方案。
注意：如果水神拿走了两把斧头a和b，(a,b)和(b,a)视为一种方案。拿走三把斧头时，(a,b,c),(b,c,a),(c,a,b),(c,b,a),(b,a,c),(a,c,b)视为一种方案。

Input

第一行是整数N，表示有N把斧头。
接下来n行升序输入N个数字Ai，表示每把斧头的价值。

Output

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

HINT

所有数据满足：Ai<=40000
题解：我们搞3个多项式A，B，C，分别代表每个物品选一个、选两个、选三个。
（这玩意好像叫做生成函数）
设当前要求价值为i的答案数。
选一把斧头的情况：答案为A[i]。
选两把斧头的情况：答案为(A[i]∗A[i]−B[i])/2$(A[i]\ast A[i]-B[i])/2$。A[i]*A[i]中包含了选相同的，所以要减去B[i]，还要除掉排列数。
选三把斧头的情况：答案为(A[i]∗A[i]∗A[i]−3∗A[i]∗B[i]+2∗C[i])/6$(A[i]\ast A[i]\ast A[i]-3\ast A[i]\ast B[i]+2\ast C[i])/6$。具体就是A[i]∗A[i]∗A[i]$A[i]\ast A[i]\ast A[i]$包含了一种选三个一样的和三种选两个 一样的，A*B包含了一种选三个一样的和一种选两个一样的。(A[i]∗A[i]∗A[i]−3∗A[i]∗B[i])$(A[i]\ast A[i]\ast A[i]-3\ast A[i]\ast B[i])$多减了两种选三个一样的，所以再加上两个C[i]。最后除掉排列数。手动容斥理解一下就好了。
具体见代码。

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=270005;
const double pi=3.141592653589793;
int n,m,x,rev[N];
long long ans;
struct complex{
    double x,y;
    complex(){
        x=y=0;
    }
    complex(double x,double y):x(x),y(y){}
    friend complex operator + (const complex &a,const complex &b){
        return complex(a.x+b.x,a.y+b.y);
    }
    friend complex operator - (const complex &a,const complex &b){
        return complex(a.x-b.x,a.y-b.y);
    }
    friend complex operator * (const complex &a,const complex &b){
        return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
    }
    friend complex operator * (const complex &a,const double &b){
        return complex(a.x*b,a.y*b);
    }
    friend complex operator / (const complex &a,const double &b){
        return complex(a.x/b,a.y/b);
    }
}a[N],b[N],c[N];
void fft(complex *a,int dft){
    for(int i=0;i<n;i++){
        if(i<rev[i]){
            swap(a[i],a[rev[i]]);
        }
    }
    for(int i=1;i<n;i<<=1){
        complex wn=complex(cos(pi/i),dft*sin(pi/i));
        for(int j=0;j<n;j+=i<<1){
            complex w=complex(1,0),x,y;
            for(int k=j;k<j+i;k++,w=w*wn){
                x=a[k];
                y=w*a[k+i];
                a[k]=x+y;
                a[k+i]=x-y;
            }
        }
    }
    if(dft==-1){
        for(int i=0;i<n;i++){
            a[i]=a[i]/n;
        }
    }
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        scanf("%d",&x);
        a[x].x++;
        b[2*x].x++;
        c[3*x].x++;
        m=max(m,3*x);
    }
    for(n=1;n<=m;n<<=1);
    for(int i=0;i<n;i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
    }
    fft(a,1);
    fft(b,1);
    fft(c,1);
    for(int i=0;i<n;i++){
        a[i]=a[i]+(a[i]*a[i]-b[i])/2+(a[i]*a[i]*a[i]-a[i]*b[i]*3+c[i]*2)/6;
    }
    fft(a,-1);
    for(int i=0;i<n;i++){
        ans=(long long)(a[i].x+0.5);
        if(ans){
            printf("%d %lld\n",i,ans);
        }
    }
    return 0;
}