[BZOJ 3771] Triple

3771. 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
样例解释
11有两种方案是4+7和5+6，其他损失值都有唯一方案，例如4=4,5=5,10=4+6,18=5+6+7.

Hint

所有数据满足：Ai<=40000

题解

假设我们没有不能重复选取而且无序计数的限制, 直接构造多项式 $A(x)=\sum\limits_{k\in S} x^k$ 然后计算 $A(x)+A(x)^2+A(x)^3$ 就好了

然而这题有点诡异...因为它要求的东西是无序对而且还不允许选两次...这就非常蠢...

好的那么我们想一想如果直接求 $A(x)+A(x)^2+A(x)^3$ 会发生什么:

对于形如 $(a,b)$ 的值, 我们会算两遍;

对于形如 $(a,b,c)$ 的值, 我们会算 $6$ 遍;

对于形如 $(a,a,b)$ 的值, 我们会算 $3$ 遍;

对于形如 $(a,a)$ 和 $(a,a,a)$ 的值, 我们会算 $1$ 遍(然而这些值不合法)

因为情况比较少, 我们手动容斥一发.

我们设 $A_1(x)$ 为原来的 $A(x)$, $A_2(x)$ 为将其中的所有项指数扩大两倍的多项式, $A_3(x)$ 同理. 这样我们就可以把两个同样的值绑定在一起计算方案.

然后我们对不同大小的组合分别计算:

大小为 $1$ 的组合, 显然就是 $A_1(x)$;

大小为 $2$ 的组合, $(a,b)$ 算了两遍, $(a,a)$ 算了一遍, 而 $(a,a)$ 不合法, 我们把它减掉. 而 $(a,a)$ 的方案数就是 $A_2(x)$, 于是这部分的答案就是 $\frac 1 2 \left(A_1(x)^2-A_2(x)\right)$

大小为 $3$ 的组合, $(a,b,c)$ 算了 $6$ 遍, $(a,a,b)$ 算了三遍, 我们想办法计算出 $(a,a,b)$ 的个数并减掉: 我们如果直接计算 $A_1(x)A_2(x)$ 的话又会把 $(a,a,a)$ 算上一遍, 所以 $(a,a,b)$ 的个数就是 $A_1(x)A_2(x)-A_3(x)$. $(a,a,a)$ 的个数就是 $A_3(x)$, 于是我们就可以得到总的方案数是 $\frac 1 6 \left(A_1(x)^3-3\left(A_1(x)A_2(x)-A_3(x)\right)-A_3(x)\right)$

于是总答案就是长这样的:

$$ A_1(x)+\frac 1 2 \left(A_1(x)^2-A_2(x)\right)+\frac 1 6 \left(A_1(x)^3-3\left(A_1(x)A_2(x)-A_3(x)\right)-A_3(x)\right) $$

代码实现

全家桶板子真好用

#include <bits/stdc++.h>

const int G=3;
const int DFT=1;
const int IDFT=-1;
const int MAXN=550000;
const int MOD=998244353;
const int INV2=(MOD+1)>>1;
const int PHI=MOD-1;

typedef std::vector<int> Poly;

Poly Sqrt(Poly);
void Read(Poly&);
Poly Inverse(Poly);
Poly Ln(const Poly&);
Poly Exp(const Poly&);
void Print(const Poly&);
void NTT(Poly&,int,int);
Poly Pow(const Poly&,int);
Poly Integral(const Poly&);
Poly Derivative(const Poly&);
Poly operator*(int,Poly);
Poly operator*(Poly,Poly);
Poly operator/(Poly,Poly);
Poly operator%(Poly,Poly);
Poly operator+(const Poly&,const Poly&);
Poly operator-(const Poly&,const Poly&);

int w[MAXN];
int rev[MAXN];

int NTTPre(int);
int Pow(int,int,int);

int main(){
    int n;
    scanf("%d",&n);
    int maxv=0;
    for(int i=0;i<n;i++){
        scanf("%d",w+i);
        maxv=std::max(maxv,w[i]);
    }
    Poly s1(maxv+1),s2(maxv*2+1),s3(maxv*3+1);
    for(int i=0;i<n;i++){
        s1[w[i]]=1;
        s2[2*w[i]]=1;
        s3[3*w[i]]=1;
    }
    Poly ans=s1+INV2*(s1*s1-s2)+Pow(6,MOD-2,MOD)*(s1*s1*s1-3*(s1*s2-s3)-s3);
    for(int i=1;i<=maxv*3;i++)
        if(ans[i])
            printf("%d %d\n",i,ans[i]);
    return 0;
}

void Read(Poly& a){
    for(auto& i:a)
        scanf("%d",&i);
}

void Print(const Poly& a){
    for(auto i:a)
        printf("%d ",i);
    puts("");
}

Poly operator*(int k,Poly a){
    for(auto& i:a)
        i=1ll*k*i%MOD;
    return a;
}

Poly Pow(const Poly& a,int k){
    Poly log=Ln(a);
    for(auto& i:log)
        i=1ll*i*k%MOD;
    return Exp(log);
}

Poly Sqrt(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,int(sqrt(a[0])));
    else{
        Poly b=a;
        b.resize((len+1)>>1);
        b=Sqrt(b);
        b.resize(len);
        Poly inv=Inverse(b);
        int bln=NTTPre(inv.size()+a.size());
        NTT(a,bln,DFT);
        NTT(inv,bln,DFT);
        for(int i=0;i<bln;i++)
            a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
        NTT(a,bln,IDFT);
        for(int i=0;i<len;i++)
            b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
        return b;
    }
}

Poly Exp(const Poly& a){
    size_t len=1;
    Poly ans(1,1),one(1,1);
    while(len<(a.size()<<1)){
        len<<=1;
        Poly b=a;
        b.resize(len);
        ans=ans*(one-Ln(ans)+b);
        ans.resize(len);
    }
    ans.resize(a.size());
    return ans;
}

Poly Ln(const Poly& a){
    Poly ans=Integral(Derivative(a)*Inverse(a));
    ans.resize(a.size());
    return ans;
}

Poly Integral(const Poly& a){
    int len=a.size();
    Poly ans(len+1);
    for(int i=1;i<len;i++)
        ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
    return ans;
}

Poly Derivative(const Poly& a){
    int len=a.size();
    Poly ans(len-1);
    for(int i=1;i<len;i++)
        ans[i-1]=1ll*a[i]*i%MOD;
    return ans;
}

Poly operator/(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans(1);
    if(n>=m){
        std::reverse(a.begin(),a.end());
        std::reverse(b.begin(),b.end());
        b.resize(n-m+1);
        ans=Inverse(b)*a;
        ans.resize(n-m+1);
        std::reverse(ans.begin(),ans.end());
    }
    return ans;
}

Poly operator%(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans;
    if(n<m)
        ans=a;
    else
        ans=a-(a/b)*b;
    ans.resize(m);
    return ans;
}

Poly operator*(Poly a,Poly b){
    int len=a.size()+b.size()-1;
    int bln=NTTPre(len);
    NTT(a,bln,DFT);
    NTT(b,bln,DFT);
    for(int i=0;i<bln;i++)
        a[i]=1ll*a[i]*b[i]%MOD;
    NTT(a,bln,IDFT);
    a.resize(len);
    return a;
}

Poly operator+(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+b[i])%MOD;
    return ans;
}

Poly operator-(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+MOD-b[i])%MOD;
    return ans;
}

Poly Inverse(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Pow(a[0],MOD-2,MOD));
    else{
        Poly b(a);
        b.resize((len+1)>>1);
        b=Inverse(b);
        int bln=NTTPre(b.size()*2+a.size());
        NTT(a,bln,DFT);
        NTT(b,bln,DFT);
        for(int i=0;i<bln;i++)
            b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
        NTT(b,bln,IDFT);
        b.resize(len);
        return b;
    }
}

void NTT(Poly& a,int len,int opt){
    a.resize(len);
    for(int i=0;i<len;i++)
        if(rev[i]>i)
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        int step=i<<1;
        int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
        for(int j=0;j<len;j+=step){
            int w=1;
            for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
                int x=a[j+k];
                int y=1ll*a[j+k+i]*w%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(opt==IDFT){
        int inv=Pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)
            a[i]=1ll*a[i]*inv%MOD;
    }
}

int NTTPre(int n){
    int bln=1,bct=0;
    while(bln<n){
        bln<<=1;
        ++bct;
    }
    for(int i=0;i<bln;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
    return bln;
}

inline int Pow(int a,int n,int p){
    int ans=1;
    while(n>0){
        if(n&1)
            ans=1ll*a*ans%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}