[BZOJ3771]Triple-NTT-生成函数-容斥原理

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

HINT

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

---

哇推了半天也没能意识到需要构造三个生成函数……

(╯‵□′)╯︵┻━┻

---

思路:
可以发现这种二元组或三元组的和可以被表示成卷积的形式。
同时对于这样的计数问题，通常会想到去使用生成函数。

那么考虑构造生成函数。
可以很轻松地构造出一个生成函数$a$，每个系数表示取一件物品，价值为对应指数的方案数。
那么有一个不考虑取到相同斧子和取斧子顺序问题的答案:$a+a^2+a^3$。

然而现实是有重复计算，于是需要考虑如何排除。

先考虑对于取两把斧子情况下，即$a^2$内的重复。
发现重复的情况为形如$(a,a)$的二元组。
考虑构造一个表示取两个相同斧子的函数$b$，排除顺序相同的方案，那么二元组情况下的方案数便是$\frac{a^2-b}{2}$。

然后是对于取三把斧子，即$a^3$内的重复。
首先可以发现形如$(a,b,b)$类的重复。
这样的重复可以被表示为$3ab$，其中系数$3$由于$a$的所在位置的方案数而存在。

然后考虑形如$(a,a,a)$类的重复，发现这样的重复在上面去除$(a,b,b)$时会被去除掉$3$次。
于是加回来两次，构造一个表示取三个相同斧子的函数$c$，减去$2c$。
考虑去除相同顺序的方案，那么最后三元组情况现的方案数是$\frac{a^3-3ab+2c}{6}$

于是使用NTT进行快速卷积，计算上式即可~

注意NTT的998244353模数会爆炸，可以考虑使用FFT(然而貌似不手写虚数类会很慢).
或者可以像本蒟蒻一样使用更大的且不会爆炸的模数2281701377，它的原根为3。

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N=40009;
const int M=N*4;
const ll md=2281701377ll;

ll n,m,l,v[N];
ll a[M],b[M],c[M],mv;
ll ta[M],tb[M],ans[M];
int rev[M];

inline ll read()
{
    ll x=0;char ch=getchar();
    while(ch<'0' || '9'<ch)ch=getchar();
    while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();
    return x;
}

inline ll qpow(ll a,ll b)
{
    ll ret=1;
    while(b)
    {
        if(b&1ll)
            ret=ret*a%md;
        a=a*a%md;
        b>>=1;
    }
    return ret;
}

inline void NTT(ll *a,int n,int flag)
{
    for(int i=0;i<n;i++)
        if(rev[i]<i)
            swap(a[rev[i]],a[i]);
    for(int h=2;h<=n;h<<=1)
    {
        ll wn=qpow(3,(md-1)/h);
        if(flag)wn=qpow(wn,md-2);
        for(int j=0;j<n;j+=h)
        {
            ll w=1ll;
            for(int k=j;k<j+(h>>1);k++)
            {
                ll x=a[k],y=w*a[k+(h>>1)]%md;
                a[k]=(x+y)%md;
                a[k+(h>>1)]=(x-y+md)%md;
                w=w*wn%md;
            }
        }
    }
    if(flag)
        for(ll i=0,inv=qpow(n,md-2);i<n;i++)
            a[i]=a[i]*inv%md;
}

int main()
{
    n=read();
    for(int i=1;i<=n;i++)
    {
        mv=max(mv,v[i]=read());
        a[v[i]]++;
        b[v[i]*2]++;
        c[v[i]*3]++;
    }

    for(m=1,l=0;m<=((mv*3));m<<=1)l++;
    for(int i=0;i<m;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));

    NTT(a,m,0);NTT(b,m,0);NTT(c,m,0);
    ll inv6=qpow(6ll,md-2);ll inv2=qpow(2ll,md-2);
    for(int i=0;i<m;i++)
    {
        ans[i]=((a[i]*a[i]%md*a[i]%md-3*a[i]%md*b[i]%md+2*c[i]%md)%md*inv6%md+md)%md;
        ans[i]+=((a[i]*a[i]%md-b[i])%md+md)%md*inv2%md;
        (ans[i]+=a[i])%=md;
    }
    NTT(ans,m,1);

    for(int i=0;i<=mv*3;i++)
        if(ans[i])
            printf("%d %lld\n",i,ans[i]);

    return 0;
}