矩阵快速幂

矩阵乘法是一种矩阵运算，满足交换律，可以写成幂的形式，所以我们可以使用矩阵快速幂来解决一些问题

先来看普通的快速幂：

快速幂

ll qpow(int a,int p){
    if(p == 0)return 1;
    ll temp = (qpow(a,p >> 1)) % M;	
    if(p % 2 == 0){
        return (temp * temp) % M;
        }
    else{
        return (temp * temp * a) % M;
        }
    }

把幂拆成 n * 2 或者n * 2 + 1的形式，避免重复运算，提高效率。因为结果通常很大，记得不要弄错MOD

矩阵乘法

设A为 的矩阵，B为 的矩阵，那么称 的矩阵C为矩阵A与B的乘积，记作 ，其中矩阵C中的第 行第 列元素可以表示为：

如下所示：

（摘自百度百科）

所以我们说：矩阵能相乘是有条件的：A矩阵的列数必须等于B的行数，得到的矩阵行数等于A的行数，列数等于B的列数

以下是矩阵乘法的代码

struct Mtrix{
   ll a[19][19];
   ll n,m;
   };
Mtrix X(Mtrix a,Mtrix b){
   Mtrix ans;
   for(ll i = 1;i <= a.n;i++)
       for(ll j = 1;j <= b.m;j++)
           ans.a[i][j] = 0;
   ans.n = a.n;
   ans.m = b.m;
   for(ll i = 1;i <= ans.n;i++){
       for(ll j = 1;j <= ans.m;j++){
           for(ll k = 1;k <= a.m;k++){
               ans.a[i][j] += (a.a[i][k] * b.a[k][j]) % M;
               ans.a[i][j] %= M;
               }
           }
       }
   return ans;
   }

矩阵快速幂

模仿一下一般快速幂，我们就可以得到矩阵快速幂的代码

Mtrix M_qpow(Mtrix a,ll p){
   if(p == 1)return a;
   Mtrix temp = M_qpow(a,p >> 1);
   if(p % 2 == 0)return X(temp,temp);
   else return X(a,X(temp,temp));
   }

我们不知道当p == 0时的矩阵值是多少，就干脆返回p == 1时为原矩阵就好

模板题：P3390 【模板】矩阵快速幂

利用矩阵快速幂优化递推问题

先看这题：

P1962 斐波那契数列

题目背景

大家都知道，斐波那契数列是满足如下性质的一个数列：

• f(1) = 1

• f(2) = 1

• f(n) = f(n-1) + f(n-2) (n ≥ 2 且 n 为整数)

题目描述

请你求出 f(n) mod 1000000007 的值。

输入输出格式

输入格式：
·第 1 行：一个整数 n

输出格式：
第 1 行： f(n) mod 1000000007 的值

说明

对于 60% 的数据： n ≤ 92

对于 100% 的数据： n在long long(INT64)范围内。

---

第一眼：水题！！

然后

题目：你以为呢？

瞟了一眼数据范围n在long long(INT64)范围内，吓个半死。

普通记忆化所需要的大量空间是肯定支持不了了，这时候就要用的矩阵快速幂了

首先依据斐波那契的定义式：

那么对于一个斐波那契矩阵：

第几项	值
n	f(n)
n + 1	f(n + 1)

我们可以这样写：

第几项	值
n	f(n) = f(n)
n + 1	f(n + 1) = f(n) + f(n - 1)

然后是不是可以由这个得到：

第几项	值
n - 1	f(n - 1)
n	f(n)

具体：

f(n)＝上面那个矩阵第一行×0＋第二行×1

f(n＋1)＝上面矩阵第一行×1＋第二行×1

怎么样，是不是特别熟悉？没错，矩阵运算！

它长这样：

0	1
1	1

这就是斐波那契的递推矩阵

再观察我们可以发现，f（n）其实就是f（1）；f（2）【纵向排列】乘以递推矩阵的（n - 1）次方，那么问题就解决了

代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<queue>
#define ll long long
using namespace std;
ll RD(){
    ll flag = 1,out = 0;char c = getchar();
    while(c < '0' || c > '9'){if(c == '-')flag = -1;c = getchar();}
    while(c >= '0' && c <= '9'){out = out * 10 + c - '0';c = getchar();}
    return flag * out;
    }
const ll M = 1000000007;
ll n;
struct Mtrix{
    ll a[19][19];
    ll n,m;
    };
Mtrix X(Mtrix a,Mtrix b){
    Mtrix ans;
    for(ll i = 1;i <= a.n;i++)
        for(ll j = 1;j <= b.m;j++)
            ans.a[i][j] = 0;
    ans.n = a.n;
    ans.m = b.m;
    for(ll i = 1;i <= ans.n;i++){
        for(ll j = 1;j <= ans.m;j++){
            for(ll k = 1;k <= a.m;k++){
                ans.a[i][j] += (a.a[i][k] * b.a[k][j]) % M;
                ans.a[i][j] %= M;
                }
            }
        }
    return ans;
    }
Mtrix M_qpow(Mtrix a,ll p){
    if(p == 1)return a;
    Mtrix temp = M_qpow(a,p >> 1);
    if(p % 2 == 0)return X(temp,temp);
    else return X(a,X(temp,temp));
    }
Mtrix a,ori;
int main(){
    n = RD();
    if(n == 1){
        cout<<1<<endl;
        return 0;
        }
    a.n = 2,a.m = 1;
    a.a[1][1] = 1,a.a[2][1] = 1;
    
    ori.n = 2,ori.m = 2;
    ori.a[1][1] = 0,ori.a[1][2] = 1,ori.a[2][1] = 1,ori.a[2][2] = 1;
    Mtrix last = M_qpow(ori,n - 1);
    Mtrix ans = X(last,a);
    cout<<ans.a[1][1] % M<<endl;
    return 0;
    }

总结

利用矩阵快速幂来求解递推是递推的一大法宝，要熟练运用