zoukankan      html  css  js  c++  java
  • 矩阵快速幂

    矩阵快速幂(入门)

    定义及模板

    对于方阵存在求幂运算的概念。对于单个数值,我们通过将指数进行二进制拆分的思想将求幂运算的复杂度降至log(n),这种思想同样可以应用到方阵的求幂运算中。事实上,方阵的快速幂与单个数值的快速幂的写法完全一致,只需要对于乘法、取模进行运算符重载即可。

    下面给出矩阵快速幂的模板

    #include <cstdio>
    #include <cstdlib>
    
    using namespace std;
    
    #define N 5+1 
    // 最大矩阵规模,根据题目修改
    
    typedef long long ll;
    const ll mod = 1e9+7; // 取模
    
    struct Matrix{
        int x;
        int y;
        ll a[N][N];
        Matrix(){};
        Matrix(int m,int n){
            x = m;
            y = n;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    a[i][j] = 0;
                }
            }
        }
    
        ll* operator[](int x){
            return a[x];
        }
    
        void ones(){
            for(int i = 1; i <= x; i++){
                a[i][i] = 1;
            }
        }
        Matrix operator *(Matrix b){
            Matrix ans(x,b.y);
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= b.y; j++){
                    for(int k = 1; k <= y; k++){
                        ans[i][j] += a[i][k] * b[k][j];
                    }
                    ans[i][j] %= mod;
                }
            }
            return ans;
        }
    
        Matrix operator %(ll md){
            Matrix ans = *this;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    ans[i][j] %= md;
                }
            }
            return ans;
        }
    
        Matrix operator ^(ll p){
            Matrix ans(x,x);
            Matrix base = *this;
            ans.ones();
            while(p){
                if(p&1){
                    ans = ans * base;
                    ans = ans % mod;
                }
                base = base * base;
                base = base % mod;
                p >>= 1;
            }
            return ans;
        }
    };
    
    
    int main(){
    
        return 0;
    }
    

    应用:矩阵加速

    矩阵加速用于线性方程的加速递推

    level 1:简单线性递推

    • e.g1 : 求Fibonacci(n)%mod,mod=1e9+7

    考虑到斐波那契数列的递推公式

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

    我们可以将数列的第n项、第n-1项、第n-2项组合成矩阵形式

    // 数列第n、n-1、n-2的组合
       [f(n)   f(n-1)   f(n-2)]
    // 数列第n+1、n、n-1的组合
       [f(n+1) f(n)     f(n-1)]
    

    这两个矩阵满足关系式子

    [left[ egin{matrix} a_{n}&a_{n-1}&a_{n-2} end{matrix} ight] * left[ egin{matrix} 1 & 1 & 0\ 1 & 0 & 1\ 0 & 0 & 0 end{matrix} ight] = left[ egin{matrix} a_{n+1}&a_{n}&a_{n-1} end{matrix} ight] ]

    进而得出

    [left[ egin{matrix} 2 & 1 & 1 end{matrix} ight] * left[ egin{matrix} 1 & 1 & 0\ 1 & 0 & 1\ 0 & 0 & 0 end{matrix} ight]^{n-2} = left[ egin{matrix} a_{n}&a_{n-1}&a_{n-2} end{matrix} ight] ]

    因此,要求fibonacci(n),通过矩阵加速可以将O(n)降至O(log(n))


    本题关系式为

    [left[ egin{matrix} 2&1&1&1 end{matrix} ight] * left[ egin{matrix} 1&0&0&0\ 0&1&1&0\ 1&0&0&1\ 0&1&0&0 end{matrix} ight]^{n-4} = left[ egin{matrix} a_{n}&a_{n-1}&a_{n-2}&a_{n-3} end{matrix} ight] ]

    #include <cstdio>
    #include <cstdlib>
    
    using namespace std;
    
    #define N 5+1 
    // 最大矩阵规模,根据题目修改
    
    typedef long long ll;
    const ll mod = 1e9+7;
    
    struct Matrix{
        int x;
        int y;
        ll a[N][N];
        Matrix(){};
        Matrix(int m,int n){
            x = m;
            y = n;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    a[i][j] = 0;
                }
            }
        }
    
        ll* operator[](int x){
            return a[x];
        }
    
        void ones(){
            for(int i = 1; i <= x; i++){
                a[i][i] = 1;
            }
        }
        Matrix operator *(Matrix b){
            Matrix ans(x,b.y);
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= b.y; j++){
                    for(int k = 1; k <= y; k++){
                        ans[i][j] += a[i][k] * b[k][j];
                    }
                    ans[i][j] %= mod;
                }
            }
            return ans;
        }
    
        Matrix operator %(ll md){
            Matrix ans = *this;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    ans[i][j] %= md;
                }
            }
            return ans;
        }
    
        Matrix operator ^(ll p){
            Matrix ans(x,x);
            Matrix base = *this;
            ans.ones();
            while(p){
                if(p&1){
                    ans = ans * base;
                    ans = ans % mod;
                }
                base = base * base;
                base = base % mod;
                p >>= 1;
            }
            return ans;
        }
    };
    
    
    int main(){
        int t;
        scanf("%d",&t);
    
        Matrix k(4,4);
        k[1][1] = 1;
        k[2][2] = 1;
        k[2][3] = 1;
        k[3][1] = 1;
        k[3][4] = 1;
        k[4][2] = 1;
    
        Matrix s(1,4);
        s[1][1] = 2;
        s[1][2] = 1;
        s[1][3] = 1;
        s[1][4] = 1;
        while(t--){
            ll n;
            scanf("%lld",&n);
            if(n <= 3){
                printf("1
    ");
            }else{
                Matrix ans = s * (k^(n-4));
                printf("%lld
    ",ans[1][1]);
            }
        }
    
        return 0;
    }
    

    level 2: 带有求和项的简单递推

    由于矩阵加速题难点在于关系矩阵的寻找,而不在于代码实现,为了节省篇幅,之后的题解不再出现具体代码实现。

    e.g1:HDU-3306 Another kind of Fibonacci

    本题的关系矩阵略难推导,应该注意到

    [egin{align} & S_{n+1} = S_{n} + A_{n+1}^2\ & A_{n+2}^2 = (xA_{n+1} + yA_{n})^2 = x^2A_{n+1}^2 + y^2A_n^2 + 2xyA_{n+1}A_n \ & A_{n+2}A_{n+1} = (xA_{n+1}+yA_{n})A_{n+1} = xA_{n+1}^2 + yA_{n+1}A_n end{align} ]

    因此可以得到

    [left[ egin{matrix} S_n & A_{n+1}^2 & A_n^2 & A_{n+1}A_n end{matrix} ight] * left[ egin{matrix} 1 & 0 & 0 & 0\ 1 & x^2 & 1 & x\ 0 & y^2 & 0 & 0\ 0 & 2xy & 0 & y end{matrix} ight] = left[ egin{matrix} S_{n+1}^2 & A_{n+2}^2 & A_{n+1}^2 & A_{n+2}A_{n+1} end{matrix} ight] ]

    本题小结

    1. 对于求和项有 (S_{n+1}-S_n = ...)
    2. 如果递推公式中出现了别的项,无法通过已有项的线性组合得出(如(A_{n+1}A_n)项无法通过(A_{n+1})(A_n)的线性组合得出),那么考虑更改递推公式或在向量中添加新项

    level 3: 打表寻找递推公式

    e.g1:luogu5004 专心OI - 跳房子

    本题不难想到二维DP的解法(虽然会爆)

    伪代码如下,其中i表示在格子上停留的次数,j表示当前走在第j个格子上

    procedure(DP):
    int dp[N][N]
    int n,m
    read n,m
    int maxstep = (n+1)/(m+1)
    
    for i:= 1 to n
        dp[1][i] = 1
    
    for i:= 2 to n
        for j:= 1 to n
            int begin = 1
            int end = j-1+m
            for k:= begin to end
                dp[i][j] += dp[i-1][k]
    
    ans = 1
    for i:= 1 to maxstep
        for j:= 1 to n
            ans += dp[i][j]
    
    return ans
    

    用这个思路来运行较小的n,m值,打表寻找规律:

    • m=1
    n = 0   m = 1   ans = 1
    n = 1   m = 1   ans = 2
    n = 2   m = 1   ans = 3
    n = 3   m = 1   ans = 5
    n = 4   m = 1   ans = 8
    n = 5   m = 1   ans = 13
    n = 6   m = 1   ans = 21
    n = 7   m = 1   ans = 34
    n = 8   m = 1   ans = 55
    n = 9   m = 1   ans = 89
    n = 10  m = 1   ans = 144
    n = 11  m = 1   ans = 233
    n = 12  m = 1   ans = 377
    n = 13  m = 1   ans = 610
    n = 14  m = 1   ans = 987
    n = 15  m = 1   ans = 1597
    n = 16  m = 1   ans = 2584
    n = 17  m = 1   ans = 4181
    n = 18  m = 1   ans = 6765
    n = 19  m = 1   ans = 10946
    n = 20  m = 1   ans = 17711
    

    可以发现m=1时ans符合斐波那契数列的增长规律

    [m = 1 o left{ egin{align} & A_{0} = 1\ & A_{1} = 2\ & A_{2} = 3\ & A_{n} = A_{n-1}+A_{n-2} (n geqslant 3) end{align} ight. ]

    • m = 2
    n = 0   m = 2   ans = 1
    n = 1   m = 2   ans = 2
    n = 2   m = 2   ans = 3
    n = 3   m = 2   ans = 4
    n = 4   m = 2   ans = 6
    n = 5   m = 2   ans = 9
    n = 6   m = 2   ans = 13
    n = 7   m = 2   ans = 19
    n = 8   m = 2   ans = 28
    n = 9   m = 2   ans = 41
    n = 10  m = 2   ans = 60
    n = 11  m = 2   ans = 88
    n = 12  m = 2   ans = 129
    n = 13  m = 2   ans = 189
    n = 14  m = 2   ans = 277
    n = 15  m = 2   ans = 406
    n = 16  m = 2   ans = 595
    n = 17  m = 2   ans = 872
    n = 18  m = 2   ans = 1278
    n = 19  m = 2   ans = 1873
    n = 20  m = 2   ans = 2745
    

    [m = 2 o left{ egin{align} & A_{0} = 1\ & A_{1} = 2\ & A_{2} = 3\ & A_{3} = 4\ & A_{n} = A_{n-1}+A_{n-3} (n geqslant 4) end{align} ight. ]

    • m = 4
    n = 0   m = 4   ans = 1
    n = 1   m = 4   ans = 2
    n = 2   m = 4   ans = 3
    n = 3   m = 4   ans = 4
    n = 4   m = 4   ans = 5
    n = 5   m = 4   ans = 6
    n = 6   m = 4   ans = 8
    n = 7   m = 4   ans = 11
    n = 8   m = 4   ans = 15
    n = 9   m = 4   ans = 20
    n = 10  m = 4   ans = 26
    n = 11  m = 4   ans = 34
    n = 12  m = 4   ans = 45
    n = 13  m = 4   ans = 60
    n = 14  m = 4   ans = 80
    n = 15  m = 4   ans = 106
    n = 16  m = 4   ans = 140
    n = 17  m = 4   ans = 185
    n = 18  m = 4   ans = 245
    n = 19  m = 4   ans = 325
    n = 20  m = 4   ans = 431
    

    [m = 4 o left{ egin{align} & A_{0} = 1\ & A_{1} = 2\ & A_{2} = 3\ & A_{3} = 4\ & A_{4} = 5\ & A_{5} = 6\ & A_{n} = A_{n-1}+A_{n-5} (n geqslant 6) end{align} ight. ]

    找到规律

    [forall m o left{ egin{align} & A_n = n+1 (0 leqslant n leqslant m+1)\ & A_{n} = A_{n-1}+A_{n-m-1} (n geqslant m+2) end{align} ight. ]

    因此有

    本题关系式为

    [left[ egin{matrix} A_n & A_{n-1} &cdots & A_{n-m} end{matrix} ight]* left[ egin{matrix} 1 & 1 & 0 & 0 &cdots & 0 & 0\ 0 & 0 & 1 & 0 &cdots & 0 & 0\ 0 & 0 & 0 & 1 & cdots & 0&0\ 0 & 0 & 0 & 0 & 1 & cdots &0\ vdots & vdots & vdots & vdots & vdots & ddots & 1& \ 1 & 0 & 0 & 0 & cdots & 0 & 0\ end{matrix} ight] = left[ egin{matrix} A_{n+1} & A_{n} cdots & A_{n-m+1} end{matrix} ight] ]

    [R_{1 imes (m+1)} * R_{(m+1) imes(m+1)} = R_{1 imes (m+1)} ]

    代码

    #include <cstdio>
    #include <cstdlib>
    using namespace std;
    
    #define N 18+1 
    // 最大矩阵规模,根据题目修改
    
    typedef long long ll;
    const ll mod = 1e9+7; // 取模
    
    struct Matrix{
        int x;
        int y;
        ll a[N][N];
        Matrix(){};
        Matrix(int m,int n){
            x = m;
            y = n;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    a[i][j] = 0;
                }
            }
        }
    
        ll* operator[](int x){
            return a[x];
        }
    
        void ones(){
            for(int i = 1; i <= x; i++){
                a[i][i] = 1;
            }
        }
        Matrix operator *(Matrix b){
            Matrix ans(x,b.y);
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= b.y; j++){
                    for(int k = 1; k <= y; k++){
                        ans[i][j] += a[i][k] * b[k][j];
                    }
                    ans[i][j] %= mod;
                }
            }
            return ans;
        }
    
        Matrix operator %(ll md){
            Matrix ans = *this;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    ans[i][j] %= md;
                }
            }
            return ans;
        }
    
        Matrix operator ^(ll p){
            Matrix ans(x,x);
            Matrix base = *this;
            ans.ones();
            while(p){
                if(p&1){
                    ans = ans * base;
                    ans = ans % mod;
                }
                base = base * base;
                base = base % mod;
                p >>= 1;
            }
            return ans;
        }
    };
    
    int dp[N][N];
    int main(){
        
        ll n,m;
        scanf("%lld%lld",&n,&m);
        Matrix a(1,m+1);
        a[1][m+1] = 1;
        for(int i = m; i >= 1; i--){
            a[1][i] = a[1][i+1]+1;
        }
    
        Matrix b(m+1,m+1);
        b[1][1] = 1;
        b[m+1][1] = 1;
        for(int i = 1; i <= m; i++){
            b[i][i+1] = 1;
        }
    
        if(n < m){
            printf("%lld
    ",a[1][m+1-n]);
        }else{
            Matrix ans = a * (b^(n-m));
            printf("%lld
    ",ans[1][1]);
        }
        // system("pause");
        return 0;
    }
    
    

    level4: 数论+DP+矩阵加速

    e.g1: P3746 [六省联考2017]组合数问题

    本题如果具备组合数递推的前置知识就不算难题,关于组合数的前置知识请看:浅谈组合数相关性质

    根据组合数的递推性质

    [C_{m}^{n} = C_{m-1}^{n} + C_{m-1}^{n-1} ]

    它其实就是DP方程,表示前m个物品选择n个物品的方案数。根据最后一个物品选或不选分别对应出(C_{m-1}^n)(C_{m-1}^{n-1})

    写成DP形式即为

    [dp[i][j] = dp[i-1][j] + dp[i-1][j-1] ]

    本题要求的为

    [sum_{i=0}^{+infin}C_{nk}^{nk+r} ]

    它所表示的是从nk个物品中选择m个物品,其中m%k = r,求对于所有的m的方案数的总和

    我们设(dp[i][j])表示从前i个物品中选择出 m % k = j 个物品的方案数

    [dp[i][0] = dp[i-1][0] + dp[i-1][k-1] \ dp[i][j] = dp[i-1][j] + dp[i-1][j-1] ]

    因此递推关系式为:

    [left[ egin{matrix} 1 & 0 & 0 & 0 & cdots & cdots & 0 end{matrix} ight] * left[ egin{matrix} 1 & 1 & 0 & 0 & 0 & cdots & 0 \ 0 & 1 & 1 & 0 & 0 & cdots & 0 \ 0 & 0 & 1 & 1 & 0 & cdots & 0 \ 0 & 0 & 0 & 1 & 1 & cdots & 0 \ vdots & vdots&vdots& vdots&vdots &cdots & 0 \ vdots & vdots&vdots& vdots&vdots &1 & 0 \ vdots & vdots&vdots& vdots&vdots &1 & 1 \ 1 & 0 & 0 & 0 & 0 & cdots & 1 end{matrix} ight]^{nk} = left[ egin{matrix} f_{nk,0} & f_{nk,1} & f_{nk,2} & cdots & cdots & f_{nk,k-1} end{matrix} ight] ]

    最后,k=1需要特判,k=1时应该为

    [left[ egin{matrix} 1 end{matrix} ight] * left[ egin{matrix} 2 end{matrix} ight]^{n} = left[ egin{matrix} f_{n,0} end{matrix} ight] ]

    最终答案为

    ans[1][r+1]; // ans: Matrix(1,k)
    
    #include <cstdio>
    #include <cstdlib>
    
    using namespace std;
    
    #define N 50+5 
    // 最大矩阵规模,根据题目修改
    
    typedef long long ll;
    ll mod; // 取模
    
    struct Matrix{
        int x;
        int y;
        ll a[N][N];
        Matrix(){};
        Matrix(int m,int n){
            x = m;
            y = n;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    a[i][j] = 0;
                }
            }
        }
    
        ll* operator[](int x){
            return a[x];
        }
    
        void ones(){
            for(int i = 1; i <= x; i++){
                a[i][i] = 1;
            }
        }
        Matrix operator *(Matrix b){
            Matrix ans(x,b.y);
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= b.y; j++){
                    for(int k = 1; k <= y; k++){
                        ans[i][j] += a[i][k] * b[k][j] % mod;
                        ans[i][j] %= mod;
                    }
                }
            }
            return ans;
        }
    
        Matrix operator %(ll md){
            Matrix ans = *this;
            for(int i = 1; i <= x; i++){
                for(int j = 1; j <= y; j++){
                    ans[i][j] %= md;
                }
            }
            return ans;
        }
    
        Matrix operator ^(ll p){
            Matrix ans(x,x);
            Matrix base = *this;
            base = base % mod;
            ans.ones();
            while(p){
                if(p&1){
                    ans = ans * base;
                    ans = ans % mod;
                }
                base = base * base;
                base = base % mod;
                p >>= 1;
            }
            return ans;
        }
    };
    
    int main(){
        ll n,p,k,r;
        scanf("%lld%lld%lld%lld",&n,&p,&k,&r);
        mod = p;
        Matrix a(1,k); // j from 0 to k-1
        a[1][1] = 1; // dp[0][0] = 1
        
        Matrix b(k,k);
        for(int i = 1; i <= k; i++){
            b[i][i] = 1;
            b[i-1][i] = 1;
        } 
        b[k][1] = 1;
        if(k == 1){
            b[1][1] = 2;
        }
        Matrix ans = a * (b^(n*k));
        printf("%lld
    ",ans[1][r+1]);
        // system("pause");
        return 0;
    }
    
    
    ---- suffer now and live the rest of your life as a champion ----
  • 相关阅读:
    连接数据库代码
    C/C++ Basicsfunction pointer
    MSMQ&Com+ Service: How to create an Com+ Service in .NetFramework
    C/C++ Basics>about #define, const
    EndpointAddress
    MSMQ Basics Transactional Messages Processing
    Thread Basics(thread synchronization&Asynchronization) part two
    Thread Basics using Timer to trigger Event at a specified internals
    FAQ about AJAXpart II
    FAQ about AJAXpart I
  • 原文地址:https://www.cnblogs.com/popodynasty/p/13798951.html
Copyright © 2011-2022 走看看