P4783 【模板】矩阵求逆

【模板】矩阵求逆

Luogu P4783

题目描述

求一个 \(N\times N\) 的矩阵的逆矩阵。答案对 \({10}^9+7\) 取模。

输入格式

第一行有一个整数 \(N\)，代表矩阵的大小；

接下来 \(N\) 行，每行 \(N\) 个整数，其中第 \(i\) 行第 \(j\) 列的数代表矩阵中的元素 \(a_{i j}\)。

输出格式

若矩阵可逆，则输出 \(N\) 行，每行 \(N\) 个整数，其中第 \(i\) 行第 \(j\) 列的数代表逆矩阵中的元素 \(b_{i j}\)，答案对 \({10}^9+7\) 取模；

否则只输出一行 No Solution。

样例 #1

样例输入 #1

3
1 2 8
2 5 6
5 1 2

样例输出 #1

718750005 718750005 968750007
171875001 671875005 296875002
117187501 867187506 429687503

样例 #2

样例输入 #2

3
3 2 4
7 2 9
2 4 3

样例输出 #2

No Solution

提示

对 \(30 \%\) 的数据有 \(N\le 100\)；
对 \(100 \%\) 的数据有 \(N\le 400\)，所有 \(0 \le a_{i j} < {10}^9 + 7\)。

Solution

假设有一个矩阵 \(A\)，如果想要计算除法，会发现我们没有定义矩阵的除法，不过会联想到数的除法，所以就想到了逆元这种方法。

矩阵的逆 \(A^{-1}\) 定义为满足 \(A^{-1}\times A=A\times A^{-1}=I\) 的矩阵（\(I\) 表示单位矩阵，即主对角线上的值都为 \(1\) 的矩阵）。并不是所有矩阵都存在逆矩阵。

假设矩阵 \(A\) 能够通过进行一系列的矩阵乘法来变成单位矩阵 \(I\) ，那么这过程中的所有矩阵的积就是需要找到的逆元。假设这过程中用到的矩阵为 \(p_1,p_2,p_3\cdots,p_k\)，那么就有 $$p_k\times p_{k-1}\times \cdots \times p_2 \times p_1\times A=I$$

也就是说需要找到合适的矩阵数列 \(p\) 来求矩阵逆。需要将矩阵逐步变成一个主对角线为 \(1\) 的单位矩阵，那么可以想到用高斯消元法求解，并且可以发现高斯消元法所用的初等行变换是可以用矩阵乘法进行表示的：

将第二行乘上一个数 \(c\)：

\[\begin{bmatrix} 1&0&0&\cdots&0\\ 0&c&0&\cdots&0\\ 0&0&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&1 \end{bmatrix} \]

倍加：

\[\begin{bmatrix} 1&0&0&\cdots&0\\ 0&1&0&\cdots&0\\ 0&c&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&1 \end{bmatrix} \]

所以就有一种做法：将原矩阵 \(A\) 和一个单位矩阵 \(I\) 左右拼接组成一个新的 \(n\times 2n\) 的矩阵，然后对左侧的矩阵 \(A\) 进行高斯消元，最后右边的单位矩阵 \(I\) 自然就乘上了数列 \(p\)，所以如果将初始矩阵表示为 \((A,I)\)，那么最终矩阵就是 \((I,A^{-1})\) 。

因为是模 \(10^9+7\)，所以除法需要用逆元，逆元用费马小定理搭配快速幂即可。整个代码的总时间复杂度为 \(\mathcal O(n^3\log p)\)，其中 \(n^3\) 来自高斯消元，\(\log p\) 来自快速幂求逆元。

Code

#include<bits/stdc++.h>
#define mem(a,b) memset(a,b,sizeof a)
#define int long long
using namespace std;
template<typename T> void read(T &k)
{
	k=0;T flag=1;char b=getchar();
	while (!isdigit(b)) {flag=(b=='-')?-1:1;b=getchar();}
	while (isdigit(b)) {k=k*10+b-48;b=getchar();}
	k*=flag;
}
const int MOD=1e9+7;
int Fpow(int a,int b)
{
	int base=a%MOD,res=1;
	while (b)
	{
		if (b&1) res=res*base%MOD;
		base=base*base%MOD,b>>=1;
	}
	return res;
}
const int _SIZE=4e2;
int a[_SIZE+5][(_SIZE<<1)+5];
int n;
void Calc()
{
	for (int i=1;i<=n;i++)
	{
		int r=i;
		for (int j=i+1;j<=n;j++)
			if (a[j][i]>a[r][i]) r=j;
		if (i!=r) swap(a[i],a[r]);
		if (!a[i][i]) {puts("No Solution");return;}
		int inv=Fpow(a[i][i],MOD-2);
		for (int k=1;k<=n;k++)
		{
			if (i==k) continue;
			int p=a[k][i]*inv%MOD;
			for (int j=i;j<=(n<<1);j++)
				a[k][j]=((a[k][j]-a[i][j]*p)%MOD+MOD)%MOD;	
		}
		for (int j=1;j<=(n<<1);j++)
			a[i][j]=(a[i][j]*inv)%MOD;
	}
	for (int i=1;i<=n;i++) {for (int j=n+1;j<=(n<<1);j++) printf("%lld ",a[i][j]); puts("");}
}
signed main()
{
	read(n);
	for (int i=1;i<=n;i++)
		for (int j=1;j<=n;j++) read(a[i][j]),a[i][i+n]=1;
	Calc();
	return 0;
}