矩阵快速幂_矩阵的快速幂-CSDN博客

本文链接：https://blog.csdn.net/qq_35527032/article/details/79509447

快速幂

快速幂算法很熟悉很容易理解，原理 $n^m = \begin{cases} (n^2)^{m/2} & m为偶数 \\(n^2)^{m/2}\cdot n & m为奇数\end{cases}$
二分可以将时间复杂度由 $O(n)$ 降到 $O(\log n)$ 。
通常m的范围非常大，常常到 $10^9$ ，所以一般给出的题中会要求取模。视情况使用32位或者64位整型。

$n^m\%d$ 的实现(见快速幂取模)：

int Quick_Power_Mod(int n,int m,int d)
{
    int ans = 1;
    n = n%d;
    while(m > 0)
    {
        if(m%2 == 1)
        {
            ans = (ans * n) % d;
        }
        m = m/2;
        n = (n*n)%d;
    }
    return ans;
}

矩阵快速幂

如果将其中的n换为矩阵的话，则为矩阵快速幂。
例题：XDOJ-1027: Feibonaqi数列
1.f(0)=0,f(1)=1
2.f(n)=2f(n-1)+f(n-2),n>1
如果直接递归的话时间会爆，设数组保存的话空间会爆。只能快速幂。
求出递推公式也不可行，因为特征根是无理数。
所以需要找出矩阵形式的递推关系，从而求出矩阵形式的显示公式：
$\left[ \begin{matrix} f(n) \\f(n-1) \end{matrix} \right] =$ $\left[\begin{matrix}2 & 1 \\1 & 0\end{matrix}\right ]$ $\times \left[ \begin{matrix} f(n-1) \\f(n-2) \end{matrix} \right]$

递推得到： $\left[ \begin{matrix} f(n) \\f(n-1) \end{matrix} \right] =$ $\left[\begin{matrix}2 & 1 \\1 & 0\end{matrix}\right ]^{n-1}$ $\times \left[ \begin{matrix} f(1) \\f(0) \end{matrix} \right]$

则使用快速幂即可解决：

//XDOJ-1027----201.3.9
//矩阵快速幂 
#include<iostream>
#include<cstdlib>
typedef long long LL;
const long long Max=1000000007;
using namespace std;
void Matrix_mul(LL A[2][2],LL B[2][2]); //矩阵相乘 
LL Matrix_pow_mod(int n);    

int main(void)
{
    int n;
    while(cin >> n)
    {
        if(n == 0 || n == 1)
        {
            cout << n << "\n";
        }
        else
        {
            cout << Matrix_pow_mod(n-1) << "\n";
        }
    }
    return 0;
}
//2*2矩阵想乘 
void Matrix_mul(LL A[2][2],LL B[2][2])
{
    LL ans[2][2] = {{0,0},{0,0}};
    for(int i = 0; i < 2 ; i++)
    {
        for(int j= 0; j < 2; j++)
        {
            for(int k=0; k < 2; k++)
            {
                ans[i][j] += A[i][k] * B[k][j];
            }
            if(ans[i][j] > Max) //大于时才取模 
            {
                ans[i][j] %= Max;
            }
        }
    }
    for(int i=0; i < 2; i++)
    {
        for(int j=0; j < 2; j++)
        {
            B[i][j] = ans[i][j];
        }   
    }   
}

//矩阵快速幂 
LL Matrix_pow_mod(int n)
{
    LL A[2][2] = {{2,1},{1,0}};
    LL Ans[2][2] = {{1,0},{0,1}};   //二阶单位阵
    while(n > 0)
    {
        if(n%2 == 1) // n & 1
        {
            Matrix_mul(A,Ans);
        }
        Matrix_mul(A,A);
        n = n/2; //n >>= 2;
    }
    return Ans[0][0]; 
}

类似题：poj-3070-Fibonacci

构造递推关系

上面的递推关系很简单，这类问题的重点就是构造矩阵递推关系。

如：POJ 3233 Matrix Power Series

考虑到 $S_k = \sum_{i=1}^kA^i = E\cdot S_{k-1}+A^k$ 或者 $S_k = A\cdot S_{k-1}+A$
可以构造出：
$\left[ \begin{matrix} S_k \\A^k \end{matrix} \right] =$ $\left[\begin{matrix}E & A \\0 & A\end{matrix}\right ]$ $\times \left[ \begin{matrix} S_{k-1} \\A^{k-1} \end{matrix} \right]$ 或者 $\left[ \begin{matrix} S_k \\E \end{matrix} \right] =$ $\left[\begin{matrix}A & A \\0 & E\end{matrix}\right ]$ $\times \left[ \begin{matrix} S_{k-1} \\E \end{matrix} \right]$

方法不一。其中A是n阶方阵，E是n阶单位阵，0是n阶零矩阵。

比如由前者就可以得出： $\left[ \begin{matrix} S_k \\A^k \end{matrix} \right] =$ $\left[\begin{matrix}E & A \\0 & A\end{matrix}\right ]^{n-1}$ $\times \left[ \begin{matrix} A \\A \end{matrix} \right]$ 或者 $\left[\begin{matrix}E & A \\0 & A\end{matrix}\right ]^n$ $=\left[\begin{matrix}E & S_k \\0 & A^k\end{matrix}\right ]$

实现：

//poj-3233----矩阵快速幂
//2018,3,10
#include<iostream>
#include<cstring>
using namespace std;
const int N = 65;

void Matrix_mul(int A[N][N],int B[N][N],int n,int m);
void Matrix_pow_mod(int n,int k,int m,int A[N][N]);
void print_matrix(int A[N][N],int n);
void input_matrix(int A[N][N],int n);

int main(void)
{
    int n,k,m;
    int A[N][N];
    while(cin >> n >> k >> m)
    {
        input_matrix(A,n);
        Matrix_pow_mod(n,k,m,A);
    }
    return 0;
}

void Matrix_mul(int A[N][N],int B[N][N],int n,int m)
{
    int res[N][N];
    memset(res,0,sizeof(res));
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j <= n; j ++)
        {
            for(int k = 1; k <= n; k ++)
            {
                if(A[i][k] || B[k][j])
                {
                    res[i][j] += A[i][k] * B[k][j];
                    res[i][j] %= m;
                }
            }   
        }
    }
    for(int i = 0; i <= n; i ++)
    {
        for(int j = 0; j <= n; j ++)
        {
            B[i][j] = res[i][j];
        }
    }
}

void Matrix_pow_mod(int n,int k,int m,int A[N][N])
{
    int B[N][N];
    int ans[N][N];
    int C[N][N];
    //矩阵初始化 
    memset(B,0,sizeof(B));
    memset(ans,0,sizeof(ans));
    memset(C,0,sizeof(C));
    for(int i = 1; i <= n; i ++)
    {
        B[i][i] = 1;
    }
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j <= n; j ++)
        {
            B[i][n+j] = A[i][j];
            B[n+i][n+j] = A[i][j];
        }
    }
    for(int i=0; i <= 2*n; i++)
    {
        ans[i][i] = 1;
    }
    //求k次幂 
    while(k > 0)
    {
        if(k & 1)
        {
            Matrix_mul(B,ans,2*n,m);
        }
        k >>= 1;
        Matrix_mul(B,B,2*n,m);
    }
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j <= n; j ++)
        {
            C[i][j] = ans[i][n+j];
        }
    }
    //print the result 
    print_matrix(C,n);
}

void print_matrix(int A[N][N],int n)
{
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j < n; j ++)
        {
            cout << A[i][j] << " ";
        }
        cout << A[i][n] << endl;
    }
}

void input_matrix(int A[N][N],int n)
{
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j <= n; j ++)
        {
            cin >> A[i][j];
        }
    }
}

优化：很明显可以看到，2n阶方阵中左半部分在乘幂中是没有变化的，所以可以优化掉大概一半的时间（这里没有实现）。一般情况下对于其中0较多的矩阵，可以先判断是否为0再乘。

注意这里的矩阵乘法采用的是一般的算法，复杂度 $O(n^3)$ ，矩阵乘法有一个 $O(n^{log7})$ 的算法: Strassen算法。采用分治策略，比较复杂暂时不会不做实现。

综上：矩阵快速幂思路与快速幂完全一致，重点是找到递推关系。

对于一些简单的递推式：

1. $f(n) = af(n-1) +bf(n-2)+c$

$\left( \begin{matrix} f(n) \\f(n-1) \\c \end{matrix} \right) =$ $\left(\begin{matrix}a & b & 1\\1 & 0 & 0 \\0 & 0 & 1\end{matrix}\right)$ $\times \left( \begin{matrix} f(n-1) \\f(n-2) \\c\end{matrix} \right)$

2. $f(n) = ac^n + bf(n-1) + d$

$\left( \begin{matrix} f(n) \\c^n \\d \end{matrix} \right) =$ $\left(\begin{matrix}b & ac & 1\\0 & c & 0 \\0 & 0 & 1\end{matrix}\right)$ $\times \left( \begin{matrix} f(n-1) \\c^{n-1} \\d\end{matrix} \right)$

类似简单矩阵快速幂模板题：
HOJ 175、HOJ 1575

例

不算那么简单的递推+矩阵快速幂：HOJ 3483 - A Very Simple Problem

题意：求

(\sum k = 1 N k x \cdot x k) mod M

$(\sum_{k=1}^Nk^x\cdot x^k)\mod M$
数据范围：

1≤N,M≤2∗109,and1≤x≤50 1 ≤ N , M ≤ 2 ∗ 10 9 , a n d 1 ≤ x ≤ 50 $1 ≤ N, M ≤ 2*10^9, and 1 ≤ x ≤ 50$

建立递推关系： $S_{n+1} = S_{n} + (n+1)^x\cdot x^{n+1}$
将 $(n+1)^x$ 二项式展开为 $\sum_{k=0}^{x}C(x,k)n^k$ ，则可以构造出如下矩阵递推式，然后使用矩阵快速幂。

|1 xC(x,0) xC(x,1) xC(x,2) ... xC(x,x)|  |Sn       | |S(n+1)           | 
|0 xC(0,0) 0       0       ... 0      |  |x^n * n^0| |x^(n+1) * (n+1)^0| 
|0 xC(1,0) xC(1,1) 0       ... 0      | *|x^n * n^1|=|x^(n+1) * (n+1)^1| 
|0 xC(2,0) xC(2,1) xC(2,2) ... 0      |  |x^n * n^2| |x^(n+1) * (n+1)^2| 
|...                                  |  |...      | |...              | 
|0 xC(x,0) xC(x,1) xC(x,2) ... xC(x,x)|  |x^n * n^x| |x^(n+1) * (n+1)^x|
//还能这样构造，我是真的佩服，只能说我太菜吧！
//出处：http://blog.csdn.net/winycg/article/details/52982723

实现：

//hoj--3483--矩阵快速幂 + 二项式定理 
//2018.3.12
#include<iostream>
#include<cstring>
using namespace std;
typedef long long LL;
const int N = 55;

struct Matrix
{
    LL mat[N][N];   
};

Matrix Matrix_mul(Matrix A,Matrix B, int n,int m); 
Matrix Matrix_pow_mod(Matrix A,int n,int k,int m);
LL C[N][N];
void Combination_num();//use Pascal identity to canculate Combination number

int main(void)
{
    Combination_num();
    int n,x,m;
    while(cin >> n >> x >> m && n != -1)
    {
        Matrix S;
        memset(S.mat,0,sizeof(S.mat));
        S.mat[1][1] = 1;
        for(int i = 2; i <= x+2; i ++)
        {
            S.mat[1][i] = (x * C[x][i-2]) % m;
        }
        for(int i = 2; i <= x+2; i ++)
        {
            for(int j = 2; j <= i; j ++)
            {
                S.mat[i][j] = (x * C[i-2][j-2]) % m;
            }
        }
        LL sum = 0;
        Matrix ans = Matrix_pow_mod(S,x+2,n-1,m);
        for(int i = 1; i <= x+2; i ++)
        {
            sum = (sum + x * ans.mat[1][i]) % m;
        }
        cout << sum << endl;
    }


    return 0;
}


Matrix Matrix_mul(Matrix A,Matrix B, int n,int m)
{
    Matrix res;
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j <= n; j ++)
        {
            res.mat[i][j] = 0;//init the matrix 
            for(int k = 1; k <= n; k ++)
            {
                res.mat[i][j] = (res.mat[i][j] + A.mat[i][k] * B.mat[k][j]) % m;
            }
        }
    }
    return res;
}

Matrix Matrix_pow_mod(Matrix A,int n,int k,int m)
{
    Matrix res;
    //initialization 
    for(int i = 1; i <= n; i ++)
    {
        for(int j = 1; j <= n; j ++)
        {
            res.mat[i][j] = 0;
        }
        res.mat[i][i] = 1;
    }
    while(k > 0)
    {
        if(k & 1)
        {
            res = Matrix_mul(A,res,n,m);
        }
        k >>= 1;
        A = Matrix_mul(A,A,n,m);
    }
    return res;
}

void Combination_num()
{
    C[0][0] = C[1][0] = C[1][1] = 1;
    for(int i = 2; i < N; i ++)
    {
        C[i][0] = C[i][i] = 1;
        for(int j = 1; j < i; j ++)
            C[i][j] = C[i-1][j] + C[i-1][j-1];//of course it's Pascal identity
    }
}