const double eps = 1e-10; double a[22 * 22][22 * 22], x[22 * 22]; //方程的左边的矩阵和等式右边的值,求解之后x存的就是结果 int equ, var; //方程数和未知数个数 int Gauss() { int i, j, k, col, max_r; for (k = 0, col = 0; k < equ && col < var; k++, col++) { max_r = k; for (i = k + 1; i < equ; i++) if (fabs(a[i][col]) > fabs(a[max_r][col])) max_r = i; if (fabs(a[max_r][col]) < eps) return 0; if (k != max_r) { for (j = col; j < var; j++) swap(a[k][j], a[max_r][j]); swap(x[k], x[max_r]); } x[k] /= a[k][col]; for (j = col + 1; j < var; j++) a[k][j] /= a[k][col]; a[k][col] = 1; for (i = 0; i < equ; i++) if (i != k) { x[i] -= x[k] * a[i][k]; for (j = col + 1; j < var; j++) a[i][j] -= a[k][j] * a[i][col]; a[i][col] = 0; } } return 1; } double P; int s[22 * 22][22 * 22]; void build() { CLR(a, 0); CLR(x, 0); FE(i, 0, 20) FE(j, i, 20) { int cur = s[i][j]; if (~cur) { a[cur][cur] = 1; if (i == 20 || j == 20) x[cur] = 0; else { int tx = min(i + 1, 20), ty = j; if (tx > ty) swap(tx, ty); int nxt = s[tx][ty]; a[cur][nxt] -= P; tx = max(i - 2, 0); ty = j; if (tx > ty) swap(tx, ty); nxt = s[tx][ty]; a[cur][nxt] -= 1 - P; x[cur] = 1; } } } } void bfs() { CLR(s, -1); queue<int> qx, qy; qx.push(0); qy.push(0); int cnt = 0; s[0][0] = cnt++; while (!qx.empty()) { int x = qx.front(); qx.pop(); int y = qy.front(); qy.pop(); int tx = min(x + 1, 20), ty = y; if (tx > ty) swap(tx, ty); if (!~s[tx][ty]) { s[tx][ty] = cnt++; qx.push(tx); qy.push(ty); } tx = max(x - 2, 0), ty = y; if (!~s[tx][ty]) { s[tx][ty] = cnt++; qx.push(tx); qy.push(ty); } } equ = var = cnt; } int main() { bfs(); while (cin >> P) { build(); Gauss(); printf("%.6lf\n", x[s[0][0]]); } return 0; }
原文地址:http://blog.csdn.net/wty__/article/details/38051001