BZOJ 3771: Triple(FFT+容斥)

题面

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为方案数。
Sample Input
4

4

5

6

7

Sample Output
4 1

5 1

6 1

7 1

9 1

10 1

11 2

12 1

13 1

15 1

16 1

17 1

18 1

样例解释

11有两种方案是4+7和5+6，其他损失值都有唯一方案，例如4=4,5=5,10=4+6,18=5+6+7.
HINT
所有数据满足：Ai<=40000

解题思路

　　(FFT)，还是比较套路的那种，但是得加一个容斥。设(f[i])表示(i)这个权值出现的次数，(g[i])表示每个物品(*2)的权值出现次数，(h[i])表示每个物品(*3)的权值出现次数。那么(f)卷(f)就是两个物品一共构成的权值，再卷一个(f)就是三个物品构成的权值，但这些多算了一部分，用容斥减去。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm> 
#include<cmath>
#define int long long
 
using namespace std;
const int MAXN = 40005*3;
const double Pi = acos(-1);
 
inline int rd(){
    int x=0,f=1;char ch=getchar();
    while(!isdigit(ch)) f=ch=='-'?0:1,ch=getchar();
    while(isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
    return f?x:-x;
}
 
inline int max(int x,int y){
    return x>y?x:y;
}
 
struct Complex{
    double x,y;
    Complex(double xx=0,double yy=0){
        x=xx;y=yy;
    }
}f[MAXN<<1],g[MAXN<<1],h[MAXN<<1];
 
Complex operator+(Complex A,Complex B){return Complex(A.x+B.x,A.y+B.y);}
Complex operator-(Complex A,Complex B){return Complex(A.x-B.x,A.y-B.y);}
Complex operator*(Complex A,Complex B){return Complex(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);}
 
int n,a[MAXN],cnt[MAXN],rev[MAXN<<1],limit=1,mx,ans[MAXN],cnt2[MAXN],cnt3[MAXN];
 
void FFT(Complex *f,int type){
    for(int i=0;i<limit;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
    Complex Wn,w,tmp;int len;
    for(int i=2;i<=limit;i<<=1){
        len=i>>1;Wn=Complex(cos(Pi/len),type*sin(Pi/len));
        for(int j=0;j<limit;j+=i){
            w=Complex(1,0);
            for(int k=j;k<j+len;k++){
                tmp=w*f[k+len];f[k+len]=f[k]-tmp;
                f[k]=f[k]+tmp;w=w*Wn;
            }
        }
    }
}
 
signed main(){
    n=rd();
    for(int i=1;i<=n;i++){
        a[i]=rd();cnt[a[i]]++;cnt2[a[i]*2]++;cnt3[a[i]*3]++;
        mx=max(mx,a[i]);ans[a[i]]++;
    }
    for(int i=1;i<=mx;i++) f[i].x=cnt[i];mx*=3;
    while(limit<=mx*2) limit<<=1;
    for(int i=0;i<limit;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?(limit>>1):0);
    FFT(f,1);
    for(int i=0;i<limit;i++) g[i]=f[i],f[i]=f[i]*f[i];
    FFT(f,-1);
    for(int i=1;i<=mx;i++) ans[i]+=((int)(f[i].x/limit+0.5)-cnt2[i])/2;
    for(int i=1;i<=mx;i++) f[i].x/=limit;
    FFT(f,1);
    for(int i=0;i<limit;i++) f[i]=f[i]*g[i];
    FFT(f,-1);
    for(int i=0;i<limit;i++) h[i].x=cnt2[i];
    FFT(h,1);
    for(int i=0;i<limit;i++) h[i]=h[i]*g[i];
    FFT(h,-1);
    for(int i=1;i<=mx;i++)
        ans[i]+=((int)(f[i].x/limit+0.5)-3*(int)(h[i].x/limit+0.5)+2*cnt3[i])/6;
    for(int i=1;i<=mx;i++) 
        if(ans[i]) printf("%lld %lld
",i,ans[i]); 
    return 0;
}