矩阵快速幂

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

先来看普通的快速幂:

快速幂

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;
    }

总结

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

posted @ 2018-07-09 12:33  Tony_Double_Sky  阅读(292)  评论(0编辑  收藏  举报