[BZOJ 3640]JC的小苹果

逆矩阵文明，，，

很显然我们可以定义一个状态\(f[i][j]\)表示当前血量为\(i\)走到\(j\)的概率，然后肥肠爆歉这个东西没法DP（可能会有的点的伤害为0，这样可以凿出来环）。考虑高消，这个东西有个好处是不可能从血量低的到血量高的状态，所以可以从大到小枚举血量，这样各层事独立的，复杂度比直接高消降低了很少。可惜复杂度为\(O(sn^3)\)（设血量为\(s\)），会T掉。

考虑转移的过程，转移时等价于解这样一个方程：

\[
\begin{aligned}
a_{11}x_{1} + a_{12}x_{2} + \ldots + a_{1n}x_{n} &= c_1\\
a_{21}x_{1} + a_{22}x_{2} + \ldots + a_{2n}x_{n} &= c_2\\
&\ldots\\
a_{n1}x_{1} + a_{n2}x_{2} + \ldots + a_{nn}x_{n} &= c_n
\end{aligned}
\]

其中的未知数\(x\)事我们要求的东西，\(c\)表示从高血量状态转移过来的概率（这个可以视作常数）。根据友矩阵那一套理论，这一系列方程等价于以下等式：

\[
\begin{bmatrix}
a_{11} & a_{12} & \ldots & a_{1n}\\
a_{21} & a_{22} & \ldots & a_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
a_{n1} & a_{n2} & \ldots & a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_{1}\\ x_{2} \\ \vdots \\ x_{n}
\end{bmatrix}
=
\begin{bmatrix}
c_1\\ c_2\\ \vdots\\ c_n
\end{bmatrix}
\]

不妨将之简记为\(AB = C\)，其中我们只有\(B\)不知道，要求的也是\(B\)。我们除一下就可以得到\(B = A^{-1}C\)，然后我们还发现每一层的\(A\)都是一样的，所以每一层的\(A^{-1}\)也都是一样的，预处理即可。这样转移部分的复杂度就变成了\(O(sn^2)\)（矩阵乘法在这里事方阵乘列向量）。

至于逆矩阵的求法，我们知道对矩阵做初等变化也就等价于乘上另一个矩阵。因此，我们将一个矩阵\(A\)用类似于高消的手段消为单位阵\(I\)，所做的初等变换也就等价于乘上\(A^{-1}\)。我们对一个单位阵\(I\)作用上一样的操作，也就等于给这个单位阵乘上了\(A^{-1}\)，这样我们就得到了\(A^{-1}\)。

代码：

/**************************************************************
    Problem: 3640
    User: danihao123
    Language: C++
    Result: Accepted
    Time:8708 ms
    Memory:13416 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <cmath>
#include <algorithm>
#include <functional>
#include <utility>
#include <vector>
#include <queue>
#include <set>
#include <map>
int n, m, hp;
const int maxn = 151, maxm = 5005;
const int maxh = 10005;
typedef double R;
typedef R Mat[maxn][maxn];
Mat D, F;
void calc_inv() {
  for(int i = 1; i <= n; i ++) {
    F[i][i] = 1;
  }
  for(int i = 1; i <= n; i ++) {
    int r = i;
    for(int j = i + 1; j <= n; j ++) {
      if(fabs(D[j][i]) > fabs(D[r][i])) {
        r = j;
      }
    }
    // assert(r > -1);
    if(r != i) {
      for(int j = 1; j <= n; j ++) {
        std::swap(D[r][j], D[i][j]);
        std::swap(F[r][j], F[i][j]);
      }
    }
    R bs = D[i][i];
    for(int j = 1; j <= n; j ++) {
      D[i][j] /= bs; F[i][j] /= bs;
    }
    for(int k = 1; k <= n; k ++) {
      if(k != i) {
        R rate = D[k][i];
        for(int j = 1; j <= n; j ++) {
          D[k][j] -= rate * D[i][j];
          F[k][j] -= rate * F[i][j];
        }
      }
    }
  }
}
void matrix_mul(const Mat &A, const Mat &B, int a, int b, int c, Mat &res) {
  static Mat C;
  for(int i = 1; i <= a; i ++) {
    for(int j = 1; j <= c; j ++) {
      C[i][j] = 0;
    }
  }
  for(int i = 1; i <= a; i ++) {
    for(int j = 1; j <= c; j ++) {
      for(int k = 1; k <= b; k ++) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
  for(int i = 1; i <= a; i ++) {
    for(int j = 1; j <= c; j ++) {
      res[i][j] = C[i][j];
    }
  }
}
 
int first[maxn], deg[maxn];
int next[maxm << 1], to[maxm << 1];
int gcnt = 0;
void add_edge(int u, int v) {
  gcnt ++;
  next[gcnt] = first[u]; first[u] = gcnt;
  to[gcnt] = v;
}
void ins_edge(int u, int v) {
  deg[u] ++;
  add_edge(u, v);
  if(u != v) {
    deg[v] ++;
    add_edge(v, u);
  }
}
 
int atk[maxn];
R f[maxh][maxn];
R solve() {
  for(int i = 1; i <= n; i ++) {
    D[i][i] = 1.00;
  }
  for(int i = 1; i < n; i ++) {
    for(int j = first[i]; j; j = next[j]) {
      int v = to[j];
      if(!atk[v]) {
        D[v][i] -= 1.00 / (R(deg[i]));
      }
    }
  }
  calc_inv();
  R ans = 0; static Mat C;
  f[hp][1] = 1.00;
  for(int i = hp; i >= 1; i --) {
    for(int j = 1; j <= n; j ++) {
      C[j][1] = f[i][j];
    }
    matrix_mul(F, C, n, n, 1, C);
    for(int j = 1; j < n; j ++) {
      for(int e = first[j]; e; e = next[e]) {
        int v = to[e];
        if(atk[v] > 0 && i - atk[v] > 0) {
          f[i - atk[v]][v] += C[j][1] / (R(deg[j]));
        }
      }
    }
    ans += C[n][1];
  }
  return ans;
}
 
int main() {
  scanf("%d%d%d", &n, &m, &hp);
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &atk[i]);
  }
  for(int i = 1; i <= m; i ++) {
    int u, v; scanf("%d%d", &u, &v);
    ins_edge(u, v);
  }
  printf("%.8lf\n", solve());
  return 0;
}

题解 高斯消元 概率与期望 bzoj 矩阵乘法 友矩阵 逆矩阵 Comments(0) 2018年6月27日 20:12

[LibreOJ 2383][HNOI2013]游走

本野蛮人竟然没做过高消期望DP，，，泪，流了下来，，，

根据期望线性性，答案就是所有边的期望被走的次数乘上边的编号的和。一条边期望经过的次数可以根据他两个端点期望经过的次数来算（但是\(n\)要特判一下），要求所有点期望走过的次数当然就可以列\(n\)个方程然后高消力。然后期望走的次数多的边编号应该小，反之亦然，所以求完每条边走的次数的期望之后就贪心一下就好力。

代码：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <utility>
#include <cmath>
const int maxn = 505;
const int maxm = maxn * maxn;
using R = double;
const R eps = 1e-9;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    return ((x < 0.00) ? -1 : 1);
  }
}
R D[maxn][maxn];
int n;
void gauss() {
#ifdef LOCAL
  for(int i = 1; i <= n; i ++) {
    for(int j = 1; j <= n + 1; j ++) {
      printf("%.3lf ", D[i][j]);
    }
    puts("");
  }
#endif
  for(int i = 1; i <= n; i ++) {
    int r = i;
    for(int j = i + 1; j <= n; j ++) {
      if(fabs(D[j][i]) > fabs(D[r][i])) {
        r = j;
      }
    }
    assert(sign(D[r][i]) != 0);
    if(r != i) {
      for(int j = 1; j <= n + 1; j ++) {
        std::swap(D[i][j], D[r][j]);
      }
    }
    for(int k = i + 1; k <= n; k ++) {
      R rate = D[k][i] / D[i][i];
      for(int j = i; j <= n + 1; j ++) {
        D[k][j] -= D[i][j] * rate;
      }
    }
  }
  for(int i = n; i >= 1; i --) {
    for(int j = i + 1; j <= n; j ++) {
      D[i][n + 1] -= D[j][n + 1] * D[i][j];
      D[i][j] = 0;
    }
    D[i][n + 1] /= D[i][i]; D[i][i] = 1;
#ifdef LOCAL
    printf("E[%d] : %.3lf\n", i, D[i][n + 1]);
#endif
  }
}

int E[maxm][2], deg[maxn];
R tms[maxm];
void add_edge(int u, int v) {
  if(u != n) D[v][u] += 1.00 / (R(deg[u]));
}

int main() {
  int m; scanf("%d%d", &n, &m);
  for(int i = 1; i <= n; i ++) {
    D[i][i] = -1;
  }
  D[1][n + 1] = -1;
  for(int i = 1; i <= m; i ++) {
    scanf("%d%d", &E[i][0], &E[i][1]);
    deg[E[i][0]] ++; deg[E[i][1]] ++;
  }
  for(int i = 1; i <= m; i ++) {
    int u = E[i][0], v = E[i][1];
    add_edge(u, v); add_edge(v, u);
  }
  gauss();
  for(int i = 1; i <= m; i ++) {
    tms[i] = 0;
    int u = E[i][0], v = E[i][1];
    if(u != n) tms[i] += D[u][n + 1] / (R(deg[u]));
    if(v != n) tms[i] += D[v][n + 1] / (R(deg[v]));
  }
  std::sort(tms + 1, tms + 1 + m, [&](const R &i, const R &j) { return i > j; });
  R ans = 0;
  for(int i = 1; i <= m; i ++) {
    ans += tms[i] * (R(i));
  }
  printf("%.3lf\n", ans);
  return 0;
}

题解 概率与期望 高斯消元 loj HNOI Comments(0) 2018年6月27日 12:59

[BZOJ 5091][Lydsy0711月赛]摘苹果

秒，秒啊.jpg

首先根据期望线性性，每个点的期望可以分开算。不妨设\(f_{i, j}\)表示\(j\)在某个第\(i\)步被走到的概率。那么每个点\(j\)的期望访问次数就是\(\sum_{i = 0}^k f_{i, j}\)。

然后考虑去求那个\(f_{i, j}\)。显然\(f_{0, j} = \frac{d_j}{2m}\)。但是考虑\(f_{1, j}\)的转移方程：

\[f_{1, j} = \sum \frac{f_{0, u}}{d_u}\]

然后展开之后发现每个\(\frac{f_{0, u}}{d_u}\)都是\(\frac{1}{2m}\)，于是乎\(f_{1, j} = \frac{d_j}{2m}\)。

如此一来，对于任意\(i\)，都有\(d_{i, j} = \frac{d_j}{2m}\)，然后就很好做了……

代码：

/**************************************************************
    Problem: 5091
    User: danihao123
    Language: C++
    Result: Accepted
    Time:532 ms
    Memory:2384 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const ll ha = 1000000007LL;
ll pow_mod(ll a, ll b) {
  if(!b) return 1LL;
  ll ans = pow_mod(a, b >> 1);
  ans = (ans * ans) % ha;
  if(1LL & b) ans = (ans * a) % ha;
  return ans;
}
ll inv(ll v) {
  return pow_mod(v, ha - 2LL);
}
 
const int maxn = 100005;
ll a[maxn], d[maxn];
int main() {
  int n, m, k; scanf("%d%d%d", &n, &m, &k);
  for(int i = 1; i <= n; i ++) {
    scanf("%lld", &a[i]);
  }
  for(int i = 1; i <= m; i ++) {
    int u, v; scanf("%d%d", &u, &v);
    d[u] ++; d[v] ++;
  }
  ll ans = 0;
  ll inv_2m = inv(2 * m);
  for(int i = 1; i <= n; i ++) {
    ans += (((d[i] * inv_2m) % ha) * a[i]) % ha;
    if(ans >= ha) ans -= ha;
  }
  ans = (ans * (ll(k))) % ha;
  printf("%lld\n", ans);
  return 0;
}

题解 BZOJ 概率与期望 Comments(1) 2018年3月23日 16:32

[BZOJ 3925][ZJOI2015]地震后的幻想乡

赛艇，赛艇.jpg
首先这个问题的本质就是让你求一个边权在\([0, 1]\)间均匀随机分布的无向图的MST的最大边边权的期望……
有一个很经典的式子：
\[E = \int_0^1 p(\geq x)\,\mathrm{d}x\]
然后考虑那个\(p\)咋整。
首先对于每个包含1的点集（我们不妨假设是从点1开始扩展）\(S\)，式子可以这么写（这里不妨用\(T\)来表示两个点集间的边的数目）：
\[p_{S, x} = \sum_{1\in S_0 \subset S} (1 - x)^{T(S_0, S - S_0)}(1 - p_{S_0, x})\]
显然答案是\(\int_0^1 p_{S, x}\,\mathrm{d}x\)，但这玩咋求……
然后我们定义一个状态\(d\)：
\[d_{S, k} = \int_0^1 (1 - x)^k p_{S, x}\,\mathrm{d}x\]
把其中的\(p\)展开，整理一下，得到：
\[
\begin{aligned}
d_{S, k}=&\sum_{1\in S_0 \subset S}(\int_0^1(1 - x)^{T(S, S - S_0) + k}\,\mathrm{d}x\\
& - \int_0^1(1 - x)^{T(S, S - S_0) + k}\,p_{S_0,\,x}\,\mathrm{d}x)
\end{aligned}
\]
里面有两大块定积分，前一块还看起来蛮好求的，但后面一块……
不就是\(d_{S_0, T(S, S - S_0) + k}\) 吗？
然后这样整个$d$就可以搞一波状压DP了。
然后观察\(d_{S, 0}\)（这里不妨假设S为全集）：
\[d_{S, 0}=\int_0^1(1 - x)^0 p_{S, x}\,\mathrm{d}x = \int_0^1 p_{S, x}\,\mathrm{d}x\]
这不就是我们要的答案吗？
至于边界处理，\(d_{\{1\}, x}\)全部搞成0就行。
代码：
/**************************************************************
    Problem: 3925
    User: danihao123
    Language: C++
    Result: Accepted
    Time:100 ms
    Memory:1272 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <bitset>
typedef double R;
const int maxn = 11;
const int maxm = 50;
int edge[maxm][2];
 
int n, m;
R d[1 << 10][maxm];
bool vis[1 << 10][maxm];
R f(int s, int k) {
  if(s == 1) return 0;
  if(vis[s][k]) return d[s][k];
  d[s][k] = 0;
  int lit_s = s >> 1;
  for(int lit_s0 = (lit_s - 1) & lit_s; ; lit_s0 = (lit_s0 - 1) & lit_s) {
    int s0 = lit_s0 * 2 + 1;
    int t = 0;
    for(int i = 1; i <= m; i ++) {
      int u = edge[i][0], v = edge[i][1];
      if(((1 << u) & s) == 0 || ((1 << v) & s) == 0) continue;
      if((((1 << u) & s0) == 0) ^ (((1 << v) & s0) == 0)) {
        t ++;
      }
    }
    int z = k + t;
    d[s][k] += 1.00 / ((double(z)) + 1.00) - f(s0, z);
    if(s0 == 1) break;
  }
  vis[s][k] = true;
  return d[s][k];
}
 
int main() {
  scanf("%d%d", &n, &m);
  for(int i = 1; i <= m; i ++) {
    int u, v; scanf("%d%d", &u, &v); u --, v --;
    edge[i][0] = u, edge[i][1] = v;
  }
  printf("%.6lf\n", f((1 << n) - 1, 0));
  return 0;
}

 

题解 数学分析 ZJOI bzoj 状压dp 概率与期望 连续型期望 Comments(1) 2018年2月13日 18:00