HDU 4870 Rating（高斯消元 ）

HDU 4870   Rating

这是前几天多校的题目，高了好久突然听旁边的大神推出来说是可以用高斯消元，一直喊着赶快敲模板，对于从来没有接触过高斯消元的我来说根本就是一头雾水，无赖之下这几天做DP，正好又做到了这个题，没办法得从头开始看，后来在网上找了别人的高斯消元的模板后发现其实也还是很好理解，就是先构造一个增广矩阵，然后化行阶梯形，最后迭代求解

首先有一个介绍高斯消元感觉过于详细的博客http://blog.csdn.net/tsaid/article/details/7329301

首先看一下这个题怎么构造这个增广矩阵，我们把所有可以达到的分数组合作为一个点，再考虑它与其他点所连的边，例如：

(300, 200)<--(1-p)----(300, 300)-----(p)----->(350, 300)

a                               b                                    c

我们就可以理解为有一条b--->a的权值为（1-p）的边，有一条b---->c的权值为p的边，那这样首先构造状态转移方程：

DP[b] = 1 + (1-p) * DP[a] + p * DP[c]

变形后得到:

DP[b] - p * DP[c] - (1-p) * DP[a] = 1;

这就是我们的增广矩阵的系数了，对于方程的一般形式Ax = B，可以理解为第b个方程中的变元的系数为A[b][b] = 1, A[b][c] = p, A[b][a] = (1-p),B[b] = 1

这样就构造出了一个（A，B）的一个增广矩阵，保存在a中

然后就是高斯消元化行阶梯行，看看代码很好理解，就只有两个操作r1<-->r2，交换两行，r2 = r2 - r1 * a   (其中a为一个常系数)

第一次写Gauss 代码比较戳，可以根据网上详细的介绍结合代码看，很容易懂的

void gauss()

{

        int col = ;

        for(int k=;k<cnt && col < cnt;k++, col ++)

        {

                    double Max = fabs(a[k][col]);

                    int  Maxr = k;

                    for(int r = k + ; r < cnt; r ++)

                            if( fabs(a[r][col])  -  Max  >  eps  )

                                    Max  =  fabs( a[Maxr = r][col] );

                    if(fabs(Max) < eps)  {  k --;  continue;  }

                    for(int c = col;c<=cnt; c ++ )

                            SWAP(a[Maxr][c],  a[k][c]  );

                    for(int r = k + ; r < cnt; r ++ ) if( fabs(a[r][col]) > eps  )

                    {

                            double tmp = a[r][col] / a[k][col];

                            for(int c = col; c <= cnt; c ++ )

                                    a[r][c] -= tmp * a[k][c];

                    }

        }

}

化行阶梯行后就是迭代求解了，由于这里一定有解，所以不需要判断无解的情况，直接算就是了

        for(int r=cnt-;r>=;r--)

         {

                     for(int c = cnt-;c>r;c--)

                             a[r][cnt] -= a[r][c] * ans[c];

                     ans[r] = a[r][cnt] / a[r][r];

         }

最后的总复杂度就是O(n^3)n接近200

另外还有一个复杂度O(n)的解法的思路（N=20）http://www.cnblogs.com/gj-Acit/p/3888390.html

 #include <map>

 #include <set>

 #include <stack>

 #include <queue>

 #include <cmath>

 #include <ctime>

 #include <vector>

 #include <cstdio>

 #include <cctype>

 #include <cstring>

 #include <cstdlib>

 #include <iostream>

 #include <algorithm>

 using namespace std;

 #define INF ((LL)100000000000000000)

 #define inf (-((LL)1<<40))

 #define lson k<<1, L, mid

 #define rson k<<1|1, mid+1, R

 #define mem0(a) memset(a,0,sizeof(a))

 #define mem1(a) memset(a,-1,sizeof(a))

 #define mem(a, b) memset(a, b, sizeof(a))

 #define FOPENIN(IN) freopen(IN, "r", stdin)

 #define FOPENOUT(OUT) freopen(OUT, "w", stdout)

 template<class T> T CMP_MIN ( T a, T b ) { return a < b;   }

 template<class T> T CMP_MAX ( T a, T b ) { return a > b;   }

 template<class T> T MAX ( T a, T b ) { return a > b ? a : b; }

 template<class T> T MIN ( T a, T b ) { return a < b ? a : b; }

 template<class T> T GCD ( T a, T b ) { return b ? GCD ( b, a % b ) : a; }

 template<class T> T LCM ( T a, T b ) { return a / GCD ( a, b ) * b;       }

 template<class T> void SWAP( T& a, T& b ) { T t = a; a = b;  b = t; }

 //typedef __int64 LL;

 typedef long long LL;

 const int MAXN = ;

 const int MAXM = ;

 const double eps = 1e-;

 int dx[] = { -, , , };

 int dy[] = {, -, , };

 double p;

 int  id[][], cnt, N = ;

 double G[][], a[][], ans[];

 struct NODE

 {

         int u, d;

         NODE(){}

         NODE(int _u, int _d):u(_u), d(_d){}

 };

 void preInit()

 {

         cnt = ;

         for(int i=;i<=N-;i++)

         {

                 for(int j = ; j <= i; j ++ )

                 {

                         id[i][j] = cnt++;

                 }

         }

         id[N][N-] = cnt ++;

 }

 void init()

 {

         preInit();

         mem0(G);mem0(a);

         for(int i=;i<=N-;i++)

         {

                 for(int j=;j<=i;j++)

                 {

                         int nx = MAX(i, j + ), ny = MIN(i, j + );

                         G[id[i][j]][id[nx][ny]] = p;

                         nx = i; ny = (j- ) >=  ? j- : ;

                         G[id[i][j]][id[nx][ny]] = -p;

                 }

         }

         for ( int i = ; i < cnt; i++ )

         {

                 a[i][i] = a[i][cnt] = 1.0;

                 if ( i == cnt- ) { a[i][cnt] = ; }

                 for ( int j = ; j < cnt; j++ ) if ( fabs ( G[i][j] ) > eps )

                             a[i][j] -= G[i][j];

         }

 }

 void gauss()

 {

         int col = ;

         for(int k=;k<cnt && col < cnt;k++, col ++)

         {

                     double Max = fabs(a[k][col]);

                     int  Maxr = k;

                     for(int r = k + ; r < cnt; r ++)

                             if( fabs(a[r][col])  -  Max  >  eps  )

                                     Max  =  fabs( a[Maxr = r][col] );

                     if(fabs(Max) < eps)  {  k --;  continue;  }

                     for(int c = col;c<=cnt; c ++ )

                             SWAP(a[Maxr][c],  a[k][c]  );

                     for(int r = k + ; r < cnt; r ++ ) if( fabs(a[r][col]) > eps  )

                     {

                             double tmp = a[r][col] / a[k][col];

                             for(int c = col; c <= cnt; c ++ )

                                     a[r][c] -= tmp * a[k][c];

                     }

         }

         for(int r=cnt-;r>=;r--)

         {

                     for(int c = cnt-;c>r;c--)

                             a[r][cnt] -= a[r][c] * ans[c];

                     ans[r] = a[r][cnt] / a[r][r];

         }

 }

 int main()

 {

 //        FOPENIN ( "in.txt" );

 //       FOPENOUT("out.txt");

         while ( ~scanf ( "%lf", &p ) )

         {

                 init();

                gauss();

                 printf("%.9lf\n", ans[]);

         }

         return ;

 }

HDU 4870 Rating（高斯消元 ）的相关教程结束。