高阶矩阵快速幂

最新推荐文章于 2023-11-02 22:51:50 发布

夜拾柒

最新推荐文章于 2023-11-02 22:51:50 发布

阅读量1.3k

点赞数 1

分类专栏：算法理解数论(math algorithm)

本文链接：https://blog.csdn.net/yeweiyang16/article/details/71512174

版权

算法理解同时被 2 个专栏收录

7 篇文章

订阅专栏

数论(math algorithm)

6 篇文章

订阅专栏

前言

所谓的高阶矩阵快速幂就是在快速幂的基础上，结合矩阵运算计算矩阵的高次方的。它可以将朴素的o(n)的时间复杂度降为o(logn)。

我们先对原理进行简单的介绍，对于快速幂来说，如果要求x^n，一般我们会将幂n 表示成 n = 2^k，那么x^n = ((x^2)^2)......，对于每个数我们只需要做k次平方运算就行了，因此我们只需要将n表示成2的幂次和。

n = 2^k1 + 2^k2 + 2^k3 + ......，那么n就可以表示成x ^ n = x^(2^k1) * x^(2^k2) * x^(2^k3)......,依次求幂就行，这里给出快速幂的算法：

typedef long long ll;

ll mod_pow(ll x, ll n, ll mod)
{
	ll res = 1;
	while(n > 0){
		if(n & 1)
			res = res * x %mod;
			//如果二进制的最低位为1，就要乘上x^(2^i)
		x = x * x % mod;
		n >>= 1;
	}
	return res;
}

有了上面的知识，我们如何快速求的矩阵的N次幂呢？

矩阵快速幂

我们可以将矩阵的幂进行离散化，进行二进制表示例如求A^156 ,我们可以将156(10) = 10011100(2)，也就是A^(156) = (A^4)*(A^8)*(A^16)*(A^128)

就可以将代码改成

while(N)
{
	if(N&1)
		res = res*A;
	n >>= 1;
	A = A*A;
}

例子

斐波那契数列是由下面的递推定义的数列

F0 = 0

F1 = 1

Fn+2 = Fn+1 + Fn

求这个数列第n项的值对10^4取余后的结果

数据范围：0 <= n <= 10^16

分析：通过逐项求解这个递推式，我们可以在O(n)的时间内求解，但是效率太低。当然我们也可以求出通项公式：

来求解，但是因为其中含有无理数，所以我们不能进行求模运算。在这种情况下，我们就要使用矩阵快速幂了

首先，我们需要将递推公式表示成矩阵形式

记这个矩阵为A，那么

只要求出A^n的值，我们就可以求出Fn了，就可以在O(logn)的时间内求出第n项的值

标答如下：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <vector>

using namespace std;

typedef vector<int> vec;
typedef vector<vec> mat;
typedef long long ll;

const int M = 10000;

//计算两个矩阵的积
mat mul(mat &A, mat &B)
{
    mat C(A.size(), vec(B[0].size()));
    for(int i=0;i< A.size(); i++)
    {
        for(int k=0;k<B.size();k++)
        {
            for(int j=0;j<B[0].size();j++)
            {
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % M;
            }
        }
    }
    return C;
}

//计算A^n
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;
}

ll n;

int main()
{
    scanf("%lld",&n);
    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\n", A[1][0]);
    return 0;
}