luoguP4783 [模板]矩阵求逆 线性代数

\(n^2\)的矩阵的逆


翻了翻题解,看到了初等矩阵这个东西,突然想起来在看线代的时候看到过....

然后又温习了一遍线性代数的知识

不妨设\(PA = E\),其中\(P\)是一堆初等矩阵的积(必须同时是行变换)

由于\(PA = E, PE = P\),因此\(P(A, E) = (E, P)\)

所以我们只要对矩阵\((A, E)\)来做初等变换

由于我们只做行变换

因此,两个分块矩阵之间互相不干扰

所以当左侧的\(A\)变化为\(E\)时,右边的\(E\)自然变成了\(P\)

复杂度\(O(n^3)\)

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

#define ri register int
#define rep(io, st, ed) for(ri io = st; io <= ed; io ++)
#define drep(io, ed, st) for(ri io = ed; io >= st; io --)

#define gc getchar
inline int read() { 
	int p = 0, w = 1; char c = gc();
	while(c > '9' || c < '0') { if(c == '-') w = -1; c = gc(); }
	while(c >= '0' && c <= '9') p = p * 10 + c - '0', c = gc();
	return p * w;
}

const int sid = 405;
const int mod = 1e9 + 7;

inline void inc(int &a, int b) { a += b; if(a >= mod) a -= mod; }
inline void dec(int &a, int b) { a -= b; if(a < 0) a += mod; }
inline int mul(int a, int b) { return 1ll * a * b % mod; }
inline int inv(int a) {
	int ret = 1, k = mod - 2;
	for( ; k; k >>= 1, a = mul(a, a))
		if(k & 1) ret = mul(ret, a);
	return ret;
}

int n;
int A[sid][sid], B[sid][sid];

inline int Guass() {
	rep(i, 1, n) {
		int pos = i; 
		rep(j, i + 1, n) if(A[j][i]) pos = j;
		if(!A[pos][i]) return 0;
		swap(A[i], A[pos]); swap(B[i], B[pos]);
		int IA = inv(A[i][i]);
		rep(j, 1, n) {
			if(i == j) continue;
			int ia = mul(A[j][i], IA);
			rep(k, 1, n) {
				if(k >= i) dec(A[j][k], mul(ia, A[i][k]));
				dec(B[j][k], mul(ia, B[i][k]));
			}
		}
	}
	rep(i, 1, n) {
		int IA = inv(A[i][i]);
		rep(j, 1, n) B[i][j] = mul(B[i][j], IA);
	}
	return 1;
}

int main() {
	n = read();
	rep(i, 1, n) rep(j, 1, n) A[i][j] = read();
	rep(i, 1, n) B[i][i] = 1;
	if(Guass()) {
		rep(i, 1, n) {
			rep(j, 1, n) printf("%d ", B[i][j]);
			printf("\n");
		}
	}
	else printf("No Solution\n");
	return 0;
}

也许下次我们可以出一道求\(AP = B\)或者\(PA = B\)\(P\)

相信能卡死一片人QAQ

posted @ 2018-12-17 21:36  remoon  阅读(229)  评论(0编辑  收藏  举报