算法笔记--矩阵及矩阵快速幂

写的不错的博客：http://www.cnblogs.com/yan-boy/archive/2012/11/29/2795294.html

优点：根据数列递推式快速计算数列an的值（当n很大时），an可以是数，也可以是矩阵

步骤：由数列递推式构造矩阵，然后用矩阵快速幂计算矩阵的幂。

构造矩阵：

对于a_n =x*a_n-1 +y*a_n-2 ,可以构造矩阵为：

[a_n ,a_n-1]=[a_n-1 ,a_n-2]*A

A=[x 1

      y 0];

矩阵快速幂模板：

（2018.8.14修改）

struct Matrix {
    int a[3][3];
    void init() {
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++)
                a[i][j] = 0;
        }
    }
    void _init() {
        init();
        for (int i = 0; i < 3; i++) a[i][i] = 1;
    }
}A, B;

Matrix mul(Matrix a, Matrix b) {
    Matrix ans;
    ans.init();
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            if(a.a[i][j]) {
                for (int k = 0; k < 3; k++) ans.a[i][k] = (ans.a[i][k] + 1LL * a.a[i][j] * b.a[j][k]) % MOD;
            }
        }
    }
    return ans;
}
Matrix q_pow(Matrix a, int k) {
    Matrix ans;
    ans._init();
    if(k <= 0) return ans;
    while(k) {
        if(k&1) ans = mul(ans, a);
        a = mul(a, a);
        k >>= 1;
    }
    return ans;
}

例题1：POJ 3070 Fibonacci

代码：

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm> 
using namespace std;
#define ll long long
#define pb push_back
#define mem(a,b) memset(a,b,sizeof(a))

const int MOD=10000;
struct Matrix
{
    ll a[2][2];
}A,res;

Matrix multiply(Matrix x,Matrix y)
{
    Matrix res;
    mem(res.a,0);
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
            for(int k=0;k<2;k++)
                res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j])%MOD;
    return res;
}

void qpow(int n)
{
    while(n)
    {
        if(n&1)res=multiply(res,A);
        A=multiply(A,A);
        n=n>>1;
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int n;
    while(cin>>n)
    {
        if(n==-1)break;
        mem(res.a,0);
        for(int i=0;i<2;i++)res.a[i][i]=1;
        A.a[0][0]=1;
        A.a[0][1]=1;
        A.a[1][0]=1;
        A.a[1][1]=0;
        qpow(n);
        cout<<res.a[0][1]<<endl; 
    }
    return 0;
} 

例题2：HDU number number number

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pb push_back
#define mem(a,b) memset(a,b,sizeof(a))

const int MOD=998244353; 
struct Matrix 
{
    ll a[2][2];
}res,A;

Matrix mul(Matrix x,Matrix y)
{
    Matrix ans;
    mem(ans.a,0);
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
            for(int k=0;k<2;k++)
                ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%MOD;
    return ans;
}

void qpow(int n)
{
    mem(res.a,0);
    for(int i=0;i<2;i++)
    res.a[i][i]=1;
    A.a[0][0]=A.a[0][1]=A.a[1][0]=1;
    A.a[1][1]=0;
    while(n)
    {
        if(n&1)res=mul(res,A);
        A=mul(A,A);
        n>>=1;
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int k;
    while(cin>>k)
    {
        int n=k*2+3;
        qpow(n);
        cout<<res.a[0][1]-1<<endl;
    }
    return 0;
} 

 例题3：POJ Matrix Power Series

方法1：

构造矩阵：

B=(A E

O E)

B^k+1=(A^k+1 E+A+A² +A³ +...+A^k 

O                 E                    )

代码1：

#include<iostream>
#include<cstring>
using namespace std;
#define ll long long
#define pb push_back
#define mem(a,b) memset(a,b,sizeof(a))

int m,n;
struct Matrix 
{
    ll a[65][65];
}res,A;

Matrix mul(Matrix x,Matrix y)
{
    Matrix ans;
    mem(ans.a,0);
    for(int i=0;i<2*n;i++)
        for(int j=0;j<2*n;j++)
            for(int k=0;k<2*n;k++)
                ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%m;
    return ans;
}

void qpow(int k)
{
    mem(res.a,0);
    for(int i=0;i<2*n;i++)res.a[i][i]=1;
    for(int i=0;i<n;i++)A.a[i][i+n]=1;
    for(int i=n;i<2*n;i++)A.a[i][i]=1;
    while(k)
    {
        if(k&1)res=mul(res,A);
        A=mul(A,A);
        k>>=1;
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int k;
    cin>>n>>k>>m;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            cin>>A.a[i][j];
    qpow(k+1);
    
    for(int i=0;i<n;i++)
    {
        for(int j=n;j<2*n;j++)
        {
            if(j!=n+i)cout<<res.a[i][j];
            else cout<<(res.a[i][j]-1+m)%m;
            if(j!=2*n-1)cout<<' ';
        } 
        cout<<endl;
    }
    return 0;
} 

方法2：

由F_n =A*F_n-1 +A构造矩阵：

(F_n ,E)=(F_n-1,E)*B

B=(A O

A E)

(F_n ,E)=(F₁,E)*B^n-1 

代码2：

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm> 
using namespace std;
#define ll long long
#define pb push_back
#define mem(a,b) memset(a,b,sizeof(a))

int n,m,k;
struct Matrix
{
    ll a[65][65];
}A,res,t;

Matrix mul(Matrix x,Matrix y)
{
    Matrix ans;
    mem(ans.a,0);
    for(int i=0;i<2*n;i++)
        for(int j=0;j<2*n;j++)
            for(int k=0;k<2*n;k++)
                ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%m;
    return ans;
}

void qpow(int k)
{
    mem(res.a,0);
    mem(A.a,0);
    for(int i=0;i<2*n;i++)res.a[i][i]=1; 
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            A.a[i][j]=A.a[i+n][j]=t.a[i][j];
    for(int j=n;j<2*n;j++)A.a[j][j]=1;
       while(k)
    {
        if(k&1)res=mul(res,A);
        A=mul(A,A);
        k=k>>1;
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    cin>>n>>k>>m;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            cin>>t.a[i][j];
    qpow(k-1);
    for(int i=0;i<n;i++)
        t.a[i][i+n]=1;
    /*for(int i=0;i<2*n;i++)
        for(int j=0;j<2*n;j++)
        {
            cout<<res.a[i][j]<<' ';
            if(j==2*n-1)cout<<endl;
        }*/
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            int ans=0;
            for(int k=0;k<2*n;k++)
            {
                ans=(ans+t.a[i][k]*res.a[k][j])%m;
            }
            cout<<ans;
            if(j!=n-1)cout<<' ';
        }
        cout<<endl;
    }
    return 0;
}

例题4：POJ Training little cats

写的不错的博客：http://blog.csdn.net/magicnumber/article/details/6217602

要用到一个对于稀疏矩阵乘法的加速技巧

代码:

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm> 
using namespace std;
#define ll long long
#define pb push_back
#define mem(a,b) memset(a,b,sizeof(a))

int n,m,k;
struct Matrix
{
    ll a[105][105];
}A,res,unit;

Matrix mul(Matrix x,Matrix y)
{
    Matrix ans;
    mem(ans.a,0);
    for(int i=0;i<n+1;i++)
        for(int j=0;j<n+1;j++)
        {
            if(x.a[i][j])
            for(int k=0;k<n+1;k++)
            ans.a[i][k]+=x.a[i][j]*y.a[j][k];
            /*for(int k=0;k<n+1;k++)
            ans.a[i][j]+=x.a[i][k]*y.a[k][j];*/
        } 
    return ans;
}

void init()
{
    mem(unit.a,0);
    for(int i=0;i<n+1;i++)
    {
        unit.a[i][i]=1;
    }
}

void qpow(int k)
{
    init();
    res=unit;
    while(k)
    {
        if(k&1)res=mul(res,A);
        A=mul(A,A);
        k>>=1;
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    while(cin>>n>>m>>k)
    {
        if(n==0&&m==0&&k==0)break; 
        int i,j;
        mem(A.a,0);
        char s;
        init();
        res=unit;
        
        while(k--)
        {
            cin>>s;
            if(s=='g')    
            {
                cin>>i;
                init();
                unit.a[i-1][n]=1;
                res=mul(unit,res);
            }
            else if(s=='e')
            {
                cin>>i;
                init();
                unit.a[i-1][i-1]=0;
                res=mul(unit,res);
            }
            else if(s=='s')
            {
                cin>>i>>j;
                init();
                unit.a[i-1][i-1]=0;
                unit.a[i-1][j-1]=1;
                unit.a[j-1][j-1]=0;
                unit.a[j-1][i-1]=1;
                res=mul(unit,res);
            }
        }
        A=res;
        
        qpow(m);
        for(int i=0;i<n;i++)
        {
            cout<<res.a[i][n];
            if(i!=n-1)cout<<' ';
            else cout<<endl;
        }
    }  
    return 0;
}

例题5：HDU - 5015

首先，如果没有最上面一行，每新增加一列都是求一遍前缀和，那么矩阵就是一个下三角全为1的矩阵

但是，这道题最上面有数，那我们矩阵稍微改一改

10 0 0 0 ... 1           23

10 1 0 0 ... 1           a1

10 1 1 0 ... 1           a2

10 1 1 1 ... 1     *    a3

......                          .

10 1 1 1 ... 1           an

0   0 0 0 ... 1            3

这个相当于算第一列，算m列就是乘m次左边的那个矩阵

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize(4)
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define pi acos(-1.0)
#define LL long long
//#define mp make_pair
#define pb push_back
#define ls rt<<1, l, m
#define rs rt<<1|1, m+1, r
#define ULL unsigned LL
#define pll pair<LL, LL>
#define pii pair<int, int>
#define piii pair<pii, int>
#define mem(a, b) memset(a, b, sizeof(a))
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define fopen freopen("in.txt", "r", stdin);freopen("out.txt", "w", stout);
//head

const int MOD = 1e7 + 7;
int a[100], n;
struct Matrix {
    int a[15][15];
    void init() {
        for (int i = 0; i <= n+1; i++) {
            for (int j = 0; j <= n+1; j++) {
                a[i][j] = 0;
            }
        }
    }
    void _init() {
        init();
        for (int i = 0; i <= n+1; i++) a[i][i] = 1;
    }
}T;
Matrix mul(Matrix A, Matrix B) {
    Matrix ans;
    ans.init();
    for (int i = 0; i <= n+1; i++) {
        for (int j = 0; j <= n+1; j++) {
            if(A.a[i][j])
                for (int k = 0; k <= n+1; k++)
                ans.a[i][k] = (ans.a[i][k] + 1LL * A.a[i][j] * B.a[j][k]) % MOD;
        }
    }
    return ans;
}
Matrix q_pow(Matrix A, int k) {
    Matrix ans;
    ans._init();
    while(k) {
        if(k&1) ans = mul(ans, A);
        A = mul(A, A);
        k >>= 1;
    }
    return ans;
}
void init() {
    T.a[0][0] = 10;
    for (int i = 1; i <= n; i++) T.a[i][0] = 10;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= i; j++)
            T.a[i][j] = 1;
    }
    for (int i = 0; i <= n+1; i++) T.a[i][n+1] = 1;
}
int main() {
    int m;
    while(~ scanf("%d %d", &n, &m)) {
        for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
        T.init();
        init();
        T = q_pow(T, m);
        int ans = T.a[n][0] * 23LL % MOD;
        for (int i = 1; i <= n; i++) ans = (ans + 1LL*T.a[n][i]*a[i]) % MOD;
        ans = (ans + 1LL*T.a[n][n+1]*3LL) % MOD;
        printf("%d
", ans);
    }
    return 0;
}