《挑战程序设计竞赛》 4.1.1 矩阵 P286

　　想写几篇挑战的感悟，也有助于自己理解这本书。但这上面大多贴的是书上的代码，主要是为了用的时候后直接复制就好了，这样就很方便了，就相当于黑盒模板了。

 　　1.线性方程组

/** rief 高斯消元法
 *
 * param 复杂度 n的3次方（n是方程数）
 * param 
 * 
eturn 返回答案数组
 *
 */     

const double EPS= 1e-8 ;

typedef vector<double> vec;
typedef vector<vec> mat;

//求解Ax=b，其中A是方阵
//当方程组无解或有无穷多解时，返回一个长度为0的数组
vec gauss_jordan(const mat& A,const vec& b){
    int n=A.size();
    mat B(n,vec(n+1));
    //把A复制给B
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++) B[i][j]=A[i][j];
    //把正在处理的未知数系数的绝对值最大的式子换到第i行
    for (int i=0;i<n;i++) B[i][n]=b[i];
    for (int i=0;i<n;i++){
        int pivot=i;
        for (int j=i;j<n;j++){
            if (abs(B[j][i])>abs(B[pivot][i])) pivot=j;
        }
        swap(B[i],B[pivot]);

        //无解或有无穷多解
        if (abs(B[i][i])<EPS) return vec();

        //把正在处理的未知数系数变为1
        for (int j=i+1;j<=n;j++) B[i][j]/=B[i][i];
        for (int j=0;j<n;j++){
            if (i!=j){
                //从第j个式子中消去第i个未知数
                for (int k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k];
            }
        }
    }
    vec x(n);
    //存放在右边的b就是答案
    for (int i=0;i<n;i++) x[i]=B[i][n];
    return x;
}

 　　2.期望值和方程组

　　看到这的时候，忘了期望怎么求了，我就去查了下。简单来说就是这样的：

期望是E(X)，且E(X) = X1*p(X1) + X2*p(X2) + …… + Xn*p(Xn)

X1,X2,X3,……,Xn为这几个数据，p(X1),p(X2),p(X3),……p(Xn)为这几个数据的概率函数。

书上说的最后的那个 +1 ，还是有些不懂，学长说是1步，有点道理。

const double EPS= 1e-9 ;

typedef vector<double> vec;
typedef vector<vec> mat;

//求解Ax=b，其中A是方阵
//当方程组无解或有无穷多解时，返回一个长度为0的数组
vec gauss_jordan(const mat& A,const vec& b){
    int n=A.size();
    mat B(n,vec(n+1));
    //把A复制给B
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++) B[i][j]=A[i][j];
    //把正在处理的未知数系数的绝对值最大的式子换到第i行
    for (int i=0;i<n;i++) B[i][n]=b[i];
    for (int i=0;i<n;i++){
        int pivot=i;
        for (int j=i;j<n;j++){
            if (abs(B[j][i])>abs(B[pivot][i])) pivot=j;
        }
        swap(B[i],B[pivot]);

        //无解或有无穷多解
        if (abs(B[i][i])<EPS) return vec();

        //把正在处理的未知数系数变为1
        for (int j=i+1;j<=n;j++) B[i][j]/=B[i][i];
        for (int j=0;j<n;j++){
            if (i!=j){
                //从第j个式子中消去第i个未知数
                for (int k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k];
            }
        }
    }
    vec x(n);
    //存放在右边的b就是答案
    for (int i=0;i<n;i++) x[i]=B[i][n];
    return x;
}


const int MAXN= 10 + 5 ;
const int MAXM= 10 + 5 ;

//输入
char grid[MAXN][MAXM];
int N,M;

bool can_goal[MAXN][MAXM]; //can_goal[x][y]为true的话，(x,y)可以到达终点
int dx[]={-1,1,0,0},dy[]={0,0,-1,1};

//搜索课到达终点的点
void dfs(int x,int y){
    can_goal[x][y]=true;
    for (int i=0;i<4;i++){
        int nx=x+dx[i],ny=y+dy[i];
        if (0<=nx&&nx<N&&0<=ny&&ny<M&&!can_goal[nx][ny]&&
            grid[nx][ny]!='#'){
            dfs(nx,ny);
        }
    }
}

//--------------------------------分割线-----------------------------------------

void solve(){
    mat A(N*M,vec(N*M,0));//初始化，全为0
    vec b(N*M,0);
    for (int x=0;x<N;x++)
        for (int y=0;y<M;y++)
            can_goal[x][y]=false;

    dfs(N-1,M-1);

    //构建矩阵
    for(int x=0;x<N;x++){
        for (int y=0;y<M;y++){

            //到达终点，或是(x,y)无法到达终点的情况
            if (x==N-1&&y==M-1||!can_goal[x][y]){
                A[x*M+y][x*M+y]=1;
                continue;
            }

            //其余情况
            int move=0;
            for (int k=0;k<4;k++){
                int nx=x+dx[k],ny=y+dy[k];
                if (0<=nx&&nx<N&&0<=ny&&ny<M&&grid[nx][ny]!='.'){
                    A[x*M+y][nx*M+ny]=-1;
                    move++;
                }
            }
            b[x*M+y]=A[x*M+y][x*M+y]=move;
        }
    }
    vec res=gauss_jordan(A,b);
    printf("%.8f
", res[0]);
}