题解：P9896 [ICPC 2018 Qingdao R] Sub-cycle Graph

· · 题解

非常好组合数学，使我推可爱的大柿子，而不是什么看起来就令人畏惧的生成函数。

抢到最优解了，于是考虑写个题解。

考虑一个“半环图”是什么样子的，不难发现其实就是一堆链和单点。于是发现我们可以按照度数把点分类，即度数为 0,1,2 的三类点。注意到度数为 0 的点看起来就很特殊，我们考虑暴力枚举它的个数。

设枚举的 0 度数点个数为 k，当前求解的图点数为 n，边数为 m。我们先假定当前 0<m<n 成立，一步一步推导整个式子。

首先，当然是要从 n 个点中挑选出 k 个度数为 0 的点，答案是：\begin{pmatrix} n\\ k\\ \end{pmatrix}。

其次，我们考虑剩余的点。不难发现所有点一定构成 n-m 个连通块，而其中有 k 个单点，剩下的链数则为 n-m-k，每条链有两个端点，则度数为 1 的点数为 2 \times (n-m-k)，度数为 2 的点个数为 2 \times m+k-n，于是我们要在剩余的点中挑选出度数为 1 的点，答案是 \begin{pmatrix} n-k\\ 2 \times (n-m-k)\\ \end{pmatrix}。

接下来，考虑将这些度数为 1 的点两两配对形成链的两端。我们想到规定每条链的一个端点，乘上另一组点的排列，每一个排列都对应着一个合法的方案。但是这样显然会算重。不难发现，对于一种方案中连有边的两个点，将其所在点集互换，仍然为一个方案，这就是算重的部分。因此，最终的答案需要将其去掉。则答案为 \begin{pmatrix} 2 \times (n-m-k)\\ n-m-k\\ \end{pmatrix}(n-m-k)! \div 2^{(n-m-k)}

最后，我们需要将剩余度数为 2 的点插入当前的链中。同样的，即使插入同一个链，其顺序也会有影响，因此考虑先确定剩余点的一个排列。剩下的就是一个球无标号盒有标号的经典问题，插板法可得这一步的答案是 \begin{pmatrix} m-1\\ n-m-k-1\\ \end{pmatrix} \times (2 \times m+k-n)！

于是根据乘法原理相乘即可得到结果为：

\begin{pmatrix} n\\ k\\ \end{pmatrix} \begin{pmatrix} n-k\\ 2 \times (n-m-k)\\ \end{pmatrix} \begin{pmatrix} 2 \times (n-m-k)\\ n-m-k\\ \end{pmatrix}(n-m-k)! \begin{pmatrix} m-1\\ n-m-k-1\\ \end{pmatrix} (2 \times m+k-n)！\div 2^{(n-m-k)}

其中，k 的取值范围可以根据 n 和 m 求得，对所有 k 的可能取值求解上式并求和即可。预处理后上式可以 O(1) 求出。

但是，我们其实可以注意到，上式如果化简，则可以消去大量冗余的阶乘。最终化简后的形式实际上如下：

\frac{(n)!(m-1)!}{(k)!(n-m-k)!(2 \times m+k-n)!(n-m-k-1)!2^{(n-m-1)}}

发现好算了不少，于是可以求得飞快。

给出代码：

#include<bits/stdc++.h>
#define R register
using namespace std;
typedef long long ll;
const int N=1e5+20,mod=1e9+7;
int jc[N],inv[N],inv2[N];
inline int pwr(int a,int b){
    int ans=1;
    for(;b;b>>=1){
        if(b&1)ans=(ll)ans*a%mod;
        a=(ll)a*a%mod;
    }
    return ans;
}
inline int calc(int i,int j,int k){
    int ans=1;
    ans=1ll*ans*jc[j-1]%mod;
    ans=1ll*ans*inv[k]%mod;
    ans=1ll*ans*inv[i-j-k]%mod;
    ans=1ll*ans*inv2[i-j-k]%mod;
    ans=1ll*ans*inv[2*j+k-i]%mod;
    ans=1ll*ans*inv[i-j-k-1]%mod;
    return ans;
}
int main()
{
    int T;
    scanf("%d",&T);
    int mx=N-10;
    jc[0]=1;
    for(R int i=1;i<=mx;i++)jc[i]=(ll)i*jc[i-1]%mod;
    inv[mx]=pwr(jc[mx],mod-2);
    inv2[mx]=pwr(pwr(2,mx),mod-2);
    for(R int i=mx-1;i>=0;i--){
        inv[i]=(ll)inv[i+1]*(i+1)%mod;
        inv2[i]=(ll)inv2[i+1]*2ll%mod;
    }
    while(T--){
        int x,y;
        scanf("%d%d",&x,&y);
        if(x==y){
            printf("%d\n",(ll)jc[x-1]*inv2[1]%mod);
            continue;
        }
        if(y>x){
            printf("%d\n",0);
            continue;
        }
        if(y==0){
            printf("%d\n",1);
            continue;
        } 
        int ans=0;
        for(int z=max(0,x-2*y);z<=x-y-1;z++){
            ans=(ans+calc(x,y,z))%mod;
        }
        ans=(ll)ans*jc[x]%mod;
        printf("%d\n",ans);
    }
    return 0;
}