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 }
时间： 2024-10-22 23:30:28