训练第四周之矩阵快速幂

矩阵快速幂 = 矩阵乘法 + 快速幂

我觉得矩阵快速幂的题目不是很难，只要通过递推找到关系矩阵A（构造任意线性递推式对应的矩阵），然后求A的k次幂乘以题目所给数据就能得到所求。

1、矩阵乘法：当第一个矩阵A的列数和另一个矩阵B的行数相等时才可以相乘，复杂度位O(n^3)

代码模板：

Mat MulMat(Mat a,Mat b)///a与b进行矩阵相乘运算,Mat为自己定义的结构体变量类型
{
    Mat c={0};
    int i,j,l;
    for(i=0;i<n;i++)//n为矩阵行数
        for(j=0;j<n;j++)
            for(l=0;l<n;l++)
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][l]*b.mat[l][j])%mod;//题目一般要mod，不加也行
     //这里之前写成c.mat[i][j]+=(a.mat[i][l]*b.mat[l][j])%mod;就wa了，估计爆数据了，以后这样写保险一点
    return c;
}

附加矩阵加法模板：

Mat AddMat(Mat a,Mat b)//a与b进行矩阵相加运算
{
    Mat c={0};
    int i,j;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
            c.mat[i][j]=(c.mat[i][j]+a.mat[i][j]+b.mat[i][j])%mod;
    }
    return c;
}

2、快速幂：如求k^n的后5位，当n的范围很小的时候可以暴力，复杂度是O(n)，但当n=1e18时就爆时间了，所以就用到快速幂

假设你已经有一个函数F(x) = k^x它可以求 k^1 直到 k^(n-1)，如何求F(n)？

当n是偶数的时候：F(n)=F(n/2)*F(n/2)

当n是奇数的时候：转化成偶数的情况 F(n) = k * F(n - 1)

这样复杂度就变成O(logn)了，好像没很大改进，但是当n为1亿的时候，O(n)为100000000，而O(logn)是26。

伪代码：

int quickmod(int a, int b){  
    ans = 1; 
    只要b不等于0：
        如果b是奇数：
              ans = (ans*a)%mod;  b--；
        如果b是偶数：
        b/=2;  
        a = ((a%mod)*(a%mod))%mod;               
return ans;  
} 

3、矩阵快速幂:就是把快速幂的数乘改为矩阵乘法就行

代码模板：

Mat fast_mod(Mat a,long long t)//t尽量用long long，不然很容易MLE，但也要看题意，有时int够了
{
    if(t==0)
        return E;//E为单位矩阵
    if(t%2)
        return MulMat(a,fast_mod(a,t-1));
    Mat p={0};
    p=fast_mod(a,t>>1);
    return MulMat(p,p);
}

复杂度：暴力的矩阵求幂是O(t*n^3)，而矩阵快速幂是O(n^3 * log t)

例题：

1、Matrix Power Series-poj3233

思路：这个题目我觉得比较难，因为题目没有给你一个固定的矩阵，就很难用递推式推出关系矩阵A。所以只能找其他关系：求Sn,t为偶数,则Sn=(E+A^n/2)* Sn/2,奇数则Sn= (E+A^n/2)*Sn/2+A^n

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
struct Mat
{
    int mat[35][35];
};
Mat E={0},A={0};
int m,n;
Mat MulMat(Mat a,Mat b)//a与b进行矩阵相乘运算
{
    Mat c={0};
    int i,j,l;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            for(l=0;l<n;l++)
            {   //这里之前写成c.mat[i][j]+=(a.mat[i][l]*b.mat[l][j])%m;就wa了，估计爆数据了，以后这样写保险一点
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][l]*b.mat[l][j])%m;
            }
            //c.mat[i][j] %= m;
            //printf("%d ",c.mat[i][j]);
        }
        //printf("\n");
    }
    return c;
}
Mat AddMat(Mat a,Mat b)//a与b进行矩阵相加运算
{
    Mat c={0};
    int i,j;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
            c.mat[i][j]=(c.mat[i][j]+a.mat[i][j]+b.mat[i][j])%m;
    }
    return c;
}
Mat fast_mod(Mat a,int t)//利用二分法求矩阵快速幂求A^t
{
    if(t==0)
        return E;
    if(t%2)
        return MulMat(a,fast_mod(a,t-1));
    Mat p={0};
    p=fast_mod(a,t>>1);
    return MulMat(p,p);
}
Mat operater(Mat a,int t)//求Sn,t为偶数,则Sn=(E+A^n/2)*Sn/2,奇数则Sn=(E+A^n/2)*Sn/2+A^n
{
    if(t==1)
        return a;
    if(t%2)
    {
        return AddMat(MulMat(AddMat(E,fast_mod(a,t/2)),operater(a,t/2)),fast_mod(a,t));
    }
    Mat p={0};
    p=operater(a,t>>1);
    return MulMat(AddMat(E,fast_mod(a,t/2)),p);
}
int main()
{
    int k,i,j;
    cin>>n>>k>>m;
    for(i=0;i<35;i++)
        E.mat[i][i]=1;//单位矩阵
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            cin>>A.mat[i][j];
    Mat Sk={0};
    Sk=operater(A,k);
    for(i=0;i<n;i++)
    {
        for(j=0;j<n-1;j++)//输出小技巧
            cout<<Sk.mat[i][j]%m<<" ";
        cout<<Sk.mat[i][n-1]%m<<endl;
    }
    return 0;
}

2、Fibonacci-poj3070

思路：这个题目比较中规中矩一些，题目给出了f[0]=0,f[1]=1,通过递推式构造关系矩阵A.mat[2][2]={1,1,1,0}，然后用矩阵快速幂求出A^n

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
struct Mat
{
    long long mat[2][2];
};
Mat A={1,1,1,0};
Mat e={1,0,0,1};
long long n;
Mat MulMat(Mat a,Mat b)
{
    Mat c={0,0,0,0};
    int i,j,t;
    for(i=0;i<2;i++)
    {
        for(j=0;j<2;j++)
        {
            for(t=0;t<2;t++)
            {
                c.mat[i][j]+=a.mat[i][t]*b.mat[t][j];
            }
            c.mat[i][j] %= 10000 ;
            //printf("%d ",c.mat[i][j]);
        }
        //printf("\n");
    }
    return c;
}
Mat fast_mod(Mat a,long long n)
{
    if(n==0)
        return e;
    if(n%2)
    {
        return MulMat(a,fast_mod(a,n-1));
    }
    else
    {
        Mat p;
        p=fast_mod(a,n>>1);
        return MulMat(p,p);
    }
}
int main()
{
    Mat f={0,0,1,0},v;
    int a[4],i;
    long long fn=0;
    while(cin>>n&&n!=-1)
    {
        v=fast_mod(A,n);
        fn=v.mat[0][1];
        if(fn<999)
        {
            cout<<fn<<endl;
            continue;
        }
        a[0]=fn%10;
        a[1]=fn/10%10;
        a[2]=fn/10/10%10;
        a[3]=fn/10/10/10%10;
        for(i=3;i>=0;i--)
        {
            if(a[i]!=0)
                break;
        }
        if(i==-1)
            cout<<"0"<<endl;
        else
        {
            for(int j=i;j>=0;j--)
                cout<<a[j];
            cout<<endl;
        }
    }
    return 0;
}

3、Number Sequence-hdu1005

思路：与第二题类似，直接给代码了

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define m 7
struct Mat
{
    int mat[2][2];
};
Mat A={0};
Mat e={1,0,0,1};
int n;
Mat MulMat(Mat a,Mat b)
{
    Mat c={0,0,0,0};
    int i,j,t;
    for(i=0;i<2;i++)
    {
        for(j=0;j<2;j++)
        {
            for(t=0;t<2;t++)
            {
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][t]*b.mat[t][j])%m;
            }
            //printf("%d ",c.mat[i][j]);
        }
        //printf("\n");
    }
    return c;
}
Mat fast_mod(Mat a,int n)
{
    if(n==0)
        return e;
    if(n%2)
    {
        return MulMat(a,fast_mod(a,n-1));
    }
    else
    {
        Mat p;
        p=fast_mod(a,n>>1);
        return MulMat(p,p);
    }
}
int main()
{
    int i,a,b;
    while(cin>>a>>b>>n&&(a!=0||b!=0||n!=0))//不看清题目，乱用&&导致4次MLE了
    {
        Mat f={0};
        A.mat[0][0]=a;
        A.mat[0][1]=b;
        A.mat[1][0]=1;
        if(n==1||n==2)
        {
            cout<<"1"<<endl;
            continue;
        }
        f=fast_mod(A,n-2);
        cout<<(f.mat[0][0]+f.mat[0][1])%m<<endl;
    }
    return 0;
}

4、 Plant-CodeForces 185A

代码：

#include<cstdio>
#include<cstring>
#include<iomanip>
#include<iostream>
#include<algorithm>
using namespace std;
struct Mat
{
    unsigned long long mat[2][2];
};
Mat A={2,1,0,4};
Mat E={1,0,0,1};
Mat p={0,0,0,0};
Mat MulMat(Mat a,Mat b)
{
    Mat c={0,0,0,0};
    int i,j,s;
    for(i=0;i<2;i++)
    {
        for(j=0;j<2;j++)
        {
            for(s=0;s<2;s++)
            {
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][s]*b.mat[s][j])%1000000007;
            }
        }
    }
    return c;
}
Mat fast_mod(long long n)//这里之前就用了int，结果MLE了两次
{
    if(n==0)
        return E;
    if(n%2)
        return MulMat(A,fast_mod(n-1));
    p=fast_mod(n>>1);
    return MulMat(p,p);
}
int main()
{
    unsigned long long n;
    unsigned long long s;
    while(cin>>n)
    {
        Mat f={0,0,0,0};
        f=fast_mod(n);
        s=(f.mat[0][1]+f.mat[0][0])%1000000007;
        cout<<s<<endl;
    }
    return 0;
}