高斯消元模板

  1 #include<stdio.h>
  2 #include<algorithm>
  3 #include<iostream>
  4 #include<string.h>
  5 #include<math.h>
  6 using namespace std;
  7 
  8 const int MAXN=50;
  9 
 10 
 11 
 12 int a[MAXN][MAXN];//增广矩阵
 13 int x[MAXN];//解集
 14 bool free_x[MAXN];//标记是否是不确定的变元
 15 
 16 
 17 
 18 /*
 19 void Debug(void)
 20 {
 21     int i, j;
 22     for (i = 0; i < equ; i++)
 23     {
 24         for (j = 0; j < var + 1; j++)
 25         {
 26             cout << a[i][j] << " ";
 27         }
 28         cout << endl;
 29     }
 30     cout << endl;
 31 }
 32 */
 33 
 34 
 35 inline int gcd(int a,int b)
 36 {
 37     int t;
 38     while(b!=0)
 39     {
 40         t=b;
 41         b=a%b;
 42         a=t;
 43     }
 44     return a;
 45 }
 46 inline int lcm(int a,int b)
 47 {
 48     return a/gcd(a,b)*b;//先除后乘防溢出
 49 }
 50 
 51 // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解,
 52 //-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数)
 53 //有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
 54 int Gauss(int equ,int var)
 55 {
 56     int i,j,k;
 57     int max_r;// 当前这列绝对值最大的行.
 58     int col;//当前处理的列
 59     int ta,tb;
 60     int LCM;
 61     int temp;
 62     int free_x_num;
 63     int free_index;
 64 
 65     for(int i=0;i<=var;i++)
 66     {
 67         x[i]=0;
 68         free_x[i]=true;
 69     }
 70 
 71     //转换为阶梯阵.
 72     col=0; // 当前处理的列
 73     for(k = 0;k < equ && col < var;k++,col++)
 74     {// 枚举当前处理的行.
 75 // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
 76         max_r=k;
 77         for(i=k+1;i<equ;i++)
 78         {
 79             if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
 80         }
 81         if(max_r!=k)
 82         {// 与第k行交换.
 83             for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
 84         }
 85         if(a[k][col]==0)
 86         {// 说明该col列第k行以下全是0了,则处理当前行的下一列.
 87             k--;
 88             continue;
 89         }
 90         for(i=k+1;i<equ;i++)
 91         {// 枚举要删去的行.
 92             if(a[i][col]!=0)
 93             {
 94                 LCM = lcm(abs(a[i][col]),abs(a[k][col]));
 95                 ta = LCM/abs(a[i][col]);
 96                 tb = LCM/abs(a[k][col]);
 97                 if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
 98                 for(j=col;j<var+1;j++)
 99                 {
100                     a[i][j] = a[i][j]*ta-a[k][j]*tb;
101                 }
102             }
103         }
104     }
105 
106   //  Debug();
107 
108     // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
109     for (i = k; i < equ; i++)
110     { // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换.
111         if (a[i][col] != 0) return -1;
112     }
113     // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵.
114     // 且出现的行数即为自由变元的个数.
115     if (k < var)
116     {
117         // 首先,自由变元有var - k个,即不确定的变元至少有var - k个.
118         for (i = k - 1; i >= 0; i--)
119         {
120             // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行.
121             // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的.
122             free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它们仍然为不确定的变元.
123             for (j = 0; j < var; j++)
124             {
125                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
126             }
127             if (free_x_num > 1) continue; // 无法求解出确定的变元.
128             // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的.
129             temp = a[i][var];
130             for (j = 0; j < var; j++)
131             {
132                 if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
133             }
134             x[free_index] = temp / a[i][free_index]; // 求出该变元.
135             free_x[free_index] = 0; // 该变元是确定的.
136         }
137         return var - k; // 自由变元有var - k个.
138     }
139     // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
140     // 计算出Xn-1, Xn-2 ... X0.
141     for (i = var - 1; i >= 0; i--)
142     {
143         temp = a[i][var];
144         for (j = i + 1; j < var; j++)
145         {
146             if (a[i][j] != 0) temp -= a[i][j] * x[j];
147         }
148         if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解.
149         x[i] = temp / a[i][i];
150     }
151     return 0;
152 }
153 int main(void)
154 {
155     freopen("in.txt", "r", stdin);
156     freopen("out.txt","w",stdout);
157     int i, j;
158     int equ,var;
159     while (scanf("%d %d", &equ, &var) != EOF)
160     {
161         memset(a, 0, sizeof(a));
162         for (i = 0; i < equ; i++)
163         {
164             for (j = 0; j < var + 1; j++)
165             {
166                 scanf("%d", &a[i][j]);
167             }
168         }
169 //        Debug();
170         int free_num = Gauss(equ,var);
171         if (free_num == -1) printf("无解!
");
172    else if (free_num == -2) printf("有浮点数解,无整数解!
");
173         else if (free_num > 0)
174         {
175             printf("无穷多解! 自由变元个数为%d
", free_num);
176             for (i = 0; i < var; i++)
177             {
178                 if (free_x[i]) printf("x%d 是不确定的
", i + 1);
179                 else printf("x%d: %d
", i + 1, x[i]);
180             }
181         }
182         else
183         {
184             for (i = 0; i < var; i++)
185             {
186                 printf("x%d: %d
", i + 1, x[i]);
187             }
188         }
189         printf("
");
190     }
191     return 0;
192 }
原文地址:https://www.cnblogs.com/zxz666/p/10726077.html