• 矩阵的快速幂


    (1)矩阵乘法

    简单的说矩阵就是二维数组,数存在里面,矩阵乘法的规则:A*B=C

    其中c[i][j]为A的第i行与B的第j列对应乘积的和,即:

     

      据说,矩阵快速幂在递推式优化上相当神奇,而且效率很高。。。

      两矩阵相乘,朴素算法的复杂度是O(N^3)。如果求一次矩阵的M次幂,按朴素的写法就是O(N^3*M)。既然是求幂,不免想到快速幂取模的算法,这里有快速幂取模的介绍,a^b %m 的复杂度可以降到O(logb)。如果矩阵相乘是不是也可以实现O(N^3 * logM)的时间复杂度呢?答案是肯定的。

     模板一

    利用快速幂的思想 根据矩阵的结合律 可以递归二分求解 
    
    struct Mat
    {
        int mat[N][N];
    };
    int n;
    Mat operator * (Mat a,Mat b)
    {
        Mat c;
        memset(c.mat,0,sizeof(c.mat));
        int i,j,k;
        for(k =0 ; k < n ; k++)
        {
            for(i = 0 ; i < n ;i++)
            {
                if(a.mat[i][k]==0) continue;//优化
                for(j = 0 ;j < n ;j++)
                {
                    if(b.mat[k][j]==0) continue;//优化
                    c.mat[i][j] = (c.mat[i][j]+(a.mat[i][k]*b.mat[k][j])%mod)%mod;
                }
            }
        }
        return c;
    }
    Mat operator ^(Mat a,int k)
    {
        Mat c;
        int i,j;
        for(i =0 ; i < n ;i++)
            for(j = 0; j < n ;j++)
            c.mat[i][j] = (i==j);
        for(; k ;k >>= 1)
        {
            if(k&1) c = c*a;
            a = a*a;
        }
        return c;
    }
    View Code

     模板二

    #include<stdio.h>
    #include<vector>
    using namespace std;
    #define inf 0x3f3f3f3f
    typedef vector<int>vec;
    typedef vector<vec>mat;
    typedef long long ll;
    const int M = 10000;
    ll n;
    mat mul(mat &A , mat &B)
    {
        mat C(A.size() , vec(B.size()));
        for(int i=0 ; i<A.size() ; i++)
        {
            for(int k=0 ; k<B.size() ; k++)
            {    
                if(A[i][k]==0)
                continue;
                for(int j=0 ; j<B[0].size() ; j++)
                {   
                    if(B[k][j]==0)
                    continue;
                    C[i][j] = (C[i][j]+A[i][k]*B[k][j])%M;
                }
            }
        }
        return C;
    }
    
    mat pow(mat A,ll n)
    {
        mat B(A.size(),vec(A.size()));
        for(int i=0 ; i<A.size() ; i++)
            B[i][i]=1;
        while(n>0)
        {
            if(n&1)
            B = mul(B,A);
            A = mul(A,A);
            n  >>= 1;
        }
        return B;
    }
    View Code

     应用性模板:

      | 1 5 1 -1 |       | F[n-1] |           | F[n]  |

      | 1 0 0 0 |    *   | F[n-2] |    =     | F[n-1] |

      | 0 1 0 0 |        |F[n-3] |             | F[n-2] |

      | 0 0 1 0 |        | F[n-4]|            | F[n-3] |

    f1=1 , f2=5 , f3=11 , f4=36 , f5=95;

    f(n)=f(n-1)+5*f(n-2)+f(n-3)-f(n-4)

    #include <cstdio>
    #include <vector>
    #include <queue>
    #include <cstring>
    #include <algorithm>
    using namespace std;
    #define mst(a,b) memset((a),(b),sizeof(a))
    #define rush() int T,scanf("%d",&T),while(T--)
     
    typedef long long ll;
    const int maxn = 4;
    const ll mod = 1000000007;
    const int INF = 0x3f3f3f3f;
     
    struct Matrix
    {
        ll temp[maxn][maxn];
    } a;
     
    void init()
    {
        for(int i=0;i<maxn;i++)
        for(int j=0;j<maxn;j++)
        {
            a.temp[i][j]=0;
        }
        a.temp[0][0]=a.temp[0][2]=1;
        a.temp[0][1]=5;
        a.temp[0][3]=-1;
        a.temp[1][0]=a.temp[2][1]=a.temp[3][2]=1;
    }
    Matrix mul(Matrix a,Matrix b)
    {
        Matrix ans;
        for (int i=0; i<maxn; i++)
            for (int j=0; j<maxn; j++)
            {
                ans.temp[i][j]=0;
                for (int k=0; k<maxn; k++)
                {
                    ans.temp[i][j]+=(a.temp[i][k]*b.temp[k][j]+mod)%mod;  //特别注意
                    ans.temp[i][j]%=mod;
                }
            }
        return ans;
    }
     
    void fun(Matrix ans,ll k)
    {
        for(int i=0; i<maxn; i++)
            for(int j=0; j<maxn; j++)
                a.temp[i][j]=(i==j);
        while(k)
        {
            if(k%2)
                a=mul(a,ans);
            ans=mul(ans,ans);
            k/=2;
        }
    }
     
    int main()
    {
        Matrix t;
        ll n;
        for(int i=0;i<maxn;i++)
        for(int j=0;j<maxn;j++)
        {
            t.temp[i][j]=0;
        }
        t.temp[0][0]=36;
        t.temp[1][0]=11;
        t.temp[2][0]=5;
        t.temp[3][0]=1;
        while(~scanf("%I64d",&n))
        {
            init();
            if(n<=4)
            {
                printf("%I64d
    ",t.temp[4-n][0]);
                continue;
            }
            fun(a,n-4);
            a=mul(a,t);
            ll ans=a.temp[0][0]%mod;
            printf("%I64d
    ",ans);
        }
        return 0;
    }
    View Code

    应用篇

    主要通过把数放到矩阵的不同位置,然后把普通递推式变成"矩阵的等比数列",最后快速幂求解递推式:

    先通过入门的题目来讲应用矩阵快速幂的套路(会这题的也可以看一下套路):

    例一:http://poj.org/problem?id=3070
    题目:斐波那契数列f(n),给一个n,求f(n)%10000,n<=1e9;

    (这题是可以用fibo的循环节去做的,不过这里不讲,反正是水题)

    矩阵快速幂是用来求解递推式的,所以第一步先要列出递推式:

     f(n)=f(n-1)+f(n-2)

    第二步是建立矩阵递推式,找到转移矩阵:

    ,简写成T * A(n-1)=A(n),T矩阵就是那个2*2的常数矩阵,而

    这里就是个矩阵乘法等式左边:1*f(n-1)+1*f(n-2)=f(n);1*f(n-1)+0*f(n-2)=f(n-1);

    这里还是说一下构建矩阵递推的大致套路,一般An与A(n-1)都是按照原始递推式来构建的,当然可以先猜一个An,主要是利用矩阵乘法凑出矩阵T,第一行一般就是递推式,后面的行就是不需要的项就让与其的相乘系数为0。矩阵T就叫做转移矩阵(一定要是常数矩阵),它能把A(n-1)转移到A(n);然后这就是个等比数列,直接写出通项:此处A1叫初始矩阵。所以用一下矩阵快速幂然后乘上初始矩阵就能得到An,这里An就两个元素(两个位置),根据自己设置的A(n)对应位置就是对应的值,按照上面矩阵快速幂写法,res[1][1]=f(n)就是我们要求的。

     AC代码:

    #include<stdio.h>
    #include<string.h>
    #define mod 10000
    struct Mat
    {
        long long  mat[2][2];
    };
    int n=2;
    Mat operator * (Mat a,Mat b)
    {
        Mat c;
        c.mat[0][0]=c.mat[0][1]=c.mat[1][0]=c.mat[1][1]=0;
        int i,j,k;
        for(k =0 ; k < n ; k++)
        {
            for(i = 0 ; i < n ;i++)
            {
                if(a.mat[i][k]==0) continue;//优化
                for(j = 0 ;j < n ;j++)
                {
                    if(b.mat[k][j]==0) continue;//优化
                    c.mat[i][j] = (c.mat[i][j]+(a.mat[i][k]*b.mat[k][j])%mod)%mod;
                }
            }
        }
        return c;
    }
    Mat operator ^(Mat a,int k)
    {
        Mat c;
        int i,j;
        for(i =0 ; i < n ;i++)
            for(j = 0; j < n ;j++)
            c.mat[i][j] = (i==j);
        for(; k ;k >>= 1)
        {
            if(k&1) c = c*a;
            a = a*a;
        }
        return c;
    }
    int main( )
    {
        long long n;
        while(scanf("%lld",&n)!=EOF)
        {
            if(n==-1)
            break;
            Mat A;
            A.mat[0][0]=1;A.mat[0][1]=1;
            A.mat[1][0]=1;A.mat[1][1]=0;
            Mat ans=A^n;
            printf("%lld
    ",ans.mat[0][1]);
        }
        return 0;
    }
    View Code

     

    #include<stdio.h>
    #include<vector>
    using namespace std;
    #define inf 0x3f3f3f3f
    typedef vector<int>vec;
    typedef vector<vec>mat;
    typedef long long ll;
    const int M = 10000;
    ll n;
    mat mul(mat &A , mat &B)
    {
        mat C(A.size() , vec(B.size()));
        for(int i=0 ; i<A.size() ; i++)
        {
            for(int k=0 ; k<B.size() ; k++)
            {    
                if(A[i][k]==0)
                continue;
                for(int j=0 ; j<B[0].size() ; j++)
                {   
                    if(B[k][j]==0)
                    continue;
                    C[i][j] = (C[i][j]+A[i][k]*B[k][j])%M;
                }
            }
        }
        return C;
    }
    
    mat pow(mat A,ll n)
    {
        mat B(A.size(),vec(A.size()));
        for(int i=0 ; i<A.size() ; i++)
            B[i][i]=1;
        while(n>0)
        {
            if(n&1)
            B = mul(B,A);
            A = mul(A,A);
            n  >>= 1;
        }
        return B;
    }
    void so( )
    {
        mat A(2,vec(2));
        A[0][0]=1;A[0][1]=1;
        A[1][0]=1;A[1][1]=0;
        A = pow(A,n);
        printf("%d
    ",A[1][0]);
    
    }
    
    int main()
    {
    
        while(scanf("%lld",&n)!=EOF)
        {
            if(n==-1)
            break;
            so();
        }
    
        return 0;
    }
    View Code

    给一些简单的递推式
    1.f(n)=a*f(n-1)+b*f(n-2)+c;(a,b,c是常数)

    2.f(n)=c^n-f(n-1) ;(c是常数)

    继续例题二:poj3233

    Description

    Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

    Input

    The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

    Output

    Output the elements of S modulo m in the same way as A is given

    这题就是求一个矩阵的和式:S(k),直接对和式建立递推:

    建立矩阵,注意此处S和A都是2*2的矩阵,E表示单位矩阵,O表示零矩阵(全是0,与其他矩阵相乘都为0),显然E,O都是2*2的

    这里转移矩阵是4*4.OVER

     

    一般这种题目都是要找递推式,这是难点.

    例三:HDU2276

    题意就是说n个灯排成环i号灯的左边是i-1号,1的左边是n号,给定初始灯的开闭状态(用1,0表示),然后每一秒都操作都是对于每个灯i,如果它左边的灯是开的就把i状态反转(0变1,1变0),问m秒都最终状态是什么,m<=1e8,n<=100;

    这题的递推式没有明说,但是每一秒操作一次其实就是一次递推,那么关键就是要找转移矩阵,而且比较是常数矩阵,这个才能用等比的性质

    我们把n个灯的状态看出一个F(n),其实就是一个n*1的01矩阵,然后考虑从上一秒的状态F(n-1)如何转移到F(n)。既然每个状态都要转移,总共n个状态,可以看出转移矩阵就是一个n*n的矩阵(每一个灯的状态都用一个1*n的矩阵来转移)

    先考虑第一个灯:影响它新状态的只有它左右的灯和它自己的状态,仔细想一下,1*n的转移中只有这两位置有用,那么其他都是0,就这样

    这里state2到staten-1都被0干掉了,只有第一个和1号左边的灯有效

    (这里不具体讲了,只有0,1的状态,自己枚举一下state1和state2,矩阵相乘后都是能正确转移的)

    其他灯的考虑都是一样的,最终:

    OVER

     

    例四:HDU 5015
    题意就是一个矩阵a,第一行依次是0,233,2333,23333......,给出矩阵的递推式:a[i][j]=a[i-1][j]+a[i][j-1],输入的是第零列,求a[n][m],n ≤ 10,m ≤ 109
    解:n很小,m很大,m大概就是作为幂次,以每一列为一个矩阵A(n-1)并且向下一列A(n)转移,那么关键就是找转移矩阵了:
    先看第一行的数233后面每次都添一个3,显然递推关系就是A(n-1)*10+3=A(n),这里多一个3,必须把列数新增一个元素3,也就是一个(n+1)*1的矩阵A(n).
    然后其他的元素转移会出现一个问题,a[i][j]=a[i-1][j]+a[i][j-1];这里a[i-1][j]依旧是新一列的元素,这是不能矩阵转移的,因为矩阵转移必须从旧的矩阵状态A(n-1)变到新状态A(n)。
    那么这里a[i-1][j]就用递归思想,沿用上一行的转移矩阵(上一行能转移出a[i-1][j]),再加上新增的转移:

     

    OVER

    当然矩阵还有很多奇葩的递推,比如这矩阵优化的dp题

    推荐一些题目:

    简单的:

    http://acm.hdu.edu.cn/showproblem.php?pid=1757

    http://acm.hdu.edu.cn/showproblem.php?pid=1575

    不简单的:

    http://acm.hdu.edu.cn/showproblem.php?pid=3483

    http://acm.hdu.edu.cn/showproblem.php?pid=2855

    http://acm.hdu.edu.cn/showproblem.php?pid=3658

    http://acm.hdu.edu.cn/showproblem.php?pid=4565

    这大牛总结的不错哇

    这个也是哇

  • 相关阅读:
    接口缓存--把接口放在redis数据库中,减少访问量
    使用vue和drf后台进行登录页面和注册页面(本文大概的疏通一下前后台是怎么交互的)
    vue导入css,js和放置html代码
    存储过程
    触发器
    视图
    pymysql模块使用
    权限管理
    多表数据查询
    单表数据查询
  • 原文地址:https://www.cnblogs.com/shuaihui520/p/9147218.html
Copyright © 2020-2023  润新知