浅谈矩阵快速幂

矩阵快速幂，从字面上来说就是矩阵之间的乘法运算关系，它分为四部分：
  一、矩阵乘法函数
  二、矩阵快速幂函数
  三、最初矩阵的求解
  四、结果的求解
第一和第二步有一些板子：

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)、i³、i²、i、1；所以他至少是一个6 x 6的矩阵，至于n等不等于6，就得看后期的求解；
我们会发现
  F(i+1) = F(i) + F(i-1) + （i+1)³ + (i+1)² + (i+1) + 1;
  F(i) = 1 * F(i-1) + 1 * F(i-2) + 1 * i³ + 1 * i² + 1 * i + 1 * 1;
  F(i-1) = 1 * F(i-1) + 0 * F(i-2) + 0 * i³ + 0 * i² + 0 * i + 0 * 1;
  (i+1)³ = 0 * F(i-1) + 0 * F(i-2) + 1 * i³ + 3 * i² + 3 * i + 1 * 1;
  (i+1)² = 0 * F(i-1) + 0 * F(i-2) + 0 * i³ + 1 * i² + 2 * i + 1 * 1;
  (i+1) = 0 * F(i-1) + 0 * F(i-2) + 0 * i³ + 0 * i² + 2 * i + 1 * 1;
  1 = 0 * F(i-1) + 0 * F(i-2) + 0 * i³ + 0 * i² + 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;
}