HDU 4870 Rating

期望$dp$，高斯消元。

$dp[x][y]$表示一个账号$x$积分，另一个$y$积分的状态下到达目标状态需要的期望次数。

假设$x>=y$，则$dp[x][y]=p*dp[x][y+1]+(1-p)*dp[x][y-2]+1$。然后高斯消元即可。

一开始以为$1000/50$是$200$，不敢写了....后来发现是$20$......

高斯消元一直掉精度......后来把矩阵每一个系数都乘上$100000$就过了......

#include<cstdio>
#include<cstring>
#include<cmath>
#include<stack>
#include<queue>
#include<algorithm>
#include<iostream>
using namespace std;
const int mod=1e9+9;
const int N=1e5+10;

const double eps=1e-8;
int const maxn = 30;
double a[maxn][maxn],ans[maxn];
bool free_x[maxn];
int equ ,var ;

int sgn(double x)
{
    return (x>eps)-(x<-eps);
}

int gauss()
{
    int i, j, k;
    int max_r;
    int col;
    double temp;
    int free_x_num;
    int free_index;

    col = 0;
    memset(free_x,true,sizeof(free_x));
    memset(ans,0,sizeof ans);
    for (k = 0; k < equ && col < var; k++, col++)
    {
        max_r = k;
        for (i = k + 1; i < equ; i++)
        {
            if (sgn(fabs(a[i][col]) - fabs(a[max_r][col]))>0) max_r = i;
        }
        if (max_r != k)
        {
            for (j = k; j < var + 1; j++) swap(a[k][j], a[max_r][j]);
        }
        if (sgn(a[k][col]) == 0 )
        {
            k--; continue;
        }
        for (i = k + 1; i < equ; i++)
        {
            if (sgn(a[i][col])!=0)
            {
                double t = a[i][col] / a[k][col];
                for (j = col; j < var + 1; j++)
                {
                    a[i][j] = a[i][j] - a[k][j] * t;
                }
            }
        }
    }
    for(i=k;i<equ;i++)
    if(sgn(a[i][col])!=0) {return 0;}
    if (k < var)
    {
        for (i = k - 1; i >= 0; i--)
        {
            free_x_num = 0;
            for (j = 0; j < var; j++)
            {
                if ( sgn(a[i][j])!=0 && free_x[j]){
                    free_x_num++, free_index = j;
                }
            }
            if(free_x_num>1)    continue;
            temp = a[i][var];
            for (j = 0; j < var; j++)
            {
                if (sgn(a[i][j])!=0 && j != free_index) temp -= a[i][j] * ans[j];
            }
            ans[free_index] = temp / a[i][free_index];
            free_x[free_index] = 0;
        }
        return var - k;
    }

    for (i = var - 1; i >= 0; i--)
    {
        temp = a[i][var];
        for (j = i + 1; j < var; j++)
        {
            if (sgn(a[i][j])!=0) temp -= a[i][j] * ans[j];
        }
        ans[i] = temp / a[i][i];
    }
    return 1;
}

double p;
double dp[25][25];

int main()
{
    while(~scanf("%lf",&p))
    {
        for(int i=19;i>=0;i--)
        {
            memset(a,0,sizeof a);
            memset(ans,0,sizeof ans);

            equ=var=i+1;

            for(int j=0;j<=i;j++)
            {
                if(j==i)
                {
                    a[j][j]+=1;
                    a[j][max(0,j-2)]+=-(1-p);

                    a[j][i+1]+=1;
                    a[j][i+1]+=p*dp[i+1][i];
                }

                else if(j==0)
                {
                    a[j][0]+=1;
                    a[j][1]+=-p;
                    a[j][0]+=-(1-p);
                    a[j][i+1]=1;
                }
                else if(j==1)
                {
                    a[j][1]+=1;
                    a[j][2]+=-p;
                    a[j][0]+=-(1-p);
                    a[j][i+1]=1;
                }
                else if(j<i)
                {
                    a[j][j]+=1;
                    a[j][j+1]+=-p;
                    a[j][j-2]+=-(1-p);
                    a[j][i+1]=1;
                }

            }

            for(int f=0;f<=20;f++)
            {
                for(int g=0;g<=20;g++)
                {
                    a[f][g]=a[f][g]*100000;
                }
            }

            gauss();
            for(int j=0;j<=i;j++) dp[i][j]=ans[j];
        }
        printf("%.6f
",dp[0][0]);
    }
    return 0;
}