hihoCoder #1047 Random Tree

题意

给出点数为 \(n\)（\(n \le 1000\)）的完全图 \(K_n\)，带边权。随机出 \(K_n\) 的一棵生成树 \(T\)。求 \(T\) 上任意两点间距离的期望。

解法

固定两点 \(u\)、\(v\)（\(u \le v\)），考虑生成树 \(T\) 上 \(u\) 到 \(v\) 的路径 \(P_{uv}\)。\(P_{uv}\) 上的边可分成三类：

1. \((u, v)\)
2. \((u, x)\)、\((y, v)\)，\(x,y \notin \\{u, v\\}\)
3. \((x,y)\)，\(x, y \notin \\{u, v\\}\)

第1类边出现在 \(P_{uv}\) 上的概率为 \(\dfrac{2}{n}\)
每个第2类边出现在 \(P_{uv}\) 上的概率为 \(\dfrac{1-\dfrac{2}{n}}{n-2}=\dfrac{1}{n}\)

考虑第3类边（对期望）的贡献。

首先应当注意到，所有第3类边出现在 \(P_{uv}\) 上是等可能的，所以我们只需要求 \(P_{uv}\) 上第三类边的数目的期望 \(E(n)\)。

用 \(f(i)\) 表示 \(K_n\) 的所有生成树中，满足「\(P_{uv}\) 上点数为 \(i\)（包括两端点 \(u\)，\(v\)）」的生成树的数目。
我们分 3 步来求 \(f(i)\)：

1. 固定 \(P_{u,v}\)，将 \(P_{uv}\) 缩成一点 \(w\)，加上余下的 \(n-i\) 个点，就得到一棵 \(n-i+1\) 个点的树 \(T'\)。

2. 将 \(w\) 的度数固定为 \(j\)，对应的生成树 \(T'\) 的数目 \(g(j)\) 的表达式为
   \begin{equation}
   g(j) = \binom{n-i-1}{j-1}(n-i)^{n-i-j} \label{E:1}
   \end{equation}
   \(\eqref{E:1}~\)式可通过 Prufer 序列与树的一一对应关系得到。

3. 与 \(w\) 相连的 \(j\) 棵子树中的每一棵，在 \(T\) 中可以连在 \(P_{uv}\) 上的 \(i\) 个点中的任意一个，所以我们得到

\[\begin{equation} \begin{aligned} f(i) &= \mathrm{A}\_{n-2}^{i-2}\sum\_{j=1}^{n-i} g(j) \cdot i^{j} \\\\ &= \mathrm{A}\_{n-2}^{i-2}\sum\_{j=1}^{n-i} \binom{n-i-1}{j-1} (n-i)^{n-i-j} \cdot i^{j} \\\\ &= \mathrm{A}\_{n-2}^{i-2}\cdot i \cdot \sum\_{j'=0}^{n-i-1} \binom{n-i-1}{j'}(n-i)^{n-i-1-j'} \cdot i^{j'} \\\\ &= \mathrm{A}\_{n-2}^{i-2} \cdot i \cdot n^{n-i-1} \label{E:2} \end{aligned} \end{equation} \]

从而

\[\begin{equation} \begin{aligned} E(n) &= \frac{\sum\limits_{i=4}^{n} f(i)(i-3)}{n^{n-2}} \\\\ &= \sum_{i=4}^{n} \frac{\mathrm{A}_{n-2}^{i-2}i(i-3)}{n^{i-1}} \end{aligned} \end{equation} \]

Implementation

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

using DB=long double;

const int N=1005;
DB res[N][N];
int a[N][N];


DB calc(int n){
    if(n<=3) return 0;
    DB pn=1;
    for(int i=1; i<=n-2; i++)
        pn*=i, pn/=n;
    // cout << pn << endl;
    DB sum=pn*(n-3);
    for(int i=n-1; i>=4; i--)
        pn*=n*i, pn/=(n-i)*(i+1), sum+=pn*(i-3);
    return sum;
}

int main(){
    // int cnt=0;
    // for(int i=0; i<=1000; i++)
    //     cnt+=fabs(t[i]-calc(i))>1e-50;
    // cout << cnt << endl;

    int n, tot=0;
    scanf("%d", &n);
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
            scanf("%d", a[i]+j), a[i][0]+=a[i][j], tot+=a[i][j];

    tot/=2;
    DB x=calc(n);
    for(int i=1; i<n; i++)
        for(int j=i+1; j<=n; j++){
            res[i][j]=(a[i][0]+a[j][0])/DB(n);
            if(n>=4)    // 注意：n=2 或 3 时，分母为 0
                res[i][j]+=x*(tot-a[i][0]-a[j][0]+a[i][j])/((n-2)*(n-3)/2);
            res[j][i]=res[i][j];
        }

    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
            printf("%.9Lf%c", res[i][j], j==n?'\n':' ');

    return 0;
}