[BZOJ3771] Triple

[BZOJ3771] 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为方案数。

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

试题分析

非常明显的生成函数。
设多项式：$$A(x)=sum_{i=1}^n x^{a_i}$$

[B(x)=sum_{i=1}^n x^{2a_i} ]

[C(x)=sum_{i=1}^n x^{3a_i} ]

以下先讨论带顺序的：
所以首先选一个的情况就是(A(x))。
选两个是(A(x)^2)，而且需要减去((a_i,a_i))这样的情况，也就是(B(x))，那么选两个：(A^2(x)-B(x))。
选三个需要减去(a_i,a_i,a_j)这样的，又需要加上((a_i,a_i,a_i))的情况，这样的情况被多减了两次，(A^3(x)-3A(x)B(x)+2C(x))。
所以最终答案为：$$frac{A^{3(x)-3A(x)B(x)+2C(x)}{6}+frac{A}2(x)-B(x)}{2}+A(x)$$
FFT计算即可。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
 
using namespace std;
#define LL long long
#define Pi 3.1415926535
 
inline int read(){
    int x=0,f=1; char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    return x*f;
}
const int MAXN = 1000010;
const int INF = 2147483600;
 
int l,r[MAXN+1],lim=1;
struct cpx{
    double a,b;
    cpx (double aa=0,double bb=0){a=aa; b=bb;}
}a[MAXN+1],b[MAXN+1],c[MAXN+1];
int rev[MAXN+1];
cpx operator + (cpx a,cpx b){return cpx(a.a+b.a , a.b+b.b);}
cpx operator - (cpx a,cpx b){return cpx(a.a-b.a , a.b-b.b);}
cpx operator / (cpx a,double b){return cpx(a.a / b , a.b / b);}
cpx operator * (cpx a,double b){return cpx(a.a * b , a.b * b);}
cpx operator * (cpx a,cpx b){return cpx(a.a*b.a - a.b*b.b , a.a*b.b + a.b*b.a);}
 
inline void FFT(cpx *A,int type){
    for(int i=0;i<lim;i++) if(rev[i]>i) swap(A[rev[i]],A[i]);
    for(int mid=1;mid<lim;mid<<=1){
        cpx Wn(1.0*cos(Pi/mid) , 1.0*type*sin(Pi/mid));
        for(int R=mid<<1,j=0;j<lim;j+=R){
            cpx w(1,0);
            for(int k=0;k<mid;k++,w=w*Wn){
                cpx x = A[k+j] , y = w*A[k+j+mid];
                A[k+j]=x+y; A[k+j+mid]=x-y;
            }
        }
    } return ;
}
int N,k[MAXN+1],Mx=-INF;
 
int main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    N=read();
    for(int i=1;i<=N;i++) {
        k[i]=read(),a[k[i]].a+=1;
        b[2*k[i]].a+=1; c[3*k[i]].a+=1;
        Mx=max(Mx,3*k[i]);
    } 
    while(lim<=Mx*3) lim<<=1,++l;
    for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1); FFT(b,1); FFT(c,1);
    for(int i=0;i<lim;i++) a[i]=(a[i]*a[i]*a[i]-3.0*a[i]*b[i]+2.0*c[i])/6.0+(a[i]*a[i]-b[i])/2.0+a[i];
    FFT(a,-1);
    for(int i=1;i<lim;i++){
        LL num=(LL)((a[i].a/lim)+0.5); // cout<<a[i].a<<" "<<lim<<"true"<<endl;
        if(num!=0) printf("%d %lld
",i,num);
    }
    return 0;
}