【数论】矩阵加速

阅读本文之前 你应该有以下知识背景

 矩阵乘法法则 && 矩阵快速幂 && 小学数学水平 && 普及组代码能力 && 普及组思维能力

---

---

矩阵加速的本质是一个用了玄学优化的递推

根据这个解释 我们可以首先知道它是一个递推     //23333

普通的递推是通过不断f[i] = f[i-x] (某种运算)f[i-y]来完成的

举个简单例子  斐波那契数列 f[x]=f[x-1]+f[x-2]

那么由于f[x]由f[x-1]+f[x-2]决定 所以我们在求f[x] 之前必须知道f[x-1]和f[x-2]

当然了普通人都是这么想的 但是到了矩阵手里这玩意就不是这样了

想像一个过程 一个矩阵乘了自己好几遍之后 再乘以一个已给出的矩阵 就是你的答案

对于斐波那契来说就是 (一个矩阵)^k * (1,1) =answer

wow听起来很疯狂但是就是矩阵加速

---

---

怎么做矩阵加速？

首先来看斐波那契的矩阵加速//注意 本人习惯使用竖着的f[x][2]

你把递推边界放在一个矩阵里{1,1}然后乘以{{0,1}，{1,1}}//这个怎么来的待会会说

得到{1,2}

你再试几组 就会发现每次都会往后走一个 即{ f[x-2],f[x-1] } --->{ f[x-1],f[x] }

你可能会说这不是比递推还慢？

但是递推能快速幂吗？ 抱歉 矩阵能

所以你要递推n次 就求一个{{0,1},{1,1}}^n 然后把它乘以{1,1}

那么问题转化为 怎么构造转移矩阵

其实也就是在每个位置上体现递推规则

[x-2]->[x-1] 所以是{0,1}

[x-1]->[x]所以是{1,1}（[x] 和[x+1]相加）

OK讲解完毕 配几道题

斐波那契：

#include<bits/stdc++.h>
#define ll long long
#define p 1000000007

using namespace std;
ll k,n,s;
struct mat{
    ll m[3][2];
}f;//ans列 
struct smat{
    ll m[3][3];
}spd,e;

smat smul(smat x,smat y){
    smat r;
    for(int i=1;i<=2;i++){
        for(int j=1;j<=2;j++){
            r.m[i][j] = 0;
        }
    }
    for(int i=1;i<=2;i++){
        for(int j=1;j<=2;j++){
            for(int k=1;k<=2;k++){
                r.m[i][j] += x.m[i][k] * y.m[k][j] %p;
                r.m[i][j] %=p;
            }
        }
    }
    return r;
}

mat mul (smat x,mat y){
    mat r;
    for(int i=1;i<=2;i++){
        r.m[i][1]=0;
    }
    for(int i=1;i<=2;i++){
        for(int k=1;k<=2;k++){
            r.m[i][1] += x.m[i][k] * y.m[k][1] %p;
            r.m[i][1] %=p;
        }
    }
    return r;
}

smat pow(smat x,ll k){
    smat ans = e;
    while (k>0){
        if(k&1){
            ans=smul(ans,x);
        }
        x=smul(x,x);
        k>>=1;
    }
    return ans;
}

int main(){
    spd.m[1][2] = spd.m[2][1] = spd.m[2][2] = 1;
    e.m[1][1] = e.m[2][2] = 1;
    f.m[2][1] = f.m[1][1] = 1;
    cin>>n;
    smat vf=pow(spd,n-2);
    f=mul(vf,f);
    cout<<f.m[2][1];
    return 0;
}

板子 P1939 【模板】矩阵加速（数列）

其实大同小异吧。。。自己看代码

#include<bits/stdc++.h>
#define ll long long
#define p 1000000007

using namespace std;
ll k,n,s;
struct mat{
    ll m[4][2];
}f;//ans列 
struct smat{
    ll m[4][4];
}spd,e;

smat smul(smat x,smat y){
    smat r;
    for(int i=1;i<=3;i++){
        for(int j=1;j<=3;j++){
            r.m[i][j] = 0;
        }
    }
    for(int i=1;i<=3;i++){
        for(int j=1;j<=3;j++){
            for(int k=1;k<=3;k++){
                r.m[i][j] += x.m[i][k] * y.m[k][j] %p;
                r.m[i][j] %=p;
            }
        }
    }
    return r;
}

mat mul (smat x,mat y){
    mat r;
    for(int i=1;i<=3;i++){
        r.m[i][1]=0;
    }
    for(int i=1;i<=3;i++){
        for(int k=1;k<=3;k++){
            r.m[i][1] += x.m[i][k] * y.m[k][1] %p;
            r.m[i][1] %=p;
        }
    }
    return r;
}

smat pow(smat x,ll k){
    smat ans = e;
    while (k>0){
        if(k&1){
            ans=smul(ans,x);
        }
        x=smul(x,x);
        k>>=1;
    }
    return ans;
}

int main(){
    spd.m[1][2] = spd.m[2][3] = spd.m[3][1] = spd.m[3][3] = 1;
    e.m[1][1] = e.m[2][2] = e.m[3][3] = 1;
    f.m[3][1] = f.m[2][1] = f.m[1][1] = 1;
    ll wow;
    cin>>wow;
    for(int i=1;i<=wow;i++){
        cin>>n;
        
        smat vf=pow(spd,n-3);
        f=mul(vf,f);
        cout<<f.m[3][1]<<endl;
        f.m[3][1] = f.m[2][1] = f.m[1][1] = 1;
    }
    return 0;
}

本文是针对有快速幂基础的

所以没有那么多解释

TAG：SIN_XIII ⑨