同余意义下的高斯消元法
以下高斯消元指的是列主元高斯消元法。
高斯消元法除了解浮点数线性方程,它还可以用于求矩阵的秩,求某些线性空间中的一组线性无关的基,解同余方程组。
求秩
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