以 杭电zhu and 772002 为例,看高斯消元如何使用
矩阵的秩等于矩阵线性无关的列的数量,自由度 = 矩阵列数 - 矩阵的秩。摘自贴吧。
线性代数,忘了。好像是通过行变换,得到一个上三角矩阵,求出某一变元的解,然后回迭求解。
下面是个例子:
#define _CRT_SECURE_ON_WARNINGS
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
#define MOD 1000000007
const int maxn = 2017; //素数最大不超过2000
typedef int Matrix[maxn][maxn];
int Rank(Matrix A, int m, int n) { //增广矩阵A,m行数(方程式个数),n列数。均从0开始
int i = 0, j = 0, k, r, u;
while (i < m && j < n) { //显然j和列有关
r = i; //选主行
for (k = i; k < m; ++k) {
if (A[k][j]) {
r = k;
break;
}
}
if (A[r][j]) {
if (r != i)
for (k = 0; k <= n; ++k)
swap(A[r][k], A[i][k]); //行交换,将第i行换为A[i][j]不为0
for (u = i + 1; u < m; ++u) //将第j列第i行一下的值全变为0。当然想改变A[u][j],u这一行都要做相同的改变(即全异或A[i][k])
if (A[u][j])
for (k = i; k <= n; ++k)
A[u][k] ^= A[i][k]; //若不是异或,普通的矩阵改为减除
++i; //移动主行
}
++j;//移动列
}
return i;//直接返回秩了,没有回带求解
}
Matrix A;
bool vis[maxn];
int prime[304]; //存素数
int getPrime(int n) { //O(n)线性素数筛
memset(vis, 1, sizeof vis);
int cnt = 0; //从0开始存,若换用注释中的,则从1开始存
for (int i = 2; i <= n; ++i) {
if (vis[i]) {
prime[cnt++] = i;//p[++cnt] = i;
}
for (int j = 0; j < cnt && (i * prime[j] <= n); ++j) { //j = 1; j <= cnt
vis[i * prime[j]] = 0;
if (i % prime[j] == 0) break;
}
}
return cnt; //返回2:n中素数的个数
}
long long pow_mod(int n) {
long long b = 2, res = 1;
while (n) {
if (n & 1) res = (res * b) % MOD;
b = (b * b) % MOD;
n >>= 1;
}
return res;
}
int main()
{
ios::sync_with_stdio(false);
int m = getPrime(2000); //303
int T, kase = 0;
cin >> T;
while (T--) {
int n, maxp = 0; //maxp 存最大素数的下标
long long x;
cin >> n;
memset(A, 0, sizeof A); //每一列表示一个数的的质因子分解
for (int i = 0; i < n; ++i) {
cin >> x;
for (int j = 0; j < m; ++j)
while (x % prime[j] == 0) {
maxp = max(maxp, j);
x /= prime[j];
A[j][i] ^= 1; //第i个元素是否含有奇数个质因子prime[j]
}
} //求出增广矩阵
int r = Rank(A, maxp + 1, n); //求秩
long long ans = pow_mod(n - r); //自由度 = 矩阵列数 - 矩阵的秩
ans = (ans - 1 + MOD) % MOD; //去掉全不选的一种情况
cout << "Case #" << ++kase << ":" << endl << ans << endl;
}
return 0;
}
POJ1830 套高斯消元的板子,加了一个判断无果的判定,会计算机实现高斯消元后,问题的关键就在于写出方程组,为何要求自由元。
#define _CRT_SECURE_ON_WARNINGS
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
#define MOD 1000000007
const int maxn = 50;
typedef int Matrix[maxn][maxn];
int Rank(Matrix A, int m, int n) {
int i = 0, j = 0, k, r, u;
while (i < m && j < n) { //第n列为矩阵右边的值
r = i; //选主行
for (k = i; k < m; ++k) {
if (A[k][j]) { //向下选,定主行
r = k;
break;
}
}
if (A[r][j]) {
if (r != i)//主行上移,i行才是真正的主子
for (k = 0; k <= n; ++k)
swap(A[r][k], A[i][k]);
for (u = i + 1; u < m; ++u)
if (A[u][j])
for (k = i; k <= n; ++k)
A[u][k] ^= A[i][k];
++i;
}
++j;
}
for (int k = i; k < m; ++k)
if (A[k][n] != 0) return -1;//左为0,等号右不为0
return i;
}
Matrix A;
long long pow_mod(int n) {
long long b = 2, res = 1;
while (n) {
if (n & 1) res = (res * b) % MOD;
b = (b * b) % MOD;
n >>= 1;
}
return res;
}
int main()
{
ios::sync_with_stdio(false);
int T, kase = 0;
cin >> T;
while (T--) {
int n;
cin >> n;
memset(A, 0, sizeof A);
for (int i = 0; i < n; ++i)
cin >> A[i][n];
int tmp;
for (int i = 0; i < n; ++i) {
cin >> tmp;
A[i][n] ^= tmp;
A[i][i] = 1;
}
int i, j;
while (cin >> i >> j) {
if (i == 0 && j == 0) break;
A[j - 1][i - 1] = 1; //第i个开关影响第j个
}
int r = Rank(A, n, n); //求秩
if (r == -1) {
puts("Oh,it's impossible~!!");
continue;
}
long long ans = pow_mod(n - r); //自由度 = 矩阵列数 - 矩阵的秩
cout << ans << endl;
}
return 0;
}