POJ2947 Widget Factory 高斯消元

  题目连接:http://poj.org/problem?id=2947

  求解形如:

    a11×x11%p+a12×x12%p+......+a1n×x1n%p=k1%p

    a21×x21%p+a22×x22%p+......+a2n×x2n%p=k2%p

    ......

    an1×xn1%p+an2×xn2%p+......+ann×xnn%p=kn%p

   直接取模求解线性方程组。在网上找了一个整数Gauss的模板,貌似网上很多都是用的这个,但这个模板在多解情况下求变元是否确定有问题,应该用高斯-约当消元求出最简的上三角矩阵:   

     2   1  -1    8                   1   0   0   2
    -3  -1   2  -11   ——>      0   1   0   3 
    -2   1   2   -3                   0   0   1  -1
那么这样就可以在O(n)的时间内判断变元是否确定。
  1 //STATUS:C++_AC_2657MS_452KB
  2 #include <functional>
  3 #include <algorithm>
  4 #include <iostream>
  5 //#include <ext/rope>
  6 #include <fstream>
  7 #include <sstream>
  8 #include <iomanip>
  9 #include <numeric>
 10 #include <cstring>
 11 #include <cassert>
 12 #include <cstdio>
 13 #include <string>
 14 #include <vector>
 15 #include <bitset>
 16 #include <queue>
 17 #include <stack>
 18 #include <cmath>
 19 #include <ctime>
 20 #include <list>
 21 #include <set>
 22 #include <map>
 23 using namespace std;
 24 //using namespace __gnu_cxx;
 25 //define
 26 #define pii pair<int,int>
 27 #define mem(a,b) memset(a,b,sizeof(a))
 28 #define lson l,mid,rt<<1
 29 #define rson mid+1,r,rt<<1|1
 30 #define PI acos(-1.0)
 31 //typedef
 32 typedef long long LL;
 33 typedef unsigned long long ULL;
 34 //const
 35 const int N=310;
 36 const int INF=0x3f3f3f3f;
 37 const int MOD=100000,STA=8000010;
 38 const LL LNF=1LL<<60;
 39 const double EPS=1e-8;
 40 const double OO=1e15;
 41 const int dx[4]={-1,0,1,0};
 42 const int dy[4]={0,1,0,-1};
 43 const int day[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};
 44 //Daily Use ...
 45 inline int sign(double x){return (x>EPS)-(x<-EPS);}
 46 template<class T> T gcd(T a,T b){return b?gcd(b,a%b):a;}
 47 template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
 48 template<class T> inline T lcm(T a,T b,T d){return a/d*b;}
 49 template<class T> inline T Min(T a,T b){return a<b?a:b;}
 50 template<class T> inline T Max(T a,T b){return a>b?a:b;}
 51 template<class T> inline T Min(T a,T b,T c){return min(min(a, b),c);}
 52 template<class T> inline T Max(T a,T b,T c){return max(max(a, b),c);}
 53 template<class T> inline T Min(T a,T b,T c,T d){return min(min(a, b),min(c,d));}
 54 template<class T> inline T Max(T a,T b,T c,T d){return max(max(a, b),max(c,d));}
 55 //End
 56 
 57 int a[N][N];//增广矩阵
 58 int x[N];//解集
 59 bool free_x[N];//标记是否是不确定的变元
 60 int n,m,k;
 61 
 62 //
 63 void Debug(int equ,int var)
 64 {
 65     int i, j;
 66     for (i = 0; i < equ; i++)
 67     {
 68         for (j = 0; j < var + 1; j++)
 69         {
 70             cout << a[i][j] << " ";
 71         }
 72         cout << endl;
 73     }
 74     cout << endl;
 75 }
 76 //
 77 
 78 inline int gcd(int a,int b)
 79 {
 80     int t;
 81     while(b!=0)
 82     {
 83         t=b;
 84         b=a%b;
 85         a=t;
 86     }
 87     return a;
 88 }
 89 inline int lcm(int a,int b)
 90 {
 91     return a/gcd(a,b)*b;//先除后乘防溢出
 92 }
 93 
 94 // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解,
 95 //-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数)
 96 //有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
 97 int Gauss(int equ,int var)
 98 {
 99     int i,j,k;
100     int max_r; // 当前这列绝对值最大的行.
101     int col; //当前处理的列
102     int ta,tb;
103     int LCM;
104     int temp;
105     int free_x_num;
106     int free_index;
107 
108     for(int i=0;i<=var;i++)
109     {
110         x[i]=0;
111         free_x[i]=true;
112     }
113 
114     //转换为阶梯阵.
115     col=0; // 当前处理的列
116     for(k = 0;k < equ && col < var;k++,col++)
117     {// 枚举当前处理的行.
118 // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
119         max_r=k;
120         for(i=k+1;i<equ;i++)
121         {
122             if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
123         }
124         if(max_r!=k)
125         {// 与第k行交换.
126             for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
127         }
128         if(a[k][col]==0)
129         {// 说明该col列第k行以下全是0了,则处理当前行的下一列.
130             k--;
131             continue;
132         }
133         for(i=0;i<equ;i++)
134         {// 枚举要删去的行.
135             if(a[i][col]!=0 && i!=k)
136             {
137                 LCM = lcm(abs(a[i][col]),abs(a[k][col]));
138                 ta = LCM/abs(a[i][col]);
139                 tb = LCM/abs(a[k][col]);
140                 if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
141                 for(j=0;j<var+1;j++)
142                 {
143                     a[i][j] = ((a[i][j]*ta-a[k][j]*tb)%7+7)%7;
144                 }
145             }
146         }
147     }
148 
149   //  Debug(equ,var);
150 
151     // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
152     for (i = k; i < equ; i++)
153     {
154         if (a[i][col] != 0) return -1;
155     }
156     // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换.
157     // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵.
158     // 且出现的行数即为自由变元的个数.
159     if (k < var)return 1;
160  /*   {
161         // 首先,自由变元有var - k个,即不确定的变元至少有var - k个.
162         for (i = k - 1; i >= 0; i--)
163         {
164             // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行.
165             // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的.
166             free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它们仍然为不确定的变元.
167             for (j = 0; j < var; j++)
168             {
169                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
170             }
171             if (free_x_num > 1) continue; // 无法求解出确定的变元.
172             // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的.
173             temp = a[i][var];
174             for (j = 0; j < var; j++)
175             {
176                 if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
177             }
178             x[free_index] = temp / a[i][free_index]; // 求出该变元.
179             free_x[free_index] = 0; // 该变元是确定的.
180         }
181         return var - k; // 自由变元有var - k个.
182     }*/
183     // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
184     // 计算出Xn-1, Xn-2 ... X0.
185     /*
186     for (i = var - 1; i >= 0; i--)
187     {
188         temp = a[i][var];
189         for (j = i + 1; j < var; j++)
190         {
191             if (a[i][j] != 0) temp -= a[i][j] * x[j];
192         }
193         if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解.
194         x[i] = temp / a[i][i];
195     }*/
196     return 0;
197 }
198 
199 map<string,int> q;
200 
201 int main()
202 {
203 //    freopen("in.txt","r",stdin);
204     int i,j,num,ans;
205     char s[10],e[10];
206     int cnt[N];
207     q["MON"]=1,q["TUE"]=2,q["WED"]=3,
208     q["THU"]=4,q["FRI"]=5,q["SAT"]=6,q["SUN"]=7;
209     while(~scanf("%d%d",&n,&m) && (n || m))
210     {
211         for(i=0;i<m;i++){
212             mem(cnt,0);
213             scanf("%d%s%s",&k,s,e);
214             while(k--){
215                 scanf("%d",&num);
216                 cnt[num-1]++;
217             }
218             for(j=0;j<n;j++)
219                 a[i][j]=cnt[j]%7;
220             if(q[s]<=q[e])a[i][n]=(q[e]-q[s]+1)%7;
221             else a[i][n]=(8-q[s]+q[e])%7;
222         }
223      //   Debug(m,n);
224 
225         ans=Gauss(m,n);
226 
227         if(ans==0){
228             for(i=0;i<n;i++){
229                 for(j=3;j<=9;j++)
230                     if(a[i][i]*j%7==a[i][n]){
231                         x[i]=j;
232                         break;
233                     }
234             }
235             printf("%d",x[0]);
236             for(i=1;i<n;i++)
237                 printf(" %d",x[i]);
238         }
239         else if(ans==-1)printf("Inconsistent data.");
240         else printf("Multiple solutions.");
241         putchar('\n');
242     }
243     return 0;
244 }
原文地址:https://www.cnblogs.com/zhsl/p/3109531.html