模板矩阵加速

Posted 2021-02-05 winlere

矩阵加速，专门用来解决一些递推的关系，其原理和矩阵运算的法则有关

由于矩阵的乘法有交换律和结合律，所以我们可以通过矩阵快速幂来快速求解递推关系，一般时间复杂度是O(nlogn)。

矩阵快速幂很简单，写一下模板就会了，但是推导单位矩阵是个难题。

一般地，我们推导单位矩阵时，有这几个步骤。

1. 确定递推初始条件
2. 确定递推式子的系数
3. 确定单位矩阵的大小，一般来说和递推式涉及到的已知量个数有关，必须考虑系数为0的原量，不能忽略。
4. 填单位矩阵，注意一下几点：
   1. 第一行才是生成新数据的力量。
   2. 所以后面的行只需要把自己表示出来就好了。（为什么是用“行”来表示已知量？根据矩阵运算法则）

一般容易犯的错误：

• 递推关系写错了
• 单位矩阵错误了
• 矩阵乘法写错了
• 答案输出错了，没有找到到底应该输出哪一个
• 矩阵的指数确定错了
• %错了
• 没开long long导致乘法爆了

下贴斐波那契数列矩阵加速版本

#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define RP(t,a,b) for(ll t=(a),edd=(b);t<=edd;t++)
const int mod=1e9+7;
typedef unsigned long long ll;
const int n=2;
ll k;
struct mtx {
    ll data[3][3];
    ll* operator [](int x){return data[x];}
    mtx() {
        memset(data,0,sizeof data);
    }
    mtx operator *(const mtx& x) {
        mtx temp;
        RP(t,0,n-1)
        RP(i,0,n-1)
        RP(k,0,n-1)
        temp.data[t][i]=(temp.data[t][i]+(data[t][k]*x.data[k][i])%mod)%mod;
        return temp;
    }
    mtx operator *=(const mtx& x) {
        return (*this)=(*this)*x;
    }
    mtx operator ^(const ll& p) {
        ll b=p;
        mtx ans,base;
        ans.unis();
        base=(*this);
        while(b) {
            if(b&1)
                ans*=base;
            base*=base;
            b>>=1;
        }
        return ans;
    }
    mtx operator ^=(const ll& p) {
        return (*this)=(*this)^p;
    }
    void unis() {
        for(ll t=0;t<n;t++)
            data[t][t]=1;
    }
    void print() {
        for(ll t=0; t<n; t++) {
            for(ll ti=0; ti<n; ti++)
                cout<<data[t][ti]<<‘ ‘;
            cout<<"
";
        }
        return;
    }
};
int T;
int main() {
    //      freopen("in.in","r",stdin);
    cin>>k;
    mtx qaq;
    qaq.data[0][1]=qaq.data[1][0]=qaq.data[1][1]=1;
    qaq^=k-1;
    ll asa=qaq.data[0][0]%mod+qaq.data[1][0]%mod;
    cout<<asa%mod<<endl;
    return 0;
}

然后又附luogu矩阵加速板子的代码

#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define RP(t,a,b) for(ll t=(a),edd=(b);t<=edd;t++)
const int mod=1e9+7;
typedef long long ll;
const int n=3;
int k;
struct mtx {
    ll data[3][3];
    mtx() {
        memset(data,0,sizeof data);
    }
    mtx operator *(const mtx& x) {
        mtx temp;
        RP(t,0,n-1)
        RP(i,0,n-1)
        RP(k,0,n-1)
        temp.data[t][i]=(temp.data[t][i]+(data[t][k]*x.data[k][i])%mod)%mod;
        return temp;
    }
    mtx operator *=(const mtx& x) {
        return (*this)=(*this)*x;
    }
    mtx operator ^(const ll& p) {
        ll b=p;
        mtx ans,base;
        ans.unis();
        base=(*this);
        while(b) {
            if(b&1)
                ans*=base;
            base*=base;
            b>>=1;
        }
        return ans;
    }
    mtx operator ^=(const ll& p) {
        return (*this)=(*this)^p;
    }
    void unis() {
        for(ll t=0;t<n;t++)
            data[t][t]=1;
    }
    void print() {
        for(ll t=0; t<n; t++) {
            for(ll ti=0; ti<n; ti++)
                cout<<data[t][ti]<<‘ ‘;
            cout<<"
";
        }
        return;
    }
};
int T;
int main() {
    //  freopen("in.in","r",stdin);
    cin>>T;
    while(T--){
        cin>>k;
        mtx qaq;
        mtx temp;
        qaq.data[0][0]=qaq.data[0][2]=qaq.data[1][0]=qaq.data[2][1]=1;
        qaq^=k;
        int pos=0;
        cout<<qaq.data[1][0]<<endl;
    }
    return 0;
}
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define RP(t,a,b) for(ll t=(a),edd=(b);t<=edd;t++)
const int mod=1e9+7;
typedef long long ll;
const int n=3;
int k;
struct mtx {
    ll data[3][3];
    mtx() {
        memset(data,0,sizeof data);
    }
    mtx operator *(const mtx& x) {
        mtx temp;
        RP(t,0,n-1)
        RP(i,0,n-1)
        RP(k,0,n-1)
        temp.data[t][i]=(temp.data[t][i]+(data[t][k]*x.data[k][i])%mod)%mod;
        return temp;
    }
    mtx operator *=(const mtx& x) {
        return (*this)=(*this)*x;
    }
    mtx operator ^(const ll& p) {
        ll b=p;
        mtx ans,base;
        ans.unis();
        base=(*this);
        while(b) {
            if(b&1)
                ans*=base;
            base*=base;
            b>>=1;
        }
        return ans;
    }
    mtx operator ^=(const ll& p) {
        return (*this)=(*this)^p;
    }
    void unis() {
        for(ll t=0;t<n;t++)
            data[t][t]=1;
    }
    void print() {
        for(ll t=0; t<n; t++) {
            for(ll ti=0; ti<n; ti++)
                cout<<data[t][ti]<<‘ ‘;
            cout<<"
";
        }
        return;
    }
};
int T;
int main() {
    //  freopen("in.in","r",stdin);
    cin>>T;
    while(T--){
        cin>>k;
        mtx qaq;
        mtx temp;
        qaq.data[0][0]=qaq.data[0][2]=qaq.data[1][0]=qaq.data[2][1]=1;
        qaq^=k;
        int pos=0;
        cout<<qaq.data[1][0]<<endl;
    }
    return 0;
}

唉我太弱了orz