hzau1200——Choosy in Food(概率DP+高斯消元)

1200: Choosy in Food

Time Limit: 1 Sec   Memory Limit: 1280 MB
Submit: 23   Solved: 2
[ Submit][ Status][ Web Board]

Description

    Xiao Ming find n trees (labeled 0...n-1) which join a circle in a magical land. Each tree
is covered with ripe berries. Xiao Ming wants to find the tree whose berries are most delicious. He starts from the tree X, walks clockwise to find the target. Each time the probability of him walking k trees is Pk% . P1+P2+..+Pk = 100

      Now we know the target of Xiao Ming is the tree Y. Require to know expectation of the number of the trees which Xiao Ming walk through

Input

    The first parameter is T, indicates the case number.

    In each case, there are four integers at the first line, N, M, X , Y (Tree labels from 0 to N-1, N < 100, M < N).  

Output

     The expectation of the number of trees Xiao Ming walk through. It will be reserved two decimal places. If Xiao Ming will never arrive the tree, output -1. 

Sample Input

2 
4 2 0 1
50 50 
4 1 0 2 
100

Sample Output

4.20
2.00

华中农业大学校赛的一个题,跟hdu4818基本上一样，只不过现在n个点已经是一个环状了，所以不用把节点抽象成2n-2个，也不须考虑方向，绕着走就行，但这题没法交，就先把我认为对的代码发了..

题意：一个人在数轴上来回走，以pi的概率走i步i∈[1, m]，给定n（数轴长度），m，e（终点），s（起点），求从s走到e经过的点数期望
 
解析：设E[x]是人从x走到e经过点数的期望值，显然对于终点有：E[e] = 0
一般的：E[x] = sum((E[x+i]+i) * p[i])(i∈[1, m]) 
(走i步经过i个点，所以是E[x+i]+i)

首先按一般思路先列出f数组表示期望
显然f[e]=0,f[i]=∑((f[x+j]+j)*p[j])
然后把第二条式子变形，则得到：
f[i]-∑(f[x+i]*p[j])=∑(j*p[j])
接下来把所有的f值作为变量，高斯消元法求解所有方程就可以了。
于是总的思路如下：
1.      读入数据（注意如果起点就是终点的话直接输出0）
2.      BFS求出所有可达点
3.      根据上一步的结果列方程（当然如果终点不可达那么直接输出-1就行了）
4.      高斯消元法解方程

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <queue>
#include <algorithm>
#include <cmath>
using namespace std;
#define M 505
#define eps 1e-8
int equ, var;
double a[M][M], x[M];//a[i][j]是系数矩阵，x[i]是解=> AX=B,B=0
//增广列和解集共用的高斯消元，也可以用a[i][M+1]为增广列的高斯消元,x单独便是解集
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;  //eps可能是1e9更好点
        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;
}


//has[x]表示人在x点时的变量号，因为我们只用可达状态建立方程，所以需要编号
int has[M*M], vis[M*M], k, e, n, m;
double p[M*M], sum;


int bfs (int u)
{
    memset (has, -1, sizeof(has));
    memset (a, 0, sizeof(a));           //忘记初始化WA勒，以后得注意
    memset (vis, 0, sizeof(vis));
    int v, i, flg = 0;
    queue<int> q;
    q.push (u);
    k = 0;
    has[u] = k++;
    while (!q.empty ())
    {
        u = q.front ();
        q.pop();
        if (vis[u]) continue;
        vis[u] = 1;
        if (u == e )//||u == n-e)  //找到终点   //终点有两个，你懂的~
        {
            a[has[u]][has[u]] = 1;//1表示有相关性 1倍的E[e]=sum()=0
            x[has[u]] = 0;
            flg = 1;
            continue;
        }
        //状态转移方程E[x] = sum ((E[x+i]+i) * p[i])
        // ----> E[x] - sum(p[i]*E[x+i]) = sum(i*p[i])等号右边增广矩阵
        a[has[u]][has[u]] = 1; //1倍的E[x]=sum(i*p[i])=sum
        x[has[u]] = sum;
        for (i = 1; i <= m; i++)
        {
            //非常重要！概率为0，该状态可能无法到达，如果还去访问并建立方程会导致无解
            if (fabs (p[i]) < eps) continue;//因为是double类型，不能直接==0，指概率为0的点  跳过
            v = (u + i) % n;
            if (has[v] == -1) has[v] = k++;//has[v]相当于E[x+i]，走m次则一个方程E[x] - sum(p[i]*E[x+i])就出来了
            a[has[u]][has[v]] -= p[i];//E[x] - sum(p[i]*E[x+i])体现减法,我觉得也就是a[][]=-p[i]
            q.push (v);
        }
    }
    return flg;
}


int main()
{
    int t, s, d, i;
    scanf ("%d", &t);
    while (t--)
    {
        scanf ("%d%d%d%d", &n, &m, &s, &e);
        //n = 2*(n-1);
        sum = 0;
        for (i = 1; i <= m; i++)
        {
            scanf ("%lf", &p[i]);
            p[i] = p[i] / 100;
            sum += p[i] * i; //走i步的期望,计算sum(i*p[i])
        }
        if (s == e)
        {
            puts ("0.00");
            continue;
        }
        //一开始向左，起点要变
        // if (d > 0)
        //s = (n - s) % n;
        if (!bfs (s))
        {
            //puts ("-1");
            printf("-1\n");
            continue;
        }
        equ = var = k;
        if(Gauss())printf ("%.2f\n", x[has[s]]);
        else printf("-1\n");
    }
    return 0;
}