高斯消元模板(+解异或方程组)
2015-08-03 10:42
92 查看
原题链接是hiho第56周:http://hihocoder.com/contest/hiho56/problem/1
用来练手高斯消元,该模板只能判断无解,多解(不能判断有几个自由变元)和小数解的情况(可以用强制转换变成整数)
高斯消元的除了解普通的方程组外,还有一种经典的问题---解异或方程组
题目链接:http://hihocoder.com/contest/hiho57/problem/1
经典问题题意:有一个矩阵放置的灯泡,状态分别为打开或者关闭。打开或者关闭当前灯泡会影响它的上下左右相邻的灯泡的状态(打开变关闭,关闭变打开)。问有没有方法使所有灯泡全开(或全关?),并输出应该操作哪些灯泡。
做法:具体的做法在hiho的链接里面有讲,不再细讲,仅贴出代码
网上有另一种模板比较详细,但是没用过,
用来练手高斯消元,该模板只能判断无解,多解(不能判断有几个自由变元)和小数解的情况(可以用强制转换变成整数)
//其中a为方程组等号左边的矩阵,b为原方程组右边的值,temp为辅助数组,val为答案,代表等号左边xi的值 double a[maxn * 2][maxn];//等号左边的矩阵 double b[maxn * 2];//等号右边的值 double temp[maxn * 2]; double val[maxn * 2]; void swap_line(int i,int j)//交换一整行,等号左边右边和原val都要交换 { memcpy(temp,a[j],sizeof(temp)); memcpy(a[j],a[i],sizeof(temp)); memcpy(a[i],temp,sizeof(temp)); swap(b[i],b[j]); } int n,m; //这里的n为等号左边矩阵的列数,如果等号左边也在a中则n要+1,把b数组去掉,放到a数组的最后一个元素后面,列数为1~n //m为行数,行数为1~m int Gauss()
{//处理上三角矩阵 for(int i = 1;i <= n;i++) { int flag = 0; for(int j = i;j <= m;j++)//从第i行开始,找到第i列不等于0的行j if(fabs(a[j][i]) > eps) { swap_line(j,i); flag = 1; break; } if(!flag)//若无法找到,则存在多个解 return 2;//解多余一个 //消除第i+1行到第m行的第i列 for(int j = i + 1;j <= m;j++) temp[j] = a[j][i] / a[i][i]; for(int j = i + 1;j <= m;j++) { for(int k = 1;k <= n;k++) a[j][k] = a[j][k] - a[i][k] * temp[j]; b[j] = b[j] - b[i] * temp[j]; } } //检查是否无解,即存在0 = x的情况 for(int i = 1;i <= m;i++) { if(fabs(b[i]) < eps)continue; int all_zero = 1;//假设第i列系数均为0 for(int j = 1;j <= n;j++) if(fabs(a[i][j]) > eps) { all_zero = 0;break; } if(all_zero) return 0;//无解 } //此时存在唯一解 //由于每一行都比前一行少一个系数,所以在m行中只有前n行有系数 //从第n行开始处理每一行的解 for(int i = n;i >= 1;i--) {//利用已经计算出的结果,将第i行中第i+1列至第n列的系数消除 for(int j = i + 1;j <= n;j++){ b[i] = b[i] - a[i][j] * val[j]; a[i][j] = 0; } val[i] = b[i] / a[i][i]; } return 1; //唯一解 //最后的答案存放在val数中,注意解可能是小数 }
高斯消元的除了解普通的方程组外,还有一种经典的问题---解异或方程组
题目链接:http://hihocoder.com/contest/hiho57/problem/1
经典问题题意:有一个矩阵放置的灯泡,状态分别为打开或者关闭。打开或者关闭当前灯泡会影响它的上下左右相邻的灯泡的状态(打开变关闭,关闭变打开)。问有没有方法使所有灯泡全开(或全关?),并输出应该操作哪些灯泡。
做法:具体的做法在hiho的链接里面有讲,不再细讲,仅贴出代码
#include<bits/stdc++.h> #define mem(a,b) memset(a,b,sizeof(a)) #define FOR(i,a,b) for(int i = a;i <= b;i++) using namespace std; typedef long long ll; const int maxn = 55; const int INF = 0x3f3f3f3f; const double eps = 1e-7; //其中a为方程组等号左边的矩阵,b为原方程组右边的值,temp为辅助数组,val为答案,代表等号左边xi的值 //其中异或方程组无解的情况是出现000000000... 最后一个是1的行,那么原异或方程组无解 int a[maxn][maxn];//等号左边的矩阵 int b[maxn];//等号右边的值 int temp[maxn]; double val[maxn]; int n,m; //这里的n为等号左边矩阵的列数,如果等号左边也在a中则n要+1,把b数组去掉,放到a数组的最后一个元素后面,列数为1~n //m为行数,行数为1~m int Gauss() {//处理上三角矩阵 for(int i = 1;i <= n;i++) { int flag = 0; int jilu; for(int j = i;j <= m;j++)//从第i行开始,找到第i列不等于0的行j if(a[j][i]) { for(int t = 1;t <= n;t++) swap(a[i][t],a[j][t]); swap(b[i],b[j]); jilu = i; flag = 1; break; } // if(!flag)//若无法找到,则存在多个解 // continue; for(int j = 1;j <= m;j++) if(a[j][i] == 1 && j != jilu) { for(int k = i;k <= n;k++) a[j][k] = a[j][k] ^ a[i][k]; b[j] = b[j] ^ b[i]; } } } char ma[7][8]; int dx[] = {0,0,1,-1}; int dy[] = {1,-1,0,0}; int main() { string s; for(int i = 1;i <= 5;i++) for(int j = 1;j <= 6;j++) scanf(" %c",&ma[i][j]); mem(a,0); for(int i = 1;i <= 5;i++) for(int j = 1;j <= 6;j++) { int now = (i - 1) * 6 + j; b[now] = (1 ^ (ma[i][j] - '0')); a[now][now] = 1; for(int k = 0;k < 4;k++) { int x = i + dx[k],y = j + dy[k]; if(x < 1 || x > 5 || y < 1 || y > 6)continue; int haha = (x - 1) * 6 + y; a[now][haha] = 1; } } n = 30,m = 30; Gauss(); int num = 0; for(int i = 1;i <= 30;i++) if(b[i]) num++; printf("%d\n",num); for(int i = 1;i <= 30;i++) { if(b[i]) { int x = (i - 1) / 6 + 1; int y = (i % 6 == 0) ? 6 : i % 6; printf("%d %d\n",x,y); } } }
网上有另一种模板比较详细,但是没用过,
#include<stdio.h> #include<algorithm> #include<iostream> #include<string.h> #include<math.h> using namespace std; const int MAXN=50; int a[MAXN][MAXN];//增广矩阵 int x[MAXN];//解集 bool free_x[MAXN];//标记是否是不确定的变元 /* void Debug(void) { int i, j; for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { cout << a[i][j] << " "; } cout << endl; } cout << endl; } */ inline int gcd(int a,int b) { int t; while(b!=0) { t=b; b=a%b; a=t; } return a; } inline int lcm(int a,int b) { return a/gcd(a,b)*b;//先除后乘防溢出 } // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解, //-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数) //有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var. int Gauss(int equ,int var) { int i,j,k; int max_r;// 当前这列绝对值最大的行. int col;//当前处理的列 int ta,tb; int LCM; int temp; int free_x_num; int free_index; for(int i=0;i<=var;i++) { x[i]=0; free_x[i]=true; } //转换为阶梯阵. col=0; // 当前处理的列 for(k = 0;k < equ && col < var;k++,col++) {// 枚举当前处理的行. // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差) max_r=k; for(i=k+1;i<equ;i++) { if(abs(a[i][col])>abs(a[max_r][col])) max_r=i; } if(max_r!=k) {// 与第k行交换. for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]); } if(a[k][col]==0) {// 说明该col列第k行以下全是0了,则处理当前行的下一列. k--; continue; } for(i=k+1;i<equ;i++) {// 枚举要删去的行. if(a[i][col]!=0) { LCM = lcm(abs(a[i][col]),abs(a[k][col])); ta = LCM/abs(a[i][col]); tb = LCM/abs(a[k][col]); if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加 for(j=col;j<var+1;j++) { a[i][j] = a[i][j]*ta-a[k][j]*tb; } } } } // Debug(); // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0). for (i = k; i < equ; i++) { // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换. if (a[i][col] != 0) return -1; } // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵. // 且出现的行数即为自由变元的个数. if (k < var) { // 首先,自由变元有var - k个,即不确定的变元至少有var - k个. for (i = k - 1; i >= 0; i--) { // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行. // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的. free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它们仍然为不确定的变元. for (j = 0; j < var; j++) { if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j; } if (free_x_num > 1) continue; // 无法求解出确定的变元. // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的. temp = a[i][var]; for (j = 0; j < var; j++) { if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j]; } x[free_index] = temp / a[i][free_index]; // 求出该变元. free_x[free_index] = 0; // 该变元是确定的. } return var - k; // 自由变元有var - k个. } // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵. // 计算出Xn-1, Xn-2 ... X0. for (i = var - 1; i >= 0; i--) { temp = a[i][var]; for (j = i + 1; j < var; j++) { if (a[i][j] != 0) temp -= a[i][j] * x[j]; } if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解. x[i] = temp / a[i][i]; } return 0; } int main(void) { freopen("in.txt", "r", stdin); freopen("out.txt","w",stdout); int i, j; int equ,var; while (scanf("%d %d", &equ, &var) != EOF) { memset(a, 0, sizeof(a)); for (i = 0; i < equ; i++) { for (j = 0; j < var + 1; j++) { scanf("%d", &a[i][j]); } } // Debug(); int free_num = Gauss(equ,var); if (free_num == -1) printf("无解!\n"); else if (free_num == -2) printf("有浮点数解,无整数解!\n"); else if (free_num > 0) { printf("无穷多解! 自由变元个数为%d\n", free_num); for (i = 0; i < var; i++) { if (free_x[i]) printf("x%d 是不确定的\n", i + 1); else printf("x%d: %d\n", i + 1, x[i]); } } else { for (i = 0; i < var; i++) { printf("x%d: %d\n", i + 1, x[i]); } } printf("\n"); } return 0; }
相关文章推荐
- C语言编程入门——链表
- Android App 启动过程分析
- IOS软键盘弹出和隐藏
- [尺取法] poj3061 Subsequence
- css3 box-sizing属性
- eclipse中新建分支
- 判断设备是否安装QQ
- ASP.NET定时执行任务 (不使用外接程序,.net内部机制实现)
- 【LeetCode】242 Valid Anagram
- zoj 1610 Count the Colors
- iOS输入框禁止输入emoji表情
- 数据结构和算法目录表
- 南邮 OJ 1025 完数
- DI和IOC
- new 与malloc 的区别
- cocos2dx3.0rc导出自定义类到lua的方法详细步骤
- 深度学习——人工神经网络再掀研究热潮
- 【转】(索引)java sql包下面的类的用法
- hihoCoder 1078 线段树的区间修改
- 第11章 散列表