Gauss 消元（模板）

Posted 2020-08-08 hk lin

  1 /* 
  2     title：Gauss消元整数解/小数解整数矩阵模板 
  3     author：lhk
  4     time: 2016.9.11
  5     没学vim的菜鸡自己手打了     
  6 */
  7 #include<cstdio>
  8 #include<iostream>
  9 #include<cstring>
 10 #include<algorithm>
 11 #include<cmath>
 12 #define clr(x) memset(x,0,sizeof(x))
 13 #define clrdown(x) memset(x,-1,sizeof(x))
 14 #define maxn 910
 15 #define maxm 910 
 16 #define mod 3 
 17 using namespace std;
 18 int A[maxn][maxm];//Gauss消元的增广矩阵 
 19 int free_x[maxm];//自由变元的位置 
 20 int x[maxm];//整数解集 
 21 double xd[maxm];//小数解集 
 22 void init(int n,int m);//矩阵初始化操作 
 23 int gauss(int n,int m);//gauss消元部分 
 24 void print(int n);//输出解集 
 25 int exgcd(int a,int b,int &x,int &y);//扩展欧几里得求逆元,对于模mod的矩阵除法需要 
 26 int lcm(int a,int b);
 27 int gcd(int a,int b);
 28 int main()
 29 {
 30     int T,n,m; 
 31     scanf("%d",&T);
 32     while(T--)
 33     {
 34         scanf("%d%d",&n,&m);
 35         init(n,m);
 36         int p=gauss(n,m); 
 37         if(p==-1)
 38         {
 39             printf("no answer\n");
 40             return 0; 
 41          } 
 42         if(p==-2)
 43         { 
 44                printf("have float answer\n");
 45             return 0; 
 46         }
 47         printf("the num of free_x is %d\n",p);
 48         print(m);
 49     }
 50     return 0;
 51 }
 52 //输出解的和,自由变元数量以及各个解 
 53 void print(int n)
 54 {
 55     //    若有小数解换为xd输出 
 56     int sum=0;
 57     for(int i=0;i<n;i++)
 58         sum+=x[i];
 59     printf("sum=%d\n",sum); 
 60     for(int i=0;i<n;i++)
 61         printf("x%d=%d\n",i+1,x[i]);
 62     return ;
 63 }
 64 //读入增广矩阵 
 65 void init(int n,int m)
 66 {
 67     clr(A);
 68     for(int i=0;i<n;i++)
 69         for(int j=0;j<=m;j++)
 70         scanf("%d",&A[i][j]); 
 71     clrdown(x);
 72     clr(free_x);
 73     return ;
 74 }
 75 int gauss(int n,int m)//输出-1是无解，-2是有小数解，>=0则是自由变元的数量和全为整数解 
 76 {
 77     int k,col,num=0,max_r,dou,max_x,LCM,ta,tb;
 78 //k为当前操作行，col为操作主元素所在列 
 79     for(k=0,col=0;k<n && col<m;k++,col++)
 80     {
 81         //若A[K][col]不为col列最大，则将k行与k+1到n-1行中A[i][col]绝对值最大的行交换
 82         max_r=k;
 83         max_x=abs(A[k][col]);
 84         for(int i=k+1;i<n;i++)
 85             if(max_x<abs(A[i][col]))
 86             {
 87                 max_x=abs(A[i][col]);
 88                 max_r=i;
 89             }
 90         if(max_r!=k)
 91         {
 92             for(int j=col;j<=m;j++)
 93                 swap(A[k][j],A[max_r][j]);
 94         }
 95         //若k到n-1行A[i][col]全为0,则主元素指向当前行下一列的元素 
 96         if(A[k][col]==0)
 97         {
 98             k--;
 99             free_x[num++]=col;
100             //自由变元为当前col 
101             continue;
102         }
103         for(int i=k+1;i<n;i++)
104         if(A[i][col])
105         {
106             LCM=lcm(abs(A[k][col]),abs(A[i][col]));
107             ta=LCM/abs(A[i][col]);
108             tb=LCM/abs(A[k][col]);
109             if(A[k][col]*A[i][col]<0) tb=-tb;//若符号不同则取反 
110             for(int j=col;j<=m;j++)
111             {
112                     A[i][j]=A[i][j]*ta-A[k][j]*tb;
113                  // A[i][j]=((A[i][j]*ta-A[k][j]*tb)%mod+mod)%mod; //模mod矩阵的消元 
114             }
115         }
116     }
117     //k行及之后若有(0,0,0……,0,a)(a!=0）的行，则无解输出-1 
118     for(int i=k;i<n;i++)
119         if(A[i][col]!=0)
120             return -1;    
121     int temp;
122     //对自由变元的赋值，可有多种方式 
123     //若为小数解则换为xd; 
124     for(int i=0;i<num;i++)
125         x[free_x[i]]=0;
126     //int  xi,yi; exgcd需要 
127     for(int i=k-1,c=m-1;i>=0;c=m-1,i--)
128     {
129         temp=A[i][m];
130         while(x[c]!=-1)
131         {
132             if(A[i][c])
133                  temp-=x[c]*A[i][c]; 
134                 //temp=((temp-(x[c]*A[i][c])%mod)%mod+mod)%mod;//模mod矩阵的回代  
135             c--;
136         }
137         if(temp%A[i][c]!=0) return -2;
138         x[c]=temp/A[i][c]; 
139         /*exgcd(A[i][c],mod,xi,yi);
140         xi=(xi%mod+mod)%mod;
141         x[c]=(temp*xi%mod+mod)%mod;*/ //模mod 矩阵的x[i]的赋值 
142     }
143     return col-k;
144 }
145 int exgcd(int a,int b,int &x,int &y)
146 {
147     if(b==0)
148     {
149         x=1;
150         y=0;
151         return a;
152     }
153     else
154     {
155         int r=exgcd(b,a%b,y,x);
156         y-=x*(a/b);
157         return r;
158     }
159 }
160 int gcd(int a,int b)
161 {
162     int c;
163     while(b!=0)
164     {
165         c=a%b;
166         a=b;
167         b=c;
168     }
169     return a;
170 }
171 int lcm(int a,int b)
172 {
173     return a/gcd(a,b)*b;
174 }