double a[maxn][maxn] , x[maxn] ; //a[i][j] 系数矩阵 ， a[i][n+1] = y[i] ， x解 int n ; //n个方程 void guass(){ int i , j , k ; double sum , rate ; for(k = 1 ; k < n ; k++){ for(i = j = k ; i <= n ; i++) j = fabs(a[i][k]) > fabs(a[j][k]) ? i : j ; for(i = k ; i <= n+1 ; i++) swap(a[j][i] , a[k][i]) ; for(i = k+1 ; i <= n ; i++) for(rate = a[i][k]/a[k][k] , j = k ; j <= n+1 ; j++) a[i][j] -= a[k][j]*rate ; } for(i = n ; i >= 1 ; i--){ for(sum = 0 , j = i+1 ; j <= n ; sum += a[i][j]*x[j],j++) ; x[i] = (a[i][n+1] - sum)/a[i][i] ; } } double Fix(double x){ // 将 -0.000 变成 0.00 ,结果输出Fix(x[]) if(fabs(x) < 1e-4) return 0 ; return x ; }
另一个版本的高斯消元
double a[Max_N][Max_N], x[Max_N],b[Max_N][Max_N]; // 系数矩阵，double型 void guass(int n){ int i,j,k;// 这里 i表示行，j表示列，k表示某行的某列（确定一个元素） for(i=0; i<n; i++) for(j=0; j<n; j++) b[i][j]=a[i][j]; for(i=0; i<n; i++) b[i][n] = x[i]; for(j=0; j<n; j++){ for(k=i=j; i<n; i++) k = fabs(b[i][j]) > fabs(b[k][j]) ? i : k ; for(i=0; i<=n; i++) swap(b[j][i],b[k][i]); //把正在处理的那一行的对角线上的系数变成1，同一行除b[j][j],该位置的系数不处理 for(i=j+1 ; i<=n; i++) b[j][i] /= b[j][j]; for(i=0; i<n; i++){ if(i != j) for(k=j+1; k<=n; k++) b[i][k] -= b[i][j]*b[j][k]; // 做行变换，被消的同一列系数不处理 } } }
zoj 3645 BiliBili
题目来源： http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=4835
题目意思很明显
给你一个方程组 让你求（x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11）
但是方程组的形式是一个二次方程组
(ai1-x1)^2 + (ai2-x2)^2 +(ai3-x1)^2 + (ai4-x2)^2 +(ai5-x1)^2 + (ai6-x2)^2 +(ai7-x1)^2 + (ai8-x2)^2 + (ai9-x2)^2 +(ai10-x1)^2 + (ai11-x2)^2 = dis ^2
于是我们可以发现 用两个方程相减就可以得到一个一次方程组。 这样操作11次得到一个新的方程组，然后用高斯消元法求解即可。
代码如下：
using namespace std ; const double pi= acos(-1); #define arc(x) (x / 180 * pi) typedef long long LL ; const double EPS = 1e-8; const int Max_N = 15; double a[Max_N][Max_N], d[Max_N],b[Max_N][Max_N]; // 系数矩阵，double型 int n; // 方程的个数 void guass(){ int i,j,k; for(j=0; j<n; j++){ for(k=i=j; i<n; i++) k = fabs(b[i][j]) > fabs(b[k][j]) ? i : k ; for(i=0; i<=n; i++) swap(b[j][i],b[k][i]); for(i=j+1 ; i<=n; i++) b[j][i] /= b[j][j]; for(i=0; i<n; i++){ if(i != j) for(k=j+1; k<=n; k++) b[i][k] -= b[i][j]*b[j][k]; } } } int main(){ int t; scanf("%d",&t); while(t--){ n=11; for(int i=0; i<=n; i++){ for(int j=0; j<n; j++) scanf("%lf",&a[i][j]); scanf("%lf",&d[i]); } memset(b,0,sizeof(b)); for(int i=1; i<=n; i++){ for(int j=0; j<n; j++){ b[i-1][j] = 2*(a[i][j] - a[i-1][j] ); b[i-1][n] += a[i][j]*a[i][j] - a[i-1][j]*a[i-1][j]; } } for(int i=0; i<n; i++) b[i][n] += d[i]*d[i] - d[i+1]*d[i+1]; guass(); for(int i=0; i<n-1; i++) printf("%.2lf ",b[i][n] + EPS); // 输出时注意这里的判断， 处理-0.0000 printf("%.2lf\n",b[n-1][n] + EPS); } return 0; }
poj 1222 EXTENDED LIGHTS OUT
题目来源： http://poj.org/problem?id=1222
分析：转化为横坐标为灯（灯0，灯1，。。。，灯29），纵坐标为开关（开关0，开关1，。。。，开关29），记为矩阵A，矩阵A与灯的初始状态无关，A（i,j)=1表示开关j对灯i有影响。设矩阵X，X是30*1的矩阵，表示开关是否需要打开。X(i，0)=1表示开关i需要打开。矩阵X是待求的。设矩阵B 30*1，表示结果。B只与灯的初始状态有关，初始时灯i是开着的，则B（i,0)=1。解方程AX=B。
代码如下：
const int Max_N = 35; int a[Max_N][Max_N]; // 系数矩阵，double型，AX = Y int n; // 方程的个数 a[i][n]存储的是Y向量 void guass(){ int i,j,k; // i表示行，j表示列 for(i=0; i<n; i++){ for(k=i; k<n; k++) //寻找不为零的 a[i][i]，通过交换 if(a[k][i]) break; for(j=0; j<=n; j++) swap(a[k][j], a[i][j]); for(j=0; j<n; j++) if(i != j && a[j][i]) // 如果a[j][i] == 0 ，则该行不需要变换 for(k=0; k<=n; k++) a[j][k] ^= a[i][k]; } } int main(){ int t,u=1; scanf("%d",&t); while(t--){ int i,j; n=30; memset(a,0,sizeof(a)); for(i=0; i<n; i++) scanf("%d",&a[i][n]); // Y value for(i=0; i<n; i++){ a[i][i]=1; if(i>5) a[i][i-6] = 1 ;// up a[i][j] 表示第i个灯被第j个开关影响 if(i<24) a[i][i+6] =1 ; // down if(i%6) a[i][i-1] = 1 ; // left if(i%6 != 5) a[i][i+1] = 1 ; //right } guass(); printf("PUZZLE #%d\n",u++); for(i=0; i<n; i++){ if(i%6 != 5) printf("%d ",a[i][n]); else printf("%d\n",a[i][n]); } } return 0; }
