高斯消元模板(转)
参考资料:http://blog.sina.com.cn/s/blog_afe6bb330101a59d.html
模板:
const int maxn=30; int a[maxn][maxn+1],x[maxn];//a是系数矩阵和增广矩阵,x是最后存放的解 // a[][maxn]中存放的是方程右面的值(bi) int equ,var;//equ是系数阵的行数,var是系数矩阵的列数(变量的个数) int free_num,ans=100000000; 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; } void swap(int &a,int &b) { int temp=a; //交换2个数 a=b; b=temp; } int Gauss() { int k,col = 0; //当前处理的列 for(k = 0; k < equ && col < var; ++k,++col) { int max_r = k; for(int i = k+1; i < equ; ++i) if(a[i][col] > a[max_r][col]) max_r = i; if(max_r != k) { for(int i = k; i < var + 1; ++i) swap(a[k][i],a[max_r][i]); } if(a[k][col] == 0) { k--; continue; } for(int i = k+1; i < equ; ++i) { if(a[i][col] != 0) { int LCM = lcm(a[i][col],a[k][col]); int ta = LCM/a[i][col], tb = LCM/a[k][col]; if(a[i][col]*a[k][col] < 0) tb = -tb; for(int j = col; j < var + 1; ++j) a[i][j] = ( (a[i][j]*ta)%2 - (a[k][j]*tb)%2 + 2 ) % 2; //a[i][j]只有0和1两种状态 } } } //上述代码是消元的过程,行消元完成 //解下来2行,判断是否无解 //注意K的值,k代表系数矩阵值都为0的那些行的第1行 for(int i = k; i < equ; ++i) if(a[i][col] != 0) return -1; // 无解返回 -1 //Debug(); //唯一解或者无穷解,k<=var //var-k==0 唯一解;var-k>0 无穷多解,自由解的个数=var-k //能执行到这,说明肯定有解了,无非是1个和无穷解的问题。 //下面这几行很重要,保证秩内每行主元非0,且按对角线顺序排列,就是检查列 for(int i = 0; i <equ; i++) if(!a[i][i]) { int j; for(j=i+1;j<=var;j++) if(a[i][j]) break; if(j == var) break; for(int k = 0; k < equ; ++k) swap(a[k][i],a[k][j]); } // ----处理保证对角线主元非0且顺序,检查列完成 // free_num=k; if (var-k>0) { //无穷多解,先枚举解,然后用下面的回带代码进行回带; //这里省略了下面的回带的代码;不管唯一解和无穷解都可以回带,只不过无穷解//回带时,默认为最后几个自由变元=0而已。 } if (var-k==0)//唯一解时 { //下面是回带求解代码,当无穷多解时,最后几行为0的解默认为0; for(int i = k-1; i >= 0; --i) //从消完元矩阵的主对角线非0的最后1行,开始往//回带 { int tmp = a[i][var] % 2; for(int j = i+1; j < var; ++j) //x[i]取决于x[i+1]--x[var]啊,所以后面的解对前面的解有影响。 if(a[i][j] != 0) tmp = ( tmp - (a[i][j]*x[j])%2 + 2 ) % 2; //if (a[i][i]==0) x[i]=tmp; //最后的空行时,即无穷解得 //else x[i] = (tmp/a[i][i]) % 2; //上面的正常解 } //回带结束了 } }