HDU4870:Rating(DP+高斯消元)

Problem Description
A little girl loves programming competition very much. Recently, she has found a new kind of programming competition named "TopTopTopCoder". Every user who has registered in "TopTopTopCoder" system will have a rating, and the initial value of rating equals to zero. After the user participates in the contest held by "TopTopTopCoder", her/his rating will be updated depending on her/his rank. Supposing that her/his current rating is X, if her/his rank is between on 1-200 after contest, her/his rating will be min(X+50,1000). Her/His rating will be max(X-100,0) otherwise. To reach 1000 points as soon as possible, this little girl registered two accounts. She uses the account with less rating in each contest. The possibility of her rank between on 1 - 200 is P for every contest. Can you tell her how many contests she needs to participate in to make one of her account ratings reach 1000 points?
 

Input
There are several test cases. Each test case is a single line containing a float number P (0.3 <= P <= 1.0). The meaning of P is described above.
 

Output
You should output a float number for each test case, indicating the expected count of contest she needs to participate in. This problem is special judged. The relative error less than 1e-5 will be accepted.
 

Sample Input
1.000000 0.814700
 

Sample Output
39.000000 82.181160

SOLUTION:

这是题解上给的描述：令(x, y)表示高分为x，低分为y的状态（x >= y），E(x, y)表示从(x, y)到达(1000, ?)的比赛场数期望。容易得到E(x, y) = P * E(x1, y1) + (1 - P) * E(x2, y2) + 1，其中，(x1, y1)表示rating上升后的状态，（x2, y2）表示rating下降后的状态。把E(1000, ?) = 0带入可以得到包含n个未知数的n个方程，n大概200多，可以高斯消元。E(0, 0)即为答案。

题解描述的很详细啊。这里因为每次的操作是涨50，或者掉100分。所以我们以50作为划分可以把1000分分为20段，会产生一共(21*20)/2 = 210种状态。

根据题意我们可以得到：dp[x][y] = p*dp[x1][y1] + (1-p)*dp[x2][y2] + 1。所以有：dp[x][y] - p*dp[x1][y1] - (1-p)*dp[x2][y2] = 1。就可以构造出系数矩阵了，然后用高斯消元求出dp[0][0]就是答案。

还有就是要注意精度，我这么写必须eps<= 10^-10。

削成对角线的形式

CODE:

#include <algorithm>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <string>
#include <queue>
#include <cmath>
#include <stack>
#include <map>
#include <set>
#define eps 1e-10
///#define M 1000100
#define LL __int64
///#define LL long long
#define INF 0x7fffffff
#define PI 3.1415926535898
#define zero(x) ((fabs(x)<eps)?0:x)

const int maxn = 240;

using namespace std;

double a[maxn][maxn];
double x[maxn];
int link[maxn][maxn];
int cnt;
int equ, var;

double Gauss()
{
    int row, col, max_r;
    row = col = 0;
    while(row < equ && col < var)
    {
        max_r = row;
        for(int i = row+1; i < equ; i++)
            if(fabs(a[i][col])-fabs(a[max_r][col]) > eps)
            max_r = i;
        if(max_r != row)
            for(int j = col; j <= var; j++)
            swap(a[row][j], a[max_r][j]);
        if(fabs(a[row][col]) < eps)
        {
            col++;
            continue;
        }
        for(int i = 0; i < equ; i++)if(i!=row)
        {
            if(fabs(a[i][col]) > eps)
            {
                double t = a[i][col]/a[row][col];
                a[i][col] = 0.0;
                for(int j = col+1; j <= var; j++)
                    a[i][j] -= a[row][j]*t;
            }
        }
        row++;
        col++;
    }
   /* for(int i = equ-1; i >= 0; i--)
    {
        if(fabs(a[i][i]) < eps) continue;
        double tmp = a[i][var];
        for(int j = i+1; j < var; j++)
            tmp -= a[i][j]*x[j];
        x[i] = tmp/a[i][i];
    }*/
   // cout<<a[0][210]<<" "<<a[0][0]<<endl;
    return a[0][210]/a[0][0] ;
}

void Del()
{
    equ = var = 210;
    cnt = 0;
    memset(link, -1,  sizeof(link));
    for(int i = 0; i < 20; i++)
        for(int j = 0; j <= i; j++)
        link[i][j] = cnt++;
}

int main()
{
    Del();
    double p;
    while(~scanf("%lf",&p))
    {
        memset(a, 0, sizeof(a));
        int x;
        for(int i = 0; i < 20; i++)
        {
            for(int j = 0; j < i; j++)
            {
                x = link[i][j];
                a[x][x] = 1;
                a[x][210] = 1;
                a[x][link[i][j+1]] -= p;
                a[x][link[i][max(0, j-2)]] -= (1-p);
            }
            x = link[i][i];
            a[x][x] = 1;
            a[x][210] = 1;
            a[x][link[i][max(0,i-2)]] -= (1-p);
            a[x][link[i+1][i]] -= p;
        }
        printf("%.6f
",Gauss());
    }
    return 0;
}