zoj3772 Calculate the Function 矩阵除法模板

原创于 2014-04-08 12:54:57 发布 · 759 阅读

1 ·

CC 4.0 BY-SA版权

ACM 同时被 2 个专栏收录

20 篇文章

订阅专栏

矩阵

0 篇文章

订阅专栏

本文详细介绍了如何使用线性代数原理和C++实现矩阵的逆运算，进而解决线性方程组的问题。通过矩阵乘法、矩阵逆运算的模板函数，实现了高效地求解任意阶方阵的逆矩阵，并应用到实际问题中。

比赛的时候用线段树水过的，然后本以为求矩阵逆元只能在double精度下进行，后来才发现图样了。

研究一番线代，终于将可逆矩阵（准确来说是方阵）的除法模板给写出来了，该方法理论上可以适用于任意N阶的矩阵求模（虽然我小题大做的用在了这题上= =）。

注意：如果对应矩阵不可逆则会异常退出，所以使用前先判断矩阵是否可逆！！

具体细节看代码。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<string>
#include<algorithm>
using namespace std;
#define MOD 1000000007
#define MAXN 2
/********矩阵乘法定义*********/
long long exgcd(long long a, long long b, long long &x, long long &y) {
	if (b == 0) {
		x = 1;
		y = 0;
		return a;
	}
	long long r = exgcd(b, a % b, y, x);
	y -= x * (a / b);
	return r;
}
long long Inv(long long a) { //若MOD是素数且gcd(a,MOD)=1可利用费马小定理直接求a^(p-2)%MOD
	long long r, x, y;
	r = exgcd(a, MOD, x, y);
	if (r == 1)
		return (x % MOD + MOD) % MOD;
	return -1;
}
struct Matrix {
	long long M[MAXN][MAXN];
	void clear() {
		memset(M, 0, sizeof(M));
	}
	void getone() {
		for (int i = 0; i < MAXN; i++) {
			for (int j = 0; j < MAXN; j++) {
				if (i == j)
					M[i][j] = 1;
				else
					M[i][j] = 0;
			}
		}
	}
	void out() {
		cout << '#' << endl;
		for (int i = 0; i < MAXN; i++) {
			for (int j = 0; j < MAXN; j++) {
				cout << M[i][j];
				if (j == MAXN - 1)
					cout << endl;
				else
					cout << ' ';
			}
		}
	}
	Matrix operator *(Matrix othr) const {
		Matrix ans;
		ans.clear();
		for (int i = 0; i < MAXN; i++) {
			for (int j = 0; j < MAXN; j++) {
				for (int k = 0; k < MAXN; k++) {
					ans.M[i][j] += M[i][k] * othr.M[k][j] % MOD;
					ans.M[i][j] %= MOD;
				}
			}
		}
		return ans;
	}
	Matrix operator *(long long x) const {
		Matrix ans;
		for (int i = 0; i < MAXN; i++) {
			for (int j = 0; j < MAXN; j++) {
				ans.M[i][j] = M[i][j] * x % MOD;
			}
		}
		return ans;
	}
	Matrix operator +(Matrix othr) const {
		Matrix ans;
		for (int i = 0; i < MAXN; i++) {
			for (int j = 0; j < MAXN; j++) {
				ans.M[i][j] = M[i][j] + othr.M[i][j];
				ans.M[i][j] %= MOD;
			}
		}
		return ans;
	}
	Matrix operator /(Matrix othr) const {
		return this->operator *(othr.Minv());
	}
	Matrix Minv() {//矩阵求逆
		long long IS[MAXN], JS[MAXN];
		long long det;
		Matrix ans;
		for (int i = 0; i < MAXN; i++)
			for (int j = 0; j < MAXN; j++)
				ans.M[i][j] = M[i][j];
		for (int k = 0; k < MAXN; k++){
			det = 0;
			for (int i = k; i < MAXN; i++)
				for (int j = k; j < MAXN; j++)
					if (abs(ans.M[i][j]) > det) {
						IS[k] = i;
						JS[k] = j;
						det = abs(ans.M[i][j]);
					}
			if (det == 0)exit(1); //异常
				if (IS[k] != k)
					for (int i = 0; i < MAXN; i++)
						swap(ans.M[k][i], ans.M[IS[k]][i]);
			if (JS[k] != k)
				for (int i = 0; i < MAXN; i++)
					swap(ans.M[i][k], ans.M[i][JS[k]]);
			ans.M[k][k] = Inv(ans.M[k][k]);
			for (int i = 0; i < MAXN; i++)
				if (i != k)
					ans.M[k][i] = ans.M[k][i] * ans.M[k][k] % MOD;
			for (int i = 0; i < MAXN; i++)
				if (i != k)
					for (int j = 0; j < MAXN; j++)
						if (j != k)
							ans.M[i][j] = ((ans.M[i][j]- ans.M[i][k] * ans.M[k][j])%MOD + MOD) % MOD;
			for (int i = 0; i < MAXN; i++)
				if (i != k)ans.M[i][k] = MOD - ans.M[i][k] * ans.M[k][k] % MOD;

		}
		for (int k = MAXN - 1; k >= 0; k--) {
			if (JS[k] != k)
				for (int i = 0; i < MAXN; i++)
					swap(ans.M[k][i], ans.M[JS[k]][i]);
			if (IS[k] != k)
				for (int i = 0; i < MAXN; i++)
					swap(ans.M[i][k], ans.M[i][IS[k]]);
		}
		return ans;
	}
};

Matrix a[100005];
int b[100005];
int main() {
	int T;
	//freopen("D://a.in","r",stdin);
	scanf("%d", &T);
	while (T--) {
		int n, m;
		scanf("%d%d", &n, &m);
		for (int i = 1; i <= n; i++) {
			scanf("%d", &b[i]);
			a[i].M[0][0] = 1;
			a[i].M[1][0] = 1;
			a[i].M[0][1] = b[i];
			a[i].M[1][1] = 0;
			if (i > 1)
				a[i] = a[i] * a[i - 1];
			//out(a[i]);
		}
		while (m--) {
			int l, r;
			scanf("%d%d", &l, &r);
			if (r - l < 2) {
				printf("%d\n", b[r]);
			} else {
				Matrix ans = a[r] / a[l + 1];
				long long res = ans.M[0][0] * b[l + 1] % MOD
						+ ans.M[0][1] * b[l] % MOD;
				printf("%lld\n", res % MOD);
			}
		}
	}
}