矩阵快速幂,从字面上来说就是矩阵之间的乘法运算关系,它分为四部分:
一、矩阵乘法函数
二、矩阵快速幂函数
三、最初矩阵的求解
四、结果的求解
第一和第二步有一些板子:
struct In{
ll m[maxn][maxn]; //定义结构体矩阵
};
In _mul(In x,In y,int n) //n × n的矩阵
{
In p; //存放 x矩阵 X y矩阵 的结果矩阵
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
p.m[i][j] = 0;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
p.m[i][j] = (p.m[i][j]+x.m[i][k]*y.m[k][j]) % mod; //关于是%mod
//还是%(mod-1)看你所求是否为质数,不是就%mod
return p;
}
In _pow(In x,ll N,int n){
In ans; //单位矩阵
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
if(i==j) ans.m[i][j]=1;
else ans.m[i][j]=0;
while(N){
if(N&1){
ans = _mul(ans,x);
}
x = _mul(x,x);
N >>= 1;
}
return ans;
}
第三步,最初矩阵的求解
在本文中就是 x.m矩阵的初始化问题:看了许多大佬的博客,他们都是直接写出这个初始化,我比较菜,就写一下求解初始化的矩阵的步骤:
以又见斐波那契数列为例:
就这个题目而言,首先求解 F(i) 会用到六个变量:F(i-1)、F(i-2)、i3、i2、i、1;所以他至少是一个6 x 6的矩阵,至于n等不等于6,就得看后期的求解;
我们会发现
F(i+1) = F(i) + F(i-1) + (i+1)3 + (i+1)2 + (i+1) + 1;
F(i) = 1 * F(i-1) + 1 * F(i-2) + 1 * i3 + 1 * i2 + 1 * i + 1 * 1;
F(i-1) = 1 * F(i-1) + 0 * F(i-2) + 0 * i3 + 0 * i2 + 0 * i + 0 * 1;
(i+1)3 = 0 * F(i-1) + 0 * F(i-2) + 1 * i3 + 3 * i2 + 3 * i + 1 * 1;
(i+1)2 = 0 * F(i-1) + 0 * F(i-2) + 0 * i3 + 1 * i2 + 2 * i + 1 * 1;
(i+1) = 0 * F(i-1) + 0 * F(i-2) + 0 * i3 + 0 * i2 + 2 * i + 1 * 1;
1 = 0 * F(i-1) + 0 * F(i-2) + 0 * i3 + 0 * i2 + 0 * i + 1 * 1;
所以很明显他是一个6×6的矩阵:
x.m的初始化:
1,1,1,1,1,1
1,0,0,0,0,0
0,0,1,3,3,1
0,0,0,1,2,1
0,0,0,0,1,1
0,0,0,0,0,1
四、结果的求解
还是以这个题为例子,F(i) 的值是有F(i-1)矩阵推导出来的,就拿F(2)来说,F(1)矩阵(记作F1[6][6]):
1,1,1,1,1,1
1,0,0,0,0,0
0,0,1,3,3,1
0,0,0,1,2,1
0,0,0,0,1,1
0,0,0,0,0,1
F(2) = 1 * F1[0][0] + 0 * F1[0][1] + 8 * F1[ 0][2]+ 4 * F1[0][3] + 2 * F1[0][4] + 1 * F1[0][5];
可以继续求,最终:
F(n) = 1 * Fn-1[0][0] + 0 * Fn-1[0][1] + 8 * Fn-1[ 0][2]+ 4 * Fn-1[0][3] + 2 * Fn-1[0][4] + 1 * Fn-1[0][5]
AC代码:
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstdio>
#define ll long long
#define mod 1000000007
using namespace std;
struct In{
ll m[6][6];
};
In _mul(In x,In y)
{
In p;
for(int i=0;i<6;i++)
for(int j=0;j<6;j++)
p.m[i][j] = 0;
for(int i=0;i<6;i++)
for(int j=0;j<6;j++)
for(int k=0;k<6;k++)
p.m[i][j] = (p.m[i][j]+x.m[i][k]*y.m[k][j])%mod;
return p;
}
In _pow(In x,ll n){
In ans;
for(int i=0;i<6;i++)
for(int j=0;j<6;j++)
if(i==j) ans.m[i][j]=1;
else ans.m[i][j]=0;
while(n){
if(n&1){
ans = _mul(ans,x);
}
x = _mul(x,x);
n >>= 1;
}
return ans;
}
int main(void)
{
ll nmb[6][6] =
{
{1,1,1,1,1,1},
{1,0,0,0,0,0},
{0,0,1,3,3,1},
{0,0,0,1,2,1},
{0,0,0,0,1,1},
{0,0,0,0,0,1}
};
ll n,t;
In res;
scanf("%lld",&t);
while(t--){
for(int i=0;i<6;i++)
for(int j=0;j<6;j++)
res.m[i][j] = nmb[i][j];
scanf("%lld",&n);
res = _pow(res,n-1);
if(n!=1&&n!=0){
ll op = (res.m[0][0]+res.m[0][2]*8+res.m[0][3]*4+res.m[0][4]*2+res.m[0][5])%mod;
printf("%lld\n",op);
}
else printf("%lld\n",n);
}
return 0;
}