比赛的时候用线段树水过的,然后本以为求矩阵逆元只能在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);
}
}
}
}