对状态转移方程的理解（【HNOI2013】游走&【hdu4035】Maze）

先来看一道例题：

【HNOI2013】游走

一个无向连通图，顶点从1编号到N，边从1编号到M。
小Z在该图上进行随机游走，初始时小Z在1号顶点，每一步小Z以相等的概率随机选 择当前顶点的某条边，沿着这条边走到下一个顶点，获得等于这条边的编号的分数。当小Z 到达N号顶点时游走结束，总分为所有获得的分数之和。
现在，请你对这M条边进行编号，使得小Z获得的总分的期望值最小。

【分析】

设$d[i]$为走上$i$号边的期望次数，则对$d[i]$从小到大排序后，有：
$$ans=\sum _{i=1}^{m}d[i]\ast (m-i+1)$$
关键是求$d[i]$。我们发现，通过一条边的次数与到达两顶点的有关。于是，设$f[i]$为通过$i$号点的期望次数，则（$d[i]$为$i$的度数）：
$$f[i]=\sum _{j\in son[i]}\frac{f[j]}{d[j]}+1$$
然而，由于图有环，转移顺序并不明确。所以说，我们把它当做一般的方程看待，其中$f[i]$、$f[j]$为未知数。
对应每一个$i\in [1,n-1]$，都有上面的方程。这样，我们就得到了$n-1$个方程组。
方程组？那就拿高斯消元乱搞一下不就好了？！

【代码】

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
const int mn = 505;
const double eps = 1e-6;
vector<int > g[mn];
double f[mn][mn], d[mn], b[mn];
struct edge{
    int a, b;
    double val;
    bool operator <(const edge x) const {return val < x.val;}
}e[mn * mn];
inline int dcmp(double x)
{
    if(fabs(x) < eps)   return 0;
    return x < 0 ? -1 : 1;
}
inline void gauss(int n)
{
    for(int i=1;i<=n;++i) {
        int p=i;
        for(int k=i+1;k<=n;++k) if(fabs(f[k][i])>fabs(f[p][i])) p=k;
        if(i!=p)
        {
            swap(b[i], b[p]);
            for(int j = 1; j <= n; j++)
                swap(f[i][j], f[p][j]);
        }
        for(int k=i+1;k<=n;++k) {
            double h=f[k][i]/f[i][i];
            b[k]-=h*b[i];
            for(int j=i;j<=n;++j) f[k][j]-=h*f[i][j];
        }
    }
    for(int i=n;i>=1;--i) {
        for(int j=i+1;j<=n;++j) b[i]-=d[j]*f[i][j];
        d[i]=b[i]/f[i][i];
    }
}
int main()
{
    int n, m, i, j, x, y;
    scanf("%d%d", &n, &m);
    for(i = 1; i <= m; i++)
        scanf("%d%d", &x, &y), g[x].push_back(y), g[y].push_back(x), e[i] = (edge) {x, y, 0};
    for(i = 1; i < n; i++)
    {
        f[i][i] = 1.0;
        int siz = g[i].size();
        for(j = 0; j < siz; j++)
        {
            int to = g[i][j];
            if(to != n) f[i][to] = -1.0 / (double)(g[to].size());
        }
    }
    b[1] = 1, gauss(n - 1);
    for(i = 1; i <= m; i++)
        e[i].val = d[e[i].a] / (double)(g[e[i].a].size()) + d[e[i].b] / (double)(g[e[i].b].size());
    sort(e + 1, e + 1 + m);
    double ans = 0;
    for(i = 1; i <= m; i++)
        ans += (m - i + 1) * e[i].val;
    printf("%.3lf\n", ans);
}

---

【hdu4035】Maze

给定一棵树，每一个节点有两个值$k_i$和$e_i$，表示在i号节点，有$k_i\%$的概率被传送至节点1，有$e_i\%$的概率走出这棵树。现在有一个人从1号节点出发，每到一个节点$i$，他会随机选择一条相邻边走，求他走出这棵树期望需要多少步。

【分析】

不难想到设$f[i]$表示从$i$出发走出迷宫的期望步数，则（$d[i]$表示$i$节点的度数）：
$$f[i]=\sum _{j\in son[i]}(f[j]+1)\ast \frac{(1-k[i]-e[i])}{d[i]}+(f[fa[i]]+1)\ast \frac{(1-k[i]-e[i])}{d[i]}+e[i]\ast 0+k[i]\ast f[1]$$
由于他可以随便乱走，所以我们依旧无法确定状态的转移顺序。所以如果不再变形，则只能用高斯消元。
然而，这道题的规模不允许我们这样搞（$n\le 10000$）。所以我们得在这个方程上搞些事情。
若按照树形dp的方式做，那么$f[fa[i]]$为未知状态，而$f[j],j\in son[i]$为已知状态。所以，我们设：
$$f[i]=a[i]\ast f[fa[i]]+b[i]\ast f[1]+c[i]$$
将i换为$j,j\in son[i]$，令$p[i]=\frac{(1-k[i]-e[i])}{d[i]}$，代入原方程：
$$f[i]=\sum _{j\in son[i]}(a[j]\ast f[i]+b[j]\ast f[1]+c[j]+1)\ast p[i]+(f[fa[i]]+1)\ast p[i]+f[1]\ast k[i]$$
令$t[i]=1-{\sum }_{j\in son[i]}a[j]\ast p[i]$，整理得：
$$f[i]=\frac{p[i]}{t[i]}f[fa[i]]+\frac{{\sum }_{j\in son[i]}b[j]\ast p[i]+k[i]}{t[i]}\ast f[1]+\frac{{\sum }_{j\in son[i]}c[j]\ast p[i]+p[i]}{t[i]}$$
所以：
$$\{\begin{array}{l}a[i]=\frac{p[i]}{t[i]}\\ b[i]=\frac{{\sum }_{j\in son[i]}b[j]\ast p[i]+k[i]}{t[i]}\\ c[i]=\frac{{\sum }_{j\in son[i]}c[j]\ast p[i]+p[i]}{t[i]}\end{array}$$
一次树形dp即可把a、b、c求出来。答案为：
$$f[1]=b[1]\ast f[1]+c[1]⟺f[1]=\frac{c[1]}{1-b[1]}$$
如果计算过程中出现了除0的情况，输出“impossible”。（若向某一个方向走的概率达到了1，则不可能走出树）
【代码】

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
const int mn = 10005;
const double eps = 1e-10;
vector<int >g[mn];
double a[mn], b[mn], c[mn], e[mn], k[mn], p[mn];
int d[mn];
bool flag;
inline void dp(int s, int fa)
{
    double tmp = 0;
    a[s] = p[s], b[s] = k[s], c[s] = 1 - k[s] - e[s];
    int m = g[s].size();
    for(int i = 0; i < m; i++)
    {
        int to = g[s][i];
        if(to != fa)
        {
            dp(to, s);
            if(flag)
                return;
            tmp += a[to] * p[s], b[s] += b[to] * p[s], c[s] += c[to] * p[s];
        }
    }
    tmp = 1.0 - tmp;
    if(fabs(tmp) < eps)
    {
        flag = 1;
        return;
    }
    a[s] /= tmp, b[s] /= tmp, c[s] /= tmp;
}
inline void init(int n)
{
    for(int i = 1; i <= n; i++)
        g[i].clear();
    flag = 0;
}
int main()
{
    int T, n, i, x, y;
    scanf("%d", &T);
    for(int tt = 1; tt <= T; tt++)
    {
        scanf("%d", &n), init(n);
        for(i = 1; i < n; i++)
            scanf("%d%d", &x, &y), g[x].push_back(y), g[y].push_back(x);
        for(i = 1; i <= n; i++)
            scanf("%lf%lf", &k[i], &e[i]), k[i] /= 100.0, e[i] /= 100.0, p[i] = 1 - k[i] - e[i];
        for(i = 1; i <= n; i++)
            d[i] = g[i].size();
        for(i = 1; i <= n; i++)
            p[i] /= (double)(d[i]);
        dp(1, 0), printf("Case %d: ", tt);
        if(flag || fabs(1 - b[1]) < eps)
            puts("impossible");
        else
            printf("%.6lf\n", c[1] / (1 - b[1]));
    }
}