同余意义下的高斯消元法
以下高斯消元指的是列主元高斯消元法。
高斯消元法除了解浮点数线性方程，它还可以用于求矩阵的秩，求某些线性空间中的一组线性无关的基，解同余方程组。
求秩
SGU 200 Cracking RSA
题意： 给你m个正的整型数，这m个数的质因子仅来自前t个质数，求它们能组成多少完全平方数。（每个数最多选一次，且不能一个数都不选）
Sample test(s)
Input
(t=)3 (m=)4
9 20 500 3
Output
3
Sample 解释： 9 = 3 ^ 2, 20 * 500 = 100 ^ 2 , 9 * 20 * 500 = 300 ^ 2是完全平方数，其余均不是。
思路：
首先可以现将它们分解，事实上，3 ^ 1 和 3 ^ 3从某种意义是一样的，因为完全平方数要求所有的不同种类质因子的个数是偶数。于是就可以用二元域来表示。
如上面的例子
2 3 5
9 0 0 0
20 0 0 1
500 0 0 1
几个数能组成完成平方数就等价于它在这样的矩阵表示中的行向量相加（二元域的加法，即异或运算）等于一个 0 向量。用高斯消元的思想化解为标准阶梯型矩阵，并且求出矩阵的秩r，并且有一定有r<=t，并且原矩阵和当前的标准阶梯型矩阵是等价的，也就是所原先m个数，与现在经过高斯消元后的新的数是完全等价的。由于每一个数，只能选取一次，而如果选择了行向量表示中非零向量的元素又行向量已经是线性无关的，则一定不能构成零向量，即完全平方数。所以要组成完成平方数就只能在行向量表示中为零向量的元素选，这样的元素总有（m-r）个，每一元素选择选或者不选，但不能全不选，所以答案就是2 ^ （m-r ） -1 。
举上面的例题为例，消元之后变为
2 3 5
20 0 0 1
9 0 0 0
1000 0 0 0
这个矩阵的秩为1，行向量表示中为零向量的元素总有（3-1）个，所以答案为2 ^（3-1）-1。
当然这道题的 m - r 最大可能等于100，即所给的数都是完全平方数，这时已经2 ^ （m-r ） -1 已经超出了long long的范围，需要手写大整数，至于求 2 ^ （m-r ） -1，再使用快速幂运算即可求出。
#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <stdlib.h>
#include <stack>
#include <vector>
#include <string.h>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;
const int MAXN=1300;
bool Isn_Prime[MAXN+5];
int prime[230];
void GetPrime(void)
{
ms(Isn_Prime);
prime[0]=0;
for(int i=2;i<=MAXN;i++)
{
if(!Isn_Prime[i]) prime[++prime[0]]=i;
for(int j=1;j<=prime[0]&&prime[j]<=MAXN/i;j++)
{
Isn_Prime[prime[j]*i]=true;
if(i%prime[j]==0) break;
}
}
}
int m,t;
int a[105][105];
int Guass(void)
{
int max_r,col,k;
for(k=col=0;k<m&&col<t;k++,col++)
{
max_r=k;
for(int i=k+1;i<m;i++)
if(abs(a[i][col])>abs(a[max_r][col]))
max_r=i;
if(a[max_r][col]==0){
k--;continue;
}
if(max_r!=k){
for(int j=col;j<t;j++)
swap(a[k][j],a[max_r][j]);
}
for(int i=k+1;i<m;i++)
if(a[i][col]){
for(int j=col;j<t;j++)
a[i][j]^=a[k][j];
}
}
return m-k;
}
struct BigInt
{
const static int mod=10000;
const static int DLEN=4;
int a[600],len;
BigInt(void){
ms(a);
len=1;
}
BigInt(int v){
ms(a);
len=0;
do{
a[len++]=v%mod;
v/=mod;
}while(v);
}
BigInt operator +(const BigInt &b)const {
BigInt res;
res.len=max(len,b.len);
for(int i=0;i<=res.len;i++)
res.a[i]=0;
for(int i=0;i<res.len;i++)
{
res.a[i]+=((i<len)?a[i]:0)+((i<b.len)?b.a[i]:0);
res.a[i+1]+=res.a[i]/mod;
res.a[i]%=mod;
}
if(res.a[res.len]>0) res.len++;
return res;
}
BigInt operator *(const BigInt &b)const {
BigInt res;
for(int i=0;i<len;i++)
{
int up=0;
for(int j=0;j<b.len;j++)
{
int tmp=a[i]*b.a[j]+res.a[i+j]+up;
res.a[i+j]=tmp%mod;
up=tmp/mod;
}
if(up)
res.a[i+b.len]=up;
}
res.len=len+b.len;
while(res.a[res.len-1]==0&&res.len>1)
res.len--;
return res;
}
void output(void)
{
printf("%d",a[len-1] );
for(int i=len-2;i>=0;i--)
printf("%04d",a[i] );
putchar('\n');
}
};
BigInt quick_pow(int k)
{
BigInt a=BigInt(2),
res=BigInt(1);
while(k){
if(k&1) res=res*a;
k>>=1;
a=a*a;
}
return res+BigInt(-1);
}
int main(int argc, char const *argv[])
{
GetPrime();
while(scanf("%d %d",&t,&m)==2){
ms(a);
for(int i=0;i<m;i++)
{
int tmp;
scanf("%d",&tmp);
for(int j=0;j<t;j++)
while(tmp%prime[j+1]==0)
{
tmp/=prime[j+1];
a[i][j]^=1;
}
}
BigInt rs=quick_pow(Guass());
rs.output();
}
return 0;
}求基
Xor
比如这么一道题，求一组最小的基
Xor
题意：xor运算符是一种二进制运算，对于两个整数a,b，如果c=a xor b，那么c二进制表示的每一位的值由a二进制表示和b二进制表示的相应位进行不进位加法得到，例如1 xor 1=0,1 xor 0 =1。
现在给出 n个数，那么用xor运算符可以得到很多不同的值，现在问第几小的数？
Example
2
1 2
那么由1,2这两个数进行xor运算可以得到3个数 1,2,3.
对于这道题，我们可以将这些数用二进制写出来
2 ^ 0 2 ^ 1 2 ^ 2
1 1 0 0
2 0 1 0
3 1 1 0
4 0 0 1
显然3可以由1 xor 2得到，也就是1,2,3他们是线性相关的，
所以我们应该对这n个数张成的空间求出一组线性无关的基，又由于求的是第几小，所以希望每个数消元后的数都尽量小，所以应该从高位开始消元 进行消元。
当然，这里的消元用的是高斯消元的那种思想，没必要把整个矩阵弄出来。
如果这道题求的是第k大，其实可以转化为第几小，应该如果求出来它的秩是n，那么总共可以组成2^n-1个数，那么其实就是第2^n-k小。
#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <stdlib.h>
#include <stack>
#include <vector>
#include <string.h>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;
LL a[10005];
int n;
//求第几小应该从最高位开始消元，这样保证了线性无关组元素和最小
void Guass(void)
{
int max_r,col,k;
for(k=0,col=60;k<n&&col>=0;col--)
{
max_r=k;
for(int i=k+1;i<n;i++)
if(a[i]&(1ll<<col)) {
max_r=i;
break;
}
if((a[max_r]&(1ll<<col))==0)
continue;
if(max_r!=k)
swap(a[k],a[max_r]);
for(int i=0;i<n;i++)
if(i!=k&&(a[i]&(1ll<<col)))
a[i]^=a[k];
k++;
}
sort(a,a+n);
n=unique(a,a+n)-a;
}
LL cal(LL k)
{
LL ans=0;
int i=0;
if(a[0]==0){
if(k==1) return 0;
k--;
i++;
}
for(;i<n&&k;k>>=1,i++)
if(k&1)
ans^=a[i];
if(i==n&&k) return -1;
return ans;
}
int main(int argc, char const *argv[])
{
int t,ti=0;
scanf("%d",&t);
while(++ti<=t){
scanf("%d",&n);
for(int i=0;i<n;i++)
scanf("%I64d",a+i);
Guass();
printf("Case #%d:\n",ti);
int q;
scanf("%d",&q);
while(q--){
LL k;
scanf("%I64d",&k);
printf("%I64d\n",cal(k) );
}
}
return 0;
}一类开关问题之解中最少的1
最简单的开关问题就是问能否达到某种状态，建立对2取模的01同余方程组，高斯消元即可。
下面考虑其扩展——求出解中最小的1的一组解
The Water Bowls
The Water Bowls
例如这样一道题，
Description
The cows have a line of 20 water bowls from which they drink. The bowls can be either right-side-up (properly oriented to serve refreshing cool water) or upside-down (a position which holds no water). They want all 20 water bowls to be right-side-up and thus use their wide snouts to flip bowls.
Their snouts, though, are so wide that they flip not only one bowl but also the bowls on either side of that bowl (a total of three or – in the case of either end bowl – two bowls).
Given the initial state of the bowls (1=undrinkable, 0=drinkable – it even looks like a bowl), what is the minimum number of bowl flips necessary to turn all the bowls right-side-up?
Input
Line 1: A single line with 20 space-separated integers
Output
Line 1: The minimum number of bowl flips necessary to flip all the bowls right-side-up (i.e., to 0). For the inputs given, it will always be possible to find some combination of flips that will manipulate the bowls to 20 0’s.
Sample Input
0 0 1 1 1 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0
Sample Output
3
Hint
Explanation of the sample:
Flip bowls 4, 9, and 11 to make them all drinkable:
0 0 1 1 1 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 [initial state]
0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 [after flipping bowl 4]
0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 [after flipping bowl 9]
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [after flipping bowl 11]
设第i个bowl是否翻转为xi，xi为真表示翻转它一次，假表示不动。
原题是否某种状态转化为零向量，事实上等价于由零向量转化到该状态（y1 , y2 , …… , y20）
所以方程就是
x1 + x2 = y1
x1 + x2 + x3 = y2
x2 + x3 + x4 = y3
……
x18 + x19 + x20 = y19
x19 + x20 = y20
在将它建立方程组后，
它就转化为求解一个对2取模的01方程组，且找解中1的个数最少的一个。
在高斯消元的过程中将自由变元的个数k求出，再进行枚举自由变元（枚举 2 ^ k 次）求解高斯消元后的方程组，然后记录最少1的个数即可。
#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <stdlib.h>
#include <stack>
#include <vector>
#include <string.h>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;
const int MAXN=25;
int g[MAXN][MAXN];
int x[MAXN];
int free_x[MAXN];
int free_num;
int Guass(void)
{
int max_r,col,k;
free_num=0;
for(k=col=0;k<20&&col<20;k++,col++)
{
max_r=k;
for(int i=k+1;i<20;i++)
if(abs(g[i][col])>abs(g[max_r][col]))
max_r=i;
if(g[max_r][col]==0){
k--;
free_x[free_num++]=col;
continue;
}
if(max_r!=k){
for(int j=col;j<21;j++)
swap(g[k][j],g[max_r][j]);
}
for(int i=k+1;i<20;i++)
if(g[i][col]){
for(int j=col;j<21;j++)
g[i][j]^=g[k][j];
}
}
for(int j=col;j<20;j++)
free_x[free_num++]=j;
for(int i=k;i<20;i++)
if(g[i][20]) return -1;
if(k<20) return 20-k;
for(int i=19;i>=0;i--)
{
x[i]=g[i][20];
for(int j=i+1;j<20;j++)
x[i]^=(g[i][j]&&x[j]);
}
return 0;
}
void solve(void)
{
int t=Guass();
if(t==0){
int ans=0;
for(int i=0;i<20;i++)
ans+=x[i];
printf("%d\n",ans );
}
else {
int ans=0x3f3f3f3f;
int tot=(1<<t);
for(int i=0;i<tot;i++)
{
int cnt=0;
for(int j=0;j<t;j++)
if(i&(1<<j)){
x[free_x[j]]=1;
cnt++;
}
else x[free_x[j]]=0;
for(int j=19-t;j>=0;j--)
{
int idx;
for(idx=j;idx<20;idx++)
if(g[j][idx]) break;
x[idx]=g[j][20];
for(int l=idx+1;l<20;l++)
if(g[j][l])
x[idx]^=x[l];
cnt+=x[idx];
}
ans=min(cnt,ans);
}
printf("%d\n",ans );
}
}
int main(int argc, char const *argv[])
{
ms(g);
for(int i=0;i<20;i++)
{
g[i][i]=1;
if(i>0) g[i-1][i]=1;//注意不要写反
if(i<19) g[i+1][i]=1;
}
for(int i=0;i<20;i++)
scanf("%d",&g[i][20]);
solve();
return 0;
}解同余方程组
这里仅讨论模数为质数，即保证一定存在逆元。
Widget Factory
Widget Factory
与实数的高斯消元法存在一点区别，就是除法变成了乘法逆元。
由于求乘法逆元比较慢，所以在比如
a00 * x0 + a01 * x1 = c0 （*）
a10 * x0 + a11 * x1 = c1 ( ** )的消去中，采取
( * )不变，（ ** ） 变为
( * ) a00 - ( * ) *a11。
但如果直接这样乘，当p很大时，式子的系数可能会溢出，所以要分别乘以lcm(a00,a11) / a11 , lcm(a00,a11) / a00
在回代的过程中，以得到的最后一个方程组为例 a * x = cc
那么x = cc * a的逆元。
其余与实数域的高斯消元原理上没有区别。
#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <stdlib.h>
#include <stack>
#include <vector>
#include <string.h>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;
const int MOD = 7;
const int MAXN=400;
int a[MAXN][MAXN];
int x[MAXN];
inline int gcd(int a,int b)
{return b?gcd(b,a%b):a;}
inline int lcm(int a,int b)
{return a/gcd(a,b)*b;}
LL inv(LL a,LL m)
{return a==1?1:(inv(m%a,m)*(m-m/a)%m);}
int Guass(int equ,int var)
{
int max_r,col,k;
for(k=col=0;k<equ&&col<var;k++,col++)
{
max_r=k;
for(int i=k+1;i<equ;i++)
if(abs(a[i][col])>abs(a[max_r][col]))
max_r=i;
if(a[max_r][col]==0){
k--;
continue;
}
if(max_r!=k)
for(int j=col;j<var+1;j++)
swap(a[k][j],a[max_r][j]);
for(int i=k+1;i<equ;i++)
{
if(a[i][col]){
int LCM=lcm(abs(a[i][col]),abs(a[k][col]));
int ta=LCM/abs(a[i][col]);
int tb=LCM/abs(a[k][col]);
if(a[i][col]*a[k][col]<0) tb=-tb;
for(int j=col;j<var+1;j++)
a[i][j]=((a[i][j]*ta-a[k][j]*tb)
%MOD+MOD)%MOD;
}
}
}
for(int i=k;i<equ;i++)
if(a[i][var]) return -1;
if(k<var) return var-k;
for(int i=var-1;i>=0;i--)
{
int tmp=a[i][var];
for(int j=i+1;j<var;j++)
if(a[i][j]){
tmp-=a[i][j]*x[j];
tmp=(tmp%MOD+MOD)%MOD;
}
x[i]=(tmp*inv(a[i][i],MOD))%MOD;
}
return 0;
}
int change(char *s)
{
if(strcmp(s,"MON")==0) return 1;
else if(strcmp(s,"TUE")==0) return 2;
else if(strcmp(s,"WED")==0) return 3;
else if(strcmp(s,"THU")==0) return 4;
else if(strcmp(s,"FRI")==0) return 5;
else if(strcmp(s,"SAT")==0) return 6;
else return 7;
}
int main(int argc, char const *argv[])
{
int n,m;
while(scanf("%d %d",&n,&m)==2&&(n||m)){
ms(a);
char str1[6],str2[6];
int k;
for(int i=0;i<m;i++)
{
scanf("%d %s %s",&k,str1,str2);
a[i][n]=((change(str2)-change(str1)+1)
%MOD+MOD)%MOD;
while(k--){
int t;
scanf("%d",&t);
t--;
a[i][t]++;
a[i][t]%=MOD;
}
}
int ans=Guass(m,n);
if(ans==0){
for(int i=0;i<n;i++)
printf("%d%c",x[i]>2?x[i]:(x[i]+7),
i==n-1?'\n':' ' );
}
else if(ans==-1) puts("Inconsistent data.");
else puts("Multiple solutions.");
}
return 0;
}误差：讨论的高斯消元是在整数p-域上的，不存在精度问题，所以就不会出现误差。
北京2014现场赛H题
Happy Matt Friends
Happy Matt Friends
题意：给出n个数，和一个数m，问从这n个数选取若干个数（每个数最多选一次，可以全不选）能通过xor运算得到的数 >=m 的方案数。
例子
n=3 m=2
1 2 3
那么
1 ^ 2 =3 >=2
2 >= 2
3 >= 2
3 ^ 1 >=2
总共四种方案。
当然，由于这道题的数据范围比较小，可以直接用动态规划的思想解决。但当数据范围达到10 ^ 18 时，就只能用高斯消元来解决，复杂度为O(N*log(max(M,ki)))。
用上面Xor的思想来解决，先高斯消元求出秩r和N-r，这N个数所能组成的数就是2 ^ r,而每个数可以有2 ^ ( N - r )种，所以只需要求出M在这 2 ^ r 个数排第几（二分查找）。
#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <stdlib.h>
#include <stack>
#include <vector>
#include <string.h>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;
const int MAXN=50;
int a[MAXN],n,nen;
int Guass(void)
{
int i,k,col,max_r;
for(k=0,col=20;k<n&&col>=0;col--)
{
max_r=k;
for(i=k+1;i<n;i++)
if(a[i]&(1<<col)){
max_r=i;
break;
}
if(a[max_r]&(1<<col)){
if(k!=max_r)
swap(a[k],a[max_r]);
for(i=0;i<n;i++)
if(i!=k&&(a[i]&(1<<col))) a[i]^=a[k];
k++;
}
else continue;
}
a[n]=0;
sort(a,a+n+1);
nen=unique(a,a+n+1)-a;
nen--;
return n-k;//相当于变元个数
}
LL cal(LL num)
{
LL res=0ll;
if(num==1) return 0;
num--;
for(int i=1;i<=nen&#num>>=1,i++)
if(num&1) res^=a[i];
return res;
}
int main(int argc, char const *argv[])
{
int t,ti=0;
scanf("%d",&t);
while(++ti<=t){
int m;
scanf("%d%d",&n,&m);
for(int i=0;i<n;i++)
scanf("%d",a+i);
int var=Guass();
LL AllN=1ll<<nen;
LL l=1,r=AllN;
while(l<r){
LL mid=l+(r-l)/2;
if(cal(mid)<m) l=mid+1;
else if(cal(mid)>m) r=mid;
else if(cal(mid)==m){
r=mid;
break;
}
}
printf("Case #%d: ",ti );
if(cal(r)>=m)
printf("%I64d\n",(AllN-r+1)*(1ll<<var) );
else puts("0");
}
return 0;
}
/*
3
3 4
1 2 3
3 3
1 4 8
3 4
1 2 8
1
8 30021
11683 21469 15437 24343 23066 6108 1892 24276
1
2 26072
20538 23197
*/
1949
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
