快速幂矩阵快速幂

Twilight Sparkle.

已于 2023-02-25 19:08:09 修改

阅读量6.6k

点赞数 29

分类专栏： ACM学习文章标签：矩阵线性代数 acm竞赛

于 2021-10-22 12:29:51 首次发布

本文链接：https://blog.csdn.net/qq_52466006/article/details/120903225

版权

ACM学习专栏收录该内容

10 篇文章 8 订阅

订阅专栏

前言：好像没啥好写的，链接可能还没有更新完

快速幂

快速幂，用于解决 $a^n$ 当n很大时的情况。通常与取模同时应用。

用最笨的方法求 $a^n$ ,即 $a*a*a*a*....$ 。时间复杂度为 $O(n)$ ,而快速幂（附带取模），可以将时间复杂度降低为 $O(logn)$ 。

利用倍增的思想，例如 $3*3*3*3*3*3*3*3$ ，等于 $(3*3)*(3*3)*(3*3)*(3*3)$ ，又等于 $((3*3)*(3*3))*((3*3)*(3*3))$ ,即 $3^8 = 9^4 = 81^2$ ,如果n为奇数， $a^n = a*a^{[n/2]}*a^{[n/2]}$ ,那么换成代码就是：

ll Powermod(ll a,ll b)
{
	ll ans=1;
	a=a%mod;
	while(b>0)
	{
		if(b%2==1)  ans=(ans*a)%mod;
		b=b/2;
		a=(a*a)%mod;
	}
	return ans;
}

其中的取模操作保证对最终答案取模结果不变，且防止超出数据类型范围。

矩阵快速幂

前置知识：矩阵乘法

设A为 $m*p$ 的矩阵，B为 $p*n$ 的矩阵，那么称 $m*n$ 的矩阵C为矩阵A与B的乘积，记作
$C = A*B$ ，其中矩阵C中的第 $i$ 行第 $j$ 列元素可以表示为：

如下所示：

（来源：搜狗百科）

矩阵乘法时间复杂度：

一个 $m*n$ 的矩阵和一个 $n*p$ 的矩阵相乘，将会得到一个 $m*p$ 的矩阵,

对于矩阵 $A(n*m)$ , $B(n*p)$ , $A$ 表示n行m列的矩阵， $A*B$ 的时间复杂度是 $O(n*m*p)$ ,代码对应如下：

//矩阵乘法
void MXMP(int m,int n,int p)//Matrix multiplication
{
    for(int i = 1;i<=m;i++)
        for(int j = 1;j<=p;j++)
            for(int k = 1;k<=m;k++)
                C[i][j] = (C[i][j] + (A[i][k]*B[k][j])%mod)%mod;//防止超界
}

矩阵快速幂

矩阵快速幂，即给定一个矩阵 $A(m*n)$ ),快速计算 $A^n$ 。一般来说，矩阵快速幂只会涉及方阵即 $A(n*n)$ ,所以这里以方阵为例。

自己写一个矩阵乘法，然后用快速幂的套进去就是矩阵快速幂。

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<queue>
#include<vector>
#include<stack>
using namespace std;

typedef long long int ll;
const int mod = (int)1e9+7;
const int MAX_N = 1e3;
int A[MAX_N][MAX_N],res[MAX_N][MAX_N];
int temp[MAX_N][MAX_N];

//A = A*B
void MXMP(int a[][MAX_N],int b[][MAX_N],int n)
{
    //重置临时矩阵temp
    for(int i = 1;i<=n;i++)
        for(int j = 1;j<=n;j++)
            temp[i][j] = 0;

    for(int i = 1;i<=n;i++)
        for(int j = 1;j<=n;j++)
            for(int k = 1;k<=n;k++)
                temp[i][j] = (temp[i][j] + (a[i][k]*b[k][j])%mod)%mod;//防止超界
    
    for(int i = 1;i<=n;i++)
        for(int j = 1;j<=n;j++)
            a[i][j] = temp[i][j];
}

//A = A**x
void PowerMod(int A[][MAX_N],int n,int x)//x为次幂，n为矩阵行/列
{
    memset(res,0,sizeof(res));
    for(int i = 1;i<=n;i++)  res[i][i] = 1;//初始化为单位矩阵
    while(x){
        if(x&1)  MXMP(res,A,n);
        MXMP(A,A,n);
        x >>= 1;
    }
    for(int i = 1;i<=n;i++)
        for(int j = 1;j<=n;j++)
            A[i][j] = res[i][j];
}


int main()
{
    //读入A矩阵
    int n,x;
    cin>>n>>x;//n行n列，x为次幂
    for(int i = 1;i<=n;i++)
        for(int j = 1;j<=n;j++)
            cin>>A[i][j];
    //_______________________________
    PowerMod(A,n,x);
    //检查A矩阵
    for(int i = 1;i<=n;i++){
        for(int j = 1;j<=n;j++){
            cout<<A[i][j]<<" ";
        }
        cout<<endl;
    }
    //________________________________


    return 0;
}