HDU 4118Time travel 概率dp+gauss消元

参考：

http://www.cnblogs.com/kuangbin/archive/2012/10/06/2713432.html

http://www.cnblogs.com/jffifa/archive/2012/09/24/2700707.html

http://blog.csdn.net/ok_again/article/details/11107719

/*
HDU 4118
题目：给出一个数轴，有一个起点和一个终点，某个人可以
走1，2，3……m步，每一种情况有一个概率，初始有一个方向，
走到头则返回，问到达终点的期望步数为多少。

比较明显的高斯求期望问题

把N个点转成2*N-2个点。
然后就是高斯消元法求解概率DP了。
*/

const double eps = 1e-9;
const int MAXN=220;
double a[MAXN][MAXN],x[MAXN];///方程的左边的矩阵和等式右边的值，求解之后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;
}


int num[MAXN];///标记能否到达
int cnt;///从s开始能到达的点数即：方程数和未知数个数，(1)转向(2)概率
///cnt从0开始
int n, N;

int  s, e, D;
int M;
double p[MAXN];

///建立方程组（5）
void build_gauss(int cnt, int M, int num[], double p[])
{
    equ = var = cnt;
    CLR(a, 0);
    CLR(x, 0);
    for (int i = 0; i < n; i++)
        if (num[i] != -1)
        {
            if (i == e || i == n - e)///终点时
            {
                a[num[i]][num[i]] = 1;
                x[num[i]] = 0;
                continue;
            }

            a[num[i]][num[i]] = 1;
            for (int j = 1; j <= M; j++)
            {
                int r = (i + j) % n;
                if (num[r] != -1)
                {
                    a[num[i]][num[r]] -= p[j];
                    x[num[i]] += j * p[j];
                }
            }
        }
}

///(1)转向(2)概率
void bfs(int s)
{
    CLR(num, -1);
    queue<int>q;
    cnt = 0;
    num[s] = cnt++;
    q.push(s);
    while (!q.empty())
    {
        int u = q.front(); q.pop();
        for (int i = 1; i <= M; i++)
        {
            if (fabs(p[i]) < eps) continue; ///（4）这点很重要，这个想到不能达到的点,????

            int r = (u + i) % n;
            if (num[r] == -1)
            {
                num[r] = cnt++;
                q.push(r);
            }
        }
    }
}

int main ()
{
    int T;
    scanf("%d", &T);
    while (T--)
    {
        RII(N, M);
        RIII(e, s, D);
        for (int i = 1; i <= M; i++)
        {
            scanf("%lf", &p[i]);
            p[i] /= 100;
        }
        if (e == s) ///（1）特判,这个特判一定需要，否则可能N==1,会被0除，RE,?????
        {
            printf("0.00\n");
            continue;
        }

        n = 2 * (N - 1);///0~(2 * N - 3)
        if (D == 1)///处理D
            s = n - s;
        bfs(s);///(2)

        if (num[e] == -1 && num[n - e] == -1)///(3)终点不能到达
        {
            printf("Impossible !\n");
            continue;
        }

        build_gauss(cnt, M, num, p);

        if (Gauss()) printf("%.2lf\n", x[num[s]]);
        else printf("Impossible !\n");///返回0，表示无解
    }
    return 0;
}


/*
HDU 4118
题目：给出一个数轴，有一个起点和一个终点，某个人可以
走1，2，3……m步，每一种情况有一个概率，初始有一个方向，
走到头则返回，问到达终点的期望步数为多少。

比较明显的高斯求期望问题

把N个点转成2*N-2个点。
然后就是高斯消元法求解概率DP了。
*/


const double eps = 1e-9;
const int MAXN=220;
double a[MAXN][MAXN + 10];///方程的左边的矩阵和最后一列即cnt列是等式右边的值，求解之后最后一列即cnt列x存的就是结果

inline int sgn(double d) {
    if (fabs(d) < eps) return 0;
    return d > 0 ? 1 : -1;
}

int Gauss(int N, int M){
    int i, j, r, c, pvt;
    double maxp;
    for (r = 0, c = 0; r < N && c < M; ++ r, ++ c) {
        for (maxp = 0, i = r; i < N; ++ i)
            if (fabs(a[i][c])>fabs(maxp)) maxp = a[pvt=i][c];
        if (sgn(maxp) == 0) {
            r--;
            continue;
        }
        if (pvt != r)
            for (j = r; j <= M; ++j) swap(a[r][j], a[pvt][j]);
        for (j = c+1; j <= M; ++j) {
            a[r][j] /= maxp;
            for (i = r+1; i < N; ++i)
                a[i][j] -= a[i][c]*a[r][j];
        }
    }
    for (i = r; i < N; ++i)
        if (sgn(a[i][M])) return -1;///无解
    if (r < M) return M-r;///自由变量个数
    for (i = M-1; i >= 0; --i)
        for (j = i+1; j < M; ++j)
            a[i][M] -= a[j][M]*a[i][j];
    return 0;///表示有唯一解
}


int num[MAXN];///标记能否到达
int cnt;///从s开始能到达的点数即：方程数和未知数个数，(1)转向(2)概率
///cnt从0开始
int n, N;

int  s, e, D;
int M;
double p[MAXN];

///建立方程组（5）
void build_gauss(int cnt, int M, int num[], double p[])///cnt和num[]， M和p[]
{
    CLR(a, 0);
    for (int i = 0; i < n; i++)
        if (num[i] != -1)
        {
            if (i == e || i == n - e)///终点时
            {
                a[num[i]][num[i]] = 1;
                a[num[i]][cnt] = 0;
                continue;
            }

            a[num[i]][num[i]] = 1;
            for (int j = 1; j <= M; j++)
            {
                int r = (i + j) % n;
                if (num[r] != -1)
                {
                    a[num[i]][num[r]] -= p[j];
                    a[num[i]][cnt] += j * p[j];
                }
            }
        }
}

///(1)转向(2)概率
void bfs(int s)
{
    CLR(num, -1);
    queue<int>q;
    cnt = 0;
    num[s] = cnt++;
    q.push(s);
    while (!q.empty())
    {
        int u = q.front(); q.pop();
        for (int i = 1; i <= M; i++)
        {
            if (fabs(p[i]) < eps) continue; ///（4）这点很重要，这个想到不能达到的点,????

            int r = (u + i) % n;
            if (num[r] == -1)
            {
                num[r] = cnt++;
                q.push(r);
            }
        }
    }
}

int main ()
{
    int T;
    scanf("%d", &T);
    while (T--)
    {
        RII(N, M);
        RIII(e, s, D);
        for (int i = 1; i <= M; i++)
        {
            scanf("%lf", &p[i]);
            p[i] /= 100;
        }
        if (e == s) ///（1）特判,这个特判一定需要，否则可能N==1,会被0除，RE,?????
        {
            printf("0.00\n");
            continue;
        }

        n = 2 * (N - 1);///0~(2 * N - 3)
        if (D == 1)///处理D
            s = n - s;
        bfs(s);///(2)

        if (num[e] == -1 && num[n - e] == -1)///(3)终点不能到达
        {
            printf("Impossible !\n");
            continue;
        }

        build_gauss(cnt, M, num, p);

        if (Gauss(cnt, cnt) == 0) printf("%.2lf\n", a[num[s]][cnt]);
        else printf("Impossible !\n");///返回0，表示无解
    }
    return 0;
}