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     else154     {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 }