矩阵快速幂

矩阵快速幂通常可以用来加速计算递推关系（因为可以使用快速幂来求矩阵的乘法）。

矩阵快速幂主要由两个函数组成：

#include<iostream>
#include<cstdio>
#include<set>
#include<map>
#include<vector>
#include<queue>
#include<stack>
#include<cmath>
#include<cstring>
#include<algorithm>

typedef vector<ll> vec;
typedef vector<vec> M;

M mul(M a,M b)  //矩阵a乘矩阵b
{
    int r=a.size();
    M c(r,vec(r,0));

    for(int i=0;i<r;i++)
    {
        for(int j=0;j<r;j++)
        {
            for(int k=0;k<r;k++)
                c[i][j]=(c[i][j]+a[i][k]*b[k][j])%mod;
        }
    }
    return c;
}

M qm(M a,int n)  //使用快速幂求矩阵a的n次方
{
    int r=a.size();
    M b(r,vec(r,0));
    for(int i=0;i<r;i++)b[i][i]=1;
    while(n)
    {
        if(n&1)b=mul(b,a);
        a=mul(a,a);
        n>>=1;
    }
    return b;
}

例题：

题目1：Queuing

代码：

#include<iostream>
#include<cstdio>
#include<vector>
#define ll long long
using namespace std;

int mod;

typedef vector<ll> vec;
typedef vector<vec> N;

N mul(N a,N b)
{
    int r=a.size(),d=a[0].size();
    N c(r,vec(d,0));
    for(int i=0;i<r;i++)
    {
        for(int j=0;j<r;j++)
        {
            ll s=0;
            for(int k=0;k<d;k++)
                c[i][j]=(c[i][j]+(a[i][k]*b[k][j])%mod)%mod;
        }
    }
    return c;
}

N qmul(N a,int n)
{
    int r=a.size(),c=a[0].size();
    N b(r,vec(c));
    for(int i=0;i<r;i++)
        for(int j=0;j<c;j++)
        b[i][j]=(i==j);
    while(n)
    {
        if(n&1)b=mul(b,a);
        n>>=1;
        a=mul(a,a);
    }
    return b;
}

int main()
{
    int n;
    while(scanf("%d%d",&n,&mod)==2)
    {
        if(n==0)printf("0\n");
        else if(n==1)printf("2\n");
        else
        {
            N a(4,vec(4,0));
            a[0][0]=a[0][2]=a[1][0]=a[2][1]=a[2][3]=a[3][1]=1;
            a=qmul(a,n-2);
            ll ans=0;
            for(int i=0;i<4;i++)
                for(int j=0;j<4;j++)
                ans+=a[i][j];
            ans%=mod;
            printf("%d\n",ans);
        }
    }
    return 0;
}

题目2：Problem of Precision

代码：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<vector>
#define ll long long
using namespace std;

const int mod=1024;

typedef vector<int> vec;
typedef vector<vec> M;

M mul(M a,M b)
{
    int r=a.size();
    M c(r,vec(r,0));

    for(int i=0;i<r;i++)
    {
        for(int j=0;j<r;j++)
        {
            for(int k=0;k<r;k++)
                c[i][j]=(c[i][j]+a[i][k]*b[k][j]+mod)%mod;
        }
    }
    return c;
}

M qm(M a,int n)
{
    int r=a.size();
    M b(r,vec(r,0));
    for(int i=0;i<r;i++)b[i][i]=1;
    while(n)
    {
        if(n&1)b=mul(b,a);
        a=mul(a,a);
        n>>=1;
    }
    return b;
}

int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        int n;
        scanf("%d",&n);

        M a(2,vec(2));
        a[0][0]=5;a[0][1]=12;
        a[1][0]=2;a[1][1]=5;

        a=qm(a,n-1);
        int an=(a[0][0]*5+a[0][1]*2);

        int ans=(2*an-1)%mod;

        printf("%d\n",ans);
    }
    return 0;
}

题目3：Chinese Rings

#include<iostream>
#include<cstdio>
#include<cmath>
#include<vector>
#define ll long long
using namespace std;

const int mod=200907;

typedef vector<ll> vec;
typedef vector<vec> M;

M mul(M a,M b)
{
    int r=a.size();
    M c(r,vec(r,0));

    for(int i=0;i<r;i++)
    {
        for(int j=0;j<r;j++)
        {
            for(int k=0;k<r;k++)
                c[i][j]=(c[i][j]+a[i][k]*b[k][j])%mod;
        }
    }
    return c;
}

M qm(M a,int n)
{
    int r=a.size();
    M b(r,vec(r,0));
    for(int i=0;i<r;i++)b[i][i]=1;
    while(n)
    {
        if(n&1)b=mul(b,a);
        a=mul(a,a);
        n>>=1;
    }
    return b;
}

int main()
{
    int n;
    while(scanf("%d",&n)&&n)
    {
        M a(3,vec(3,0));
        a[0][1]=2;
        a[0][0]=a[0][2]=a[1][0]=a[2][2]=1;
        a=qm(a,n-1);

        ll ans=(a[0][0]+a[0][2]+mod)%mod;

        printf("%lld\n",ans);
    }
    return 0;
}

题目4：Kiki & Little Kiki 2

#include<iostream>
#include<cstdio>
#include<cmath>
#include<vector>
#include<cstring>
#define ll long long
using namespace std;

const int mod=2;

typedef vector<int> vec;
typedef vector<vec> M;

char s[120];

M mul(M a,M b)
{
    int r=a.size();
    M c(r,vec(r,0));

    for(int i=0; i<r; i++)
    {
        for(int j=0; j<r; j++)
        {
            for(int k=0; k<r; k++)
                c[i][j]=(c[i][j]+a[i][k]*b[k][j])%mod;
        }
    }
    return c;
}

M qm(M a,int n)
{
    int r=a.size();
    M b(r,vec(r,0));
    for(int i=0; i<r; i++)b[i][i]=1;
    while(n)
    {
        if(n&1)b=mul(b,a);
        a=mul(a,a);
        n>>=1;
    }
    return b;
}

int main()
{
    int m;
    while(scanf("%d%s",&m,s)==2)
    {
        int l=strlen(s);
        M a(l,vec(l,0));
        a[0][0]=a[0][l-1]=1;
        for(int i=1;i<l;i++)a[i][i]=a[i][i-1]=1;
        a=qm(a,m);
        for(int i=0;i<l;i++)
        {
            int ans=0;
            for(int j=0;j<l;j++)
                ans+=(s[j]-'0')*a[i][j];
            printf("%d",ans%2);
        }
        printf("\n");
    }
    return 0;
}

解这类题的关键是用矩阵表示递推关系。

 注意：注意取摸以及矩阵元素的数据类型（分析是否会爆int）