高斯消元模版

  1. 云栖社区>
  2. 博客>
  3. 正文

高斯消元模版

angel_kitty 2017-07-05 10:00:00 浏览907
展开阅读全文

这模版敲了我俩个小时+写注释,参考自kuangbin!

两百行的大模拟,累死了QAQ

下面附上模版!

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

 

网友评论

登录后评论
0/500
评论
angel_kitty
+ 关注