长沙理工大学第十二届ACM大赛-重现赛 G

题目描述

vigoss18 辞职成功终于逃出了公司，但是没过太久，公司就发现vigoss18 的所作所为，于是派人来把他抓

回去。

vigoss18 必须一直跑路，躲避公司的围捕。可以抽象的看成一个有向图，图中可能存在重边和自环。

刚开始他站在位置1，每单位时间vigoss18 必须从目前站的位置，等概率选择一条边然后移动到对应的节

点上去或者不动（如果当前节点有t条边，则有1/(t+1)的概率选择一条边移动或者原地不动），可以认为每次需

要花费1 单位时间。

他就这样一直跑一直跑，过了很长很长的时间...

公司把你派出来寻找vigoss18，如果能抓到他，你将能升官发财赢取白富美走向人生巅峰。

但是你精力有限，不是太走的开身，所以写了一个程序，来计算vigoss18 在每个位置的概率，可以认为过

了很长时间以后，vigoss18 在每个位置的概率是收敛的。所以你需要告诉上司，他最可能在哪个位置(概率

最大的那个位置)。

你的上司并不想知道过程，他只想知道结果，所以你只需要告诉他这个概率最大是多少即可。

输入描述:

多组输入，保证绝大部分为小数据。
每组输入第一行n m(1<=n<=100,1<=m<=10000)，表示n个点m条有向边。
接下来m行，每行u v(1<=u,v<=n)，表示有一条有向边从u连向v

输出描述:

算出vigoss18在所有位置的概率，并输出其中的最大值即可。 
你的答案与标准答案的误差应保持在1e-6以内。

示例1

输入

3 3
1 2
2 3
3 1

输出

0.333333333

题解

$dp$。

$dp[i][j]$表示走$i$步，停在$j$的概率，$dp[i][j]$可以从$dp[i-1][*]$得到，拿矩阵快速幂跑跑就能得到比较高的精度。

#include <bits/stdc++.h>
using namespace std;
 
const int maxn = 100 + 10;
int g[maxn][maxn];
int out[maxn];
int n, m;
 
struct M {
  int r, c;
  double a[maxn][maxn];
};
 
M mul(const M &a, const M &b) {
  M res;
  res.r = a.r;
  res.c = b.c;
  for(int i = 1; i <= res.r; i ++) {
    for(int j = 1; j <= res.c; j ++) {
      res.a[i][j] = 0;
      for(int k = 1; k <= a.c; k ++) {
        res.a[i][j] = res.a[i][j] + a.a[i][k] * b.a[k][j];
      }
    }
  }
  return res;
}
 
int main() {
  while(~scanf("%d%d", &n, &m)) {
    for(int i = 1; i <= n; i ++) {
      out[i] = 0;
      for(int j = 1; j <= n; j ++) {
        g[i][j] = 0;
      }
    }
    while(m --) {
      int u, v;
      scanf("%d%d", &u, &v);
      g[u][v] ++;
      out[u] ++;
    }
 
    M A;
    A.r = n;
    A.c = n;
    for(int i = 1; i <= n; i ++) {
      for(int j = 1; j <= n; j ++) {
        if(i == j) A.a[i][j] = 1.0;
        else A.a[i][j] = 0.0;
      }
    }
 
    M B;
    B.r = 1;
    B.c = n;
    for(int j = 1; j <= n; j ++) {
      if(j == 1) B.a[1][j] = 1.0;
      else B.a[1][j] = 0.0;
    }
 
    M C;
    C.r = n;
    C.c = n;
    for(int j = 1; j <= n; j ++) {
      for(int i = 1; i <= n; i ++) {
        C.a[i][j] = 0.0;
        // i -> j
        if(i == j) {
          C.a[i][j] = 1.0 / (out[i] + 1);
        } else {
          C.a[i][j] = 1.0 * g[i][j] / (out[i] + 1);
        }
      }
    }
 
    int b = 0x7FFFFFFF;
    while(b) {
      if(b & 1) B = mul(B, C);
      b = b / 2;
      C = mul(C, C);
    }
    A = mul(A, B);
 
    double ans = 0.0;
    for(int j = 1; j <= n; j ++) {
      ans = max(ans, A.a[1][j]);
    }
    printf("%.8f
", ans);
 
  }
  return 0;
}