《挑战程序设计竞赛》 4.1.1 矩阵 P286
想写几篇挑战的感悟,也有助于自己理解这本书。但这上面大多贴的是书上的代码,主要是为了用的时候后直接复制就好了,这样就很方便了,就相当于黑盒模板了。
1.线性方程组
/** \brief 高斯消元法
*
* \param 复杂度 n的3次方(n是方程数)
* \param
* \return 返回答案数组
*
*/
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\n", res[0]); }