2014 Multi-University Training Contest 1 - HDU4870

本文介绍了一种利用高斯消元方法求解特定概率期望问题的技术。通过建立递推方程并转化为线性方程组，采用高斯消元算法高效地解决了存在自反馈的概率转移问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

高斯消元+概率期望基本上是照着标程写的，发现原来用的高斯消元模板解不出来以后就用这个吧

推导过程：

对于这样的状态表示从当前rating x 和 y 递推到的期望。

因此我们很容易想到递推方程：

$\begin{cases} D(x,y) = p\cdot D(max(x, y + 1), min(x, y + 1)) + (1-p)\cdot D(x,y-1) \\D(20,x)=0 ~(x = 20) \end{cases}$

虽然能写出递推方程但是存在能推回自身的情况，所以不好做动态规划和搜索（当然搜索可以得到答案但是会漏很多情况而且会爆内存）。

所以我们可以得到一个方程组：

$\begin{cases} -D(x,y)+p \cdot D(x_1,y_1)+(1-p) \cdot D(x_2,y_2)=-1 \\ -D(x,y)'+p \cdot D(x_1,y_1)'+(1-p) \cdot D(x_2,y_2)'=-1 \\ -D(x,y)''+p \cdot D(x_1,y_1)''+(1-p) \cdot D(x_2,y_2)''=-1 \\~~~~~~~.~.~.~.~.~. \end{cases}$

可以得到增广矩阵

$\begin{pmatrix} &-1 &... &p &... &q &... &-1 \\ &-1 &... &p &... &q &... &-1 \\ &-1 &... &p &... &q &... &-1 \\ &... &... &... &... &...&...&... \end{pmatrix}$

然后进行高斯消元得到所有未知数的解 X0 为所求 D(0,0)

#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<string>
#include<vector>
#include<map>
#include <cmath>
#include<cstring>
#include <cassert>
#define maxn 220
#define INF (1e9)
#define eps 1e-10
typedef long long ll;
using namespace std;
double a[maxn][maxn], b[maxn];
int r[maxn], c[maxn];
void gauss_elimination(int n, double a[maxn][maxn], int r[maxn], int c[maxn]) {
    for (int i = 0; i < n; ++i) {
        r[i] = c[i] = i;
    }
    for (int k = 0; k < n; ++k) {
        int ii = k, jj = k;
        for (int i = k; i < n; ++i) {
            for (int j = k; j < n; ++j) {
                if (fabs(a[i][j]) > fabs(a[ii][jj])) {
                    ii = i;
                    jj = j;
                }
            }
        }
        swap(r[k], r[ii]);
        swap(c[k], c[jj]);
        for (int i = 0; i < n; ++i) {
            swap(a[i][k], a[i][jj]);
        }
        for (int j = 0; j < n; ++j) {
            swap(a[k][j], a[ii][j]);
        }
        if (fabs(a[k][k]) < eps) {
            continue;
        }
        for (int i = k + 1; i < n; ++i) {
            a[i][k] = a[i][k] / a[k][k];
            for (int j = k + 1; j < n; ++j) {
                a[i][j] -= a[i][k] * a[k][j];
            }
        }
    }
}
void solve(int n, double a[maxn][maxn], int r[maxn], int c[maxn], double b[maxn]) {
    static double x[maxn];
    for (int i = 0; i < n; ++i) {
        x[i] = b[r[i]];
    }
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < i; ++j) {
            x[i] -= a[i][j] * x[j];
        }
    }
    for (int i = n - 1; i >= 0; --i) {
        for (int j = n - 1; j > i; --j) {
            x[i] -= a[i][j] * x[j];
        }
        if (fabs(a[i][i]) >= eps) {
            x[i] /= a[i][i];
        }  else assert(fabs(x[i]) < eps);
    }
    for (int i = 0; i < n; ++i) {
        b[c[i]] = x[i];
    }
}
int pp = 0;
map<pair<int,int>, int> Ma;
void dfs(int x, int y) { // x >= y
    if (x != 20 && !Ma.count(make_pair(x,y))) {
        Ma[make_pair(x,y)] = pp++;
        //cout << x << " " << y << " " << pp << endl;
        if (min(y+1,20) >= x) {
            dfs(y+1,x);
        }
        else if (min(y+1,20) < x) {
            dfs(x, y+1);
        }
        dfs(x, max(y-2,0));
    }
}
int doit(int x, int y) {
    int Max = max(x, y);
    int Min = min(x, y);
    if (Max == 20) return -1;
    return Ma[make_pair(Max,Min)];
}
int main() {
    double p;
    dfs(0,0);
    while (scanf("%lf", &p) != EOF) {
        memset(a, 0, sizeof(a));
        for (__typeof(Ma.begin()) it = Ma.begin(); it != Ma.end(); ++ it) {
            int x = (it->first).first;
            int y = (it->first).second;
            int from = it->second;
            int to = doit(x, min(y+1,20));
            if (to != -1) a[from][to] = p;
            to = doit(x, max(y-2,0));
            if (to != -1) a[from][to] = 1-p;
            a[from][from] -= 1;
            b[from] = -1;
        }
        gauss_elimination(pp,a,r,c);
        solve(pp,a,r,c,b);
        printf("%.6lf\n", b[0]);
    }
}