随机函数
题解
首先，我们考虑怎么算这个概率。
显然一种方法是对于每一对
(
f
a
,
f
b
)
(f_a,f_b)
(fa?,fb?)计算在它们中选择一对
(
x
,
y
)
(x,y)
(x,y)使得
f
a
(
x
)
=
f
b
(
y
)
f_a(x)=f_b(y)
fa?(x)=fb?(y)的概率，记作
p
a
,
b
p_{a,b}
pa,b?。
由于每对数之间是独立的，所以最后成立的概率为
p
a
,
b
n
p_{a,b}^n
pa,bn?。
我们可以每次取枚举映射到
i
i
i的数在两边有多少个，求出它的概率，在变成其的
n
n
n次方。
这样的做法是
O
(
n
7
)
O\left(n^7\right)
O(n7)的，不太好优化，因为它只能单个求出来计算
n
n
n次方，不好合并。
当然，你可以通过枚举整数拆分，做到拆分数平方级别的复杂度，也可以过。
我们看我们上面的做法是通过枚举函数
(
a
,
b
)
(a,b)
(a,b)，计算所有二分图的概率和，我们得考虑别的方面，比如说枚举二分图，计算合法函数的概率。
显然，每对
(
x
i
,
y
i
)
(x_i,y_i)
(xi?,yi?)就相当于是二分图上的一条边，二分图上有
n
n
n条边，允许有重边。
我们记
d
p
x
,
y
,
i
,
c
dp_{x,y,i,c}
dpx,y,i,c?表示左边用了
x
x
x个点右边用了
y
y
y个点，总共用了
i
i
i条不同的边，产生了
c
c
c个连通块的方案数。
显然，我们可以比较轻松的计算出答案为:
A
n
s
=
∑
x
,
y
,
i
,
c
S
(
n
,
i
)
i
!
(
m
x
)
(
m
y
)
d
p
x
,
y
,
i
,
c
m
?
2
m
+
x
+
y
?
c
Ans=\sum_{x,y,i,c}S(n,i)i!\binom{m}{x}\binom{m}{y}dp_{x,y,i,c}m^{-2m+x+y-c}
Ans=x,y,i,c∑?S(n,i)i!(xm?)(ym?)dpx,y,i,c?m?2m+x+y?c但我们又要怎么快速地算出我们的
d
p
dp
dp值呢？
我们可以通过类似背包的形式来计算。
我们先记
t
x
,
y
,
e
t_{x,y,e}
tx,y,e?表示左边
x
x
x个点，右边
y
y
y个点，用了
e
e
e条边，形成不带孤点的有标号二分图方案数，显然，这可以通过容斥得到，枚举孤点的数量，时间复杂度
O
(
n
5
)
O\left(n^5\right)
O(n5)。
记
f
(
x
,
y
,
e
)
f(x,y,e)
f(x,y,e)表示左边
x
x
x个点，右边
y
y
y个点，用了
e
e
e条边形成联通有标号二分图的方案数，这个也可以用容斥，枚举左边第一个点在哪个联通二分图内得到，时间复杂度
O
(
n
6
)
O\left(n^6\right)
O(n6)。
之后，我们就可以求出上面的那个
d
p
dp
dp值了。但是上面的哪个
d
p
dp
dp的状态是
n
4
n^4
n4的，暴力背包会达到
O
(
n
7
)
O\left(n^7\right)
O(n7)。
不过我们可以发现在计算答案的式子里与
c
c
c有关的只有
m
?
c
m^{-c}
m?c一项，我们不妨在背包的过程中将这个
m
?
c
m^{-c}
m?c加进去，就没必要在维护这一项了。
我们每次枚举在最小标号在已有二分图最小标号前面的连通图，将这个连通图加入，只需枚举其他点与已有二分图的位置关系即可。
这样就可以实现
O
(
n
6
)
O\left(n^6\right)
O(n6)的背包计算。
最后通过
d
p
dp
dp值计算我们的答案。
总时间复杂度 O ( n 6 ) O\left(n^6\right) O(n6)，看起来不太能过，不过由于存在大量的无用状态没必要转移，简单剪枝一下即可通过 n ? 40 n\leqslant 40 n?40的部分。
源码
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#include<bits/stdc++.h>
using namespace std;
#define MAXN 1055
#define lowbit(x) (x&-x)
#define reg register
#define pb push_back
#define mkpr make_pair
#define fir first
#define sec second
typedef long long LL;
typedef unsigned long long uLL;
typedef long double Ld;
typedef pair<int,int> pii;
const LL INF=0x3f3f3f3f3f3f3f3f;
const int mo=998244353;
const int mod=1e6+7;
const int inv2=499122177;
const int jzm=2333;
const int zero=100;
const int orG=3,ivG=332748118;
const double Pi=acos(-1.0);
const double eps=1e-9;
template<typename _T>
_T Fabs(_T x){return x<0?-x:x;}
template<typename _T>
void read(_T &x){
_T f=1;x=0;char s=getchar();
while(s>'9'||s<'0'){if(s=='-')f=-1;s=getchar();}
while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}
x*=f;
}
template<typename _T>
void print(_T x){if(x<0)putchar('-'),print(-x);if(x>9)print(x/10);putchar(x%10+'0');}
LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);}
int add(int x,int y,int p){return x+y<p?x+y:x+y-p;}
void Add(int &x,int y,int p){x=add(x,y,p);}
int qkpow(int a,int s,int p){int t=1;while(s){if(s&1)t=1ll*t*a%p;a=1ll*a*a%p;s>>=1;}return t;}
int n,m,P,fac[4005],inv[4005],ff[4005],Cm[45],c[45][45],im[85];
int t[45][45][45],f[45][45][45],g[45][45][45],ans,S[45][45];
void init(){
fac[0]=fac[1]=inv[0]=inv[1]=ff[1]=Cm[0]=1;
for(int i=2;i<=4000;i++)
fac[i]=1ll*i*fac[i-1]%P,
ff[i]=1ll*(P-P/i)*ff[P%i]%P,
inv[i]=1ll*inv[i-1]*ff[i]%P;
for(int i=1;i<=n;i++)
Cm[i]=1ll*(m-i+1)*ff[i]%P*Cm[i-1]%P;
for(int i=0;i<=n;i++){
c[i][0]=c[i][i]=1;
for(int j=1;j<i;j++)
c[i][j]=add(c[i-1][j-1],c[i-1][j],P);
}
S[0][0]=1;
for(int i=1;i<=n;i++)
for(int j=1;j<=i;j++)
S[i][j]=add(S[i-1][j-1],1ll*j*S[i-1][j]%P,P);
}
int C(int x,int y){
if(x<0||y<0||x<y)return 0;
return 1ll*fac[x]*inv[y]%P*inv[x-y]%P;
}
int main(){
read(n);read(m);read(P);init();
for(int x=1;x<=n;x++)
for(int y=1;y<=n;y++){
int d=max(x,y);for(int i=d;i<=n;i++)t[x][y][i]=C(x*y,i);
for(int j=0;j<=x;j++)for(int k=0;k<=y&&j+k<x+y;k++){
const int tmp=P-1ll*c[x][j]*c[y][k]%P;
for(int i=d;i<=n;i++)Add(t[x][y][i],1ll*tmp*t[j][k][i]%P,P);
}
for(int i=x+y-1;i<=n;i++)f[x][y][i]=t[x][y][i];
for(int j=1;j<=x;j++)for(int k=1;k<=y&&j+k<x+y;k++){
const int tmp=P-1ll*c[x-1][j-1]*c[y][k]%P;
for(int i=x+y-1;i<=n;i++)for(int r=j+k-1;r<=i;r++)
Add(f[x][y][i],1ll*tmp*f[j][k][r]%P*t[x-j][y-k][i-r]%P,P);
}
}
g[0][0][0]=1;const int ivm=qkpow(m,P-2,P);
im[0]=1;for(int i=1;i<=n+n;i++)im[i]=1ll*ivm*im[i-1]%P;
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)for(int k=i+j-1;k<=n;k++)f[i][j][k]=1ll*im[i+j-1]*f[i][j][k]%P;
for(int x=0;x<=n;x++)for(int y=0;y<=n;y++)
for(int i=1;x+i<=n;i++)for(int j=1;y+j<=n;j++){
const int tmp=1ll*c[x+i-1][i-1]*c[y+j][j]%P;
for(int z=max(x,y);z<=n;z++)for(int k=i+j-1;z+k<=n;k++)
Add(g[x+i][y+j][z+k],1ll*tmp*g[x][y][z]%P*f[i][j][k]%P,P);
}
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=max(i,j);k<=n;k++)
Add(ans,1ll*S[n][k]*fac[k]%P*Cm[i]%P*Cm[j]%P*g[i][j][k]%P,P);
ans=1ll*qkpow(qkpow(m,P-2,P),n+n,P)*ans%P;
printf("%d\n",ans);
return 0;
}
