51Nod1252: 桥与海港

推荐游戏:Poly Bridge

Steam 上的 Poly Bridge 太贵了买不起。

传送门

题目大意

有 $n$ 个带编号的岛屿,把这 $n$ 个岛屿分成 $m$ 个区域,每个岛屿属于这 $m$ 个区域之一,每个区域至少包含一个岛屿。

现在在岛屿之间建造桥梁,要求桥两端的岛屿必须属于同一区域。所以,恰有 $n-m$ 座桥来维护每个区域内岛屿的连通性。

至于区域之间,计划在每个区域选一个岛屿建造一个海港,要求这个岛屿不能直接和多于 $k$ 座桥相连。求方案数模 $x$ 的值。

多组数据,其中数据组数 $\le 10$,$2\le n\le 10^9,1\le m,k\le \min(n,100),1< x < 2^{31}$。

说实话我觉得 Bridge Constructor Portal 更好玩,虽然没有船开但是可以整顿公司搭顺风车的行为

题解

我们先搞清楚这道题到底是要我们干嘛。

题目真的叙述得太长了,不想认真看。【我一开始甚至看成每个岛屿都只能与 $k$ 个桥相连。

我用 OI 的语言说一遍题面吧:

带标号的 $n$ 个点组成 $m$ 棵有根树,每棵树的根节点的度数不超过 $k$。求方案数模 $x$ 的值。

这样就简洁了。

就像是一道标准的图论 + 组合数学一样。

但是好难啊,不会,怎么办呢?

还是一样的,先简化问题。


prufer 编码

我们先考虑 $n$ 个点形成 $m$ 个有根树怎么做。

还是不会。。。

那 $n$ 个点,一棵有根树呢?

prufer 编码。

如果你还不知道 prufer 编码,这里推荐一篇博客:树的计数 + prufer序列与Cayley公式 学习笔记

好了,现在你一定知道 prufer 编码了。我们不关心它的证明,只需要知道它的工作方式就行了。

那么基础的基础,$n$ 个点一棵无根树的方案数为 $n^{n-2}$。

有根树呢?

随便选一个点作为根,有 $n$ 种方案,共 $n^{n-2}\times n=n^{n-1}$。

森林呢?

。。。

不会。

考虑森林就是一棵树删去一些边,直接乘上 $\mathrm C_{n-1}^{m-1}$?

显然是错的。

因为对同一个森林,变成树不止一种方法,并没有一一对应。

所以这里留一个坑:带标号无根森林。

我只会有根森林的计数。

这里我们就要换一种方式来理解有根树的方案数了。

我们看一看公式 $n^{n-1}$,就像是 prufer 编码多了一位一样,对吧?

那么我们就把这一位在编码中补出来。

首先是对叶子结点的定义改为:非根的度数为 $1$ 的节点。

然后在编码时,将结束条件改为只剩下根节点。

这样,我们就得到了有根树对应的一个编码。

如何还原呢?

同无根树的操作,只是最后点集里会剩下一个点,这个点就是根节点。

那么有根森林能不能这样搞呢?

不能。自己试一试就知道了,还原的时候可能把别的树的根节点给连到另一棵树上去。

所以我们需要变化一下。

我们已经知道了,$set$ 中剩下的一定是每棵树的根,那为了防止它被连走,为什么不先把它们拿出来呢?

这样子再去试,诶,好像可以!

所以步骤变成了:先选 $m$ 个点作为根,再来编码。

所以有根森林的方案数为 $\mathrm C_{n}^{m}\times n^{n-m}$。

然后去写发现是错的。。。

emmm.

好吧可能我的确很菜。【其实一样就该看出来不对了,比如 $m=1$。

那怎么搞,这会是真的不会了啊!

适当打表发现,诶?好像是 $\mathrm C_{n-1}^{m-1} \times n^{n-m}$。

其实可以发现编码的最后一位一定是根,没得选的,一定得拿走,只能从剩下的 $n-1$ 个点选 $m-1$ 个根,所以是 $\mathrm C_{n-1}^{m-1}$。


神奇 DP

好了,我们现在已经会 $n$ 个点 $m$ 棵树的有根森林了,是时候来解决根的度数不大于 $k$ 了。

妈的还是不会。。。

要是能只管根和第一层的节点就好了,这样看着都要舒服,没有下面那堆糟心的玩意儿。

那就来看看把根节点取走图变成了什么吧。

现在这个图有 $n-m$ 个节点,$\sum_{i=1}^{m}{son_i}$ 棵树,树根为每个根节点的儿子的森林。好像,没限制?

设 $p=\sum_{i=1}^{m}{son_i}$,方案数为 $\mathrm C_{n-m-1}^{p-1}\times (n-m)^{n-m-p}$。

诶,真舒爽。

然后就是把这些树根插回那些根节点上去了。

这个子问题分得真是精髓,两个问题看上去似乎都能解决诶。

设 $dp_{i,j}$ 表示前 $i$ 个根节点连上了 $j$ 个树根的方案数。

$$dp_{i,j}=\sum_{x=0}^{k}{dp_{i-1,j-x}\times\mathrm C_{j}^{k}}$$

答案为 $dp_{m,p}$。

但是好像时间复杂度。。。

$\mathrm O(m\times m\times k\times k=10^8)$。

无法接受。

那只能优化 DP 方程了。

考虑第 $j$ 个节点,它可以连到前 $i$ 个根中的任意一个。什么时候不合法呢?就是这个树根还连了 $k$ 个节点的时候。

那把那部分减掉不就好了吗?某一个根节点连了 $k$ 个,不就是其它 $i-1$ 个根节点连了 $j-k-1$ 个吗。

$$dp_{i,j}=i\times(dp_{i,j-1}-dp_{i-1,j-k-1}\times\mathrm C_{j-1}^{k})$$

时间复杂度 $\mathrm O(m\times m\times k=10^6)$,ok 了。


大数组合

最后的最后,我们只需要把答案乘上一个 $\mathrm C_{n}^{m}$,这道题就算完结啦!

n <= 1e9

emmm.

Lucas?

x < 2^31

花Q!

花Q,花Q!

花Q!

fnndp.

啊再见再见,这道题不适合我这种不会 CRT 的人。。。

贴一个 Lucas + CRT 的博客:中国剩余定理与扩展 Lucas定理与扩展 学习笔记

像我这么蠢的人自然是不会这种高精难的玩意的,我们需要更亲民的方法。

首先我们可以发现,这个组合数虽然难算,但是它的项数并不多,就是取模有点麻烦。

为什么取模会影响到计算呢?因为枚举到的分子可能包含模数,就变成 $0$ 了,但是之后的分母又有可能把它的一些因数消掉,导致答案并不是 $0$。

emmm.

那我们可以把模数的因数单独算啊,剩下的反正也没影响,直接 exgcd 把逆元乘进去就好了啊。

又到了一题解一度的第三节出代码环节。

//C(n,m)%mod
//pm[i],ct:模数的质因数及个数
//pc[i]:质因数出现次数
//inv():逆元
//fp():快速幂
for(int i=1;i<=m;i++){
    int a=n+1-i,b=i;
    for(int j=0;j<ct;j++){
        while(a%pm[j]==0)pc[j]++,a/=pm[j];
        while(b%pm[j]==0)pc[j]--,b/=pm[j];
    }
    (ans*=a*inv(b))%=mod;
}
for(int i=0;i<ct;i++)
    if(pc[i])
        (ans*=fp(pm[i],pc[i]))%=mod;

诶你看这不就解决了吗。时间复杂度 $\mathrm O(m\times \log{x}\times \log{n})$,随便跑。

Lucas + CRT 我不知道时间复杂度是多少,应该比这个快,估计是 $m\rightarrow \log{m}$。【看着就贼复杂,谁想写啊。


顺便贴一张图。


代码

/*"-Aria on the Planets-"*/
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double D;
typedef pair<int,int> Pii;
#define mp(a,b) make_pair(a,b)
#define fir first
#define sec second
inline int min(int _a,int _b){return _a<_b?_a:_b;}
inline int max(int _a,int _b){return _a>_b?_a:_b;}
template<class _T>inline void rd(_T &_a){int _f=0,_ch=getchar();_a=0;while(_ch<'0'||_ch>'9'){if(_ch=='-')_f=1;_ch=getchar();}while(_ch>='0'&&_ch<='9')_a=(_a<<1)+(_a<<3)+_ch-'0',_ch=getchar();if(_f)_a=-_a;}
const int inf=0x3f3f3f3f;
const D eps=1e-8;
const int N=46342,M=105;
int mod,ct,isp[N],p[N],pm[10],pc[10],pw[M*M],dp[M][M*M];LL c[M*M];
void eg(int a,int b,int &x,int &y){b?(eg(b,a%b,y,x),y-=a/b*x):(x=1,y=0);}
int inv(int a){int x,y;eg(a,mod,x,y);return x;}
int fp(LL a,int b){LL r=1;for(;b;b>>=1,(a*=a)%=mod)if(b&1)(r*=a)%=mod;return r;}
void ud(int a,int b,LL &tmp){for(int i=0;i<ct;i++){while(a%pm[i]==0)pc[i]++,a/=pm[i];while(b%pm[i]==0)pc[i]--,b/=pm[i];}(tmp*=1ll*a*inv(b)%mod)%=mod;}
void ud(int a,int b,LL &tmp,LL &x){ud(a,b,tmp);x=tmp;for(int i=0;i<ct;i++)if(pc[i])(x*=fp(pm[i],pc[i]))%=mod;}
int main(){
    //clock_t start=clock();
    //freopen("test.in","r",stdin);
    //freopen("test.out","w",stdout);
    int cnt=0;for(int i=2;i<N;i++){if(!isp[i])p[cnt++]=i;for(int j=0;j<cnt&&i*p[j]<N;j++){isp[i*p[j]]=1;if(i%p[j])break;}}
    int T;rd(T);while(T--){
        int n,m,q,mx;LL ans=0;rd(n);rd(m);rd(q);rd(mod);if(n==m){printf("%d\n",1%mod);continue;}mx=min(n-m,m*q);
        LL tmp=mod;ct=0;for(int i=0;i<cnt&&p[i]*p[i]<=tmp;i++)if(tmp%p[i]==0)for(pm[ct++]=p[i];tmp%p[i]==0;tmp/=p[i]);if(tmp>1)pm[ct++]=tmp;
        memset(c,0,q*sizeof(int));memset(pc,0,sizeof pc);c[q]=1;tmp=1;for(int i=q+1;i<mx;i++)ud(i,i-q,tmp,c[i]);
        dp[0][0]=1;for(int i=1;i<=m;i++){dp[i][0]=1;for(int j=1;j<=mx;j++)dp[i][j]=1ll*i*(dp[i][j-1]-(j>q?c[j-1]*dp[i-1][j-1-q]%mod:0))%mod;}
        pw[mx]=fp(n-m,n-m-mx);for(int i=mx-1;~i;i--)pw[i]=1ll*pw[i+1]*(n-m)%mod;
        memset(pc,0,sizeof pc);c[0]=1;for(int i=tmp=1;i<mx;i++)ud(n-m-i,i,tmp,c[i]);
        for(int i=1;i<=mx;i++)(ans+=1ll*c[i-1]*dp[m][i]%mod*pw[i])%=mod;
        memset(pc,0,sizeof pc);for(int i=tmp=1;i<=m;i++)ud(n+1-i,i,tmp);(ans*=tmp)%=mod;for(int i=0;i<ct;i++)if(pc[i])(ans*=fp(pm[i],pc[i]))%=mod;
        printf("%lld\n",(ans+mod)%mod);
    }
    //printf("\n%dms",(int)((D)(clock()-start)/CLOCKS_PER_SEC*1000));
    return 0;
}