矩阵快速幂以及其优化【华东交大课程】

于 2022-03-20 00:12:40 发布
阅读量420 收藏 2
点赞数 2
分类专栏： 算法入门 文章标签： c++ 算法 矩阵
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/kicamon/article/details/123606049
版权

矩阵快速幂以及其优化【华东交大课程】

快速幂基础：C++快速幂_Kicamon的博客-CSDN博客

矩阵快速幂就是在快速幂的基础上结合矩阵运算的用法，其用途较为广泛，可以很大程度上优化代码。

一、矩阵快速幂的实现

在矩阵快速幂中，我们利用类来存储矩阵及其运算法则。

我们先来看看代码的具体实现，再来讲解其语法。

struct Matrix
{
    ll a[N][N];
    Matrix()//构造函数
    {
        memset(a,0,sizeof a);
    }
    
    void init()//单元矩阵
    {
        for(int i = 0;i <= n;i++)
            a[i][i] = 1;
    }
    
    Matrix operator*(const Matrix& B)const//运算符*的重构
    {
        Matrix C;
        for(int i = 0;i < n;i++)
            for(int k = 0;k < n;k++)
                for(int j = 0;j < n;j++)
                    C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
        return C;
    }
    
    Matrix operator^(const int& t)const
    {
        Matrix A = (*this),res;
        res.init();
        int p = t;
        while(p)
        {
            if(p & 1)res = res * A;
            A = A * A;
            p >>= 1;
        }
        return res;
    }
};

  • 该结构体使用到的陌生语法有：构造函数、运算符重载和this指针，下边我们来看看着两个语法：

    • 构造函数：

      构造函数的特点是在结构体或者类中，函数名与结构体名称相同且其声明前方没有返回值。它会在结构体创建的时候自动调用，完成对结构体的初始化。与之对应的是析构函数，作用是对结构体的清除，感兴趣的同学可以自行了解。

    • 运算符重载：

      • 运算发重载的含义：指的是将已有的运算符赋予新的功能，并通过识别使用是的数据类型和数目来判断采用哪种操作。c/c++中有许多运算符其实已经经过了重构，例如*，在将其作用于地址的时候，它的作用是解引用；将其用于两个数字的时候，它的作用是得到两数的乘积。运算符重载的存在，使得代码可以变得更加简洁易懂。
      • 在这里，我们将 * 运算符重载出求两矩阵的积的作用，将 ^ 运算符重载出求矩阵的次方的作用。当然，在 ^ 中我们使用了重载之后的 * 运算符。

    • this指针

      • this指针存在于成员函数中，是一个用来指向对象本身地址的指针。例如上面代码中，this指针将该指针所在的结构体本身赋值给了结构体A。

二、例题讲解

1、利用矩阵快速幂推导斐波那契数列

在这之前，我们推导斐波那契数列的方法是这样的：

f[1] = f[2] = 1;
for(int i = 3;i <= n;i++)
    f[i] = f[i - 1] + f[i - 2];

但是，学了矩阵快速幂之后，我们可以利用另一种方法来推导它：

在这里插入图片描述
所以我们可以利用矩阵快速幂先得到状态转换矩阵的n次方，再将其于初始矩阵相乘

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 2, MOD = 1e9 + 7;

struct Matrix
{
    ll a[2][2];
    Matrix()
    {
        memset(a, 0, sizeof a);
    }

    void init()
    {
        for (int i = 0; i < 2; i++)
            a[i][i] = 1;
    }

    Matrix operator*(const Matrix& B)const
    {
        Matrix C;
        for (int i = 0; i < 2; i++)
            for (int k = 0; k < 2; k++)
                for (int j = 0; j < 2; j++)
                    C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
        return C;
    }

    Matrix operator^(const int& t)const
    {
        Matrix A = (*this), res;
        res.init();
        int p = t;
        while (p)
        {
            if (p & 1)res = res * A;
            A = A * A;
            p >>= 1;
        }
        return res;
    }
};

Matrix bace;

int main()
{
    bace.a[0][0] = 0, bace.a[0][1] = 1;
    bace.a[1][0] = 1, bace.a[1][1] = 1;
    int n;
    cin >> n;
    if (n == 1 || n == 2)puts("1");
    else
    {
        Matrix ans = bace ^ n;
        cout << ans.a[0][1] << endl;//f[n] = ans.a[0][0] * f[0] + ans.a[0][1] * f[1],f[0] = 0,f[1] = 1
    }
    return 0;
}

2、求Sn

我们先来看看式子的推导过程：

在这里插入图片描述

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 2, MOD = 1e9 + 7;

struct Matrix
{
    ll a[2][2];
    Matrix()
    {
        memset(a, 0, sizeof a);
    }

    void init()
    {
        for (int i = 0; i < 2; i++)
            a[i][i] = 1;
    }

    Matrix operator*(const Matrix& B)const
    {
        Matrix C;
        for (int i = 0; i < 2; i++)
            for (int k = 0; k < 2; k++)
                for (int j = 0; j < 2; j++)
                    C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
        return C;
    }

    Matrix operator^(const int& t)const
    {
        Matrix A = (*this), res;
        res.init();
        int p = t;
        while (p)
        {
            if (p & 1)res = res * A;
            A = A * A;
            p >>= 1;
        }
        return res;
    }
};

Matrix bace;

int main()
{
    int a, b, n, m;
    cin >> a >> b >> n >> m;
    bace.a[0][0] = a, bace.a[0][1] = b;
    bace.a[1][0] = 1, bace.a[1][1] = a;
    Matrix t = bace ^ n;
    int Xn = t.a[0][0];
    cout << (Xn * 2) % MOD << endl;
    return 0;
}

3、梦想之弧

先看看推导过程：

在这里插入图片描述

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int MOD = 1e9 + 7;

int A0, AX, AY, B0, BX, BY;

struct Matrix
{
	int a[6][6];

	Matrix()
	{
		memset(a, 0, sizeof a);
	}

	void init()
	{
		for (int i = 0; i < 5; i++)
			a[i][i] = 1;
	}

	Matrix operator*(const Matrix& B)const
	{
		Matrix C;
		for (int i = 0; i < 5; i++)
			for (int k = 0; k < 5; k++)
				for (int j = 0; j < 5; j++)
					C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
		return C;
	}

	Matrix operator^(const int& t)const
	{
		Matrix A = (*this), res;
		res.init();
		int p = t;
		while (p)
		{
			if (p & 1)res = res * A;
			A = A * A;
			p >>= 1;
		}
		return res;
	}
};

Matrix bace, once;

int main()
{
	int n;
	cin >> n;
	cin >> A0 >> AX >> AY;
	cin >> B0 >> BX >> BY;

	bace.a[0][0] = AX * BX, bace.a[0][4] = 1;
	bace.a[1][0] = AX * BY, bace.a[1][1] = AX;
	bace.a[2][0] = AY * BX, bace.a[2][2] = BX;
	bace.a[3][0] = AY * BY, bace.a[3][1] = AY, bace.a[3][2] = BY, bace.a[3][3] = 1;
	bace.a[4][4] = 1;

	int a1 = A0 * AX + AY, b1 = B0 * BX + BY;
	once.a[0][0] = a1 * b1, once.a[0][1] = a1, once.a[0][2] = b1, once.a[0][3] = 1, once.a[0][4] = A0 * B0;

	Matrix ans = once * (bace ^ (n - 1));
	cout << ans.a[0][4] << endl;
	return 0;
}

三、矩阵乘法的优化

在上边我们用到的矩阵乘法是三重循环的，时间复杂度为O(n3)，在进行较大的矩阵运算时便会显得吃力。下边我们介绍两种方法来对矩阵的乘法进行优化。

1、省去为0的运算

Matrix operator*(const Matrix& B)const
{
    Matrix C;
    for(int i = 0;i < n;i++)
    {
        for(int k = 0;k < n;k++)
        {
            if(!a[i][k])continue;//为0时乘积为0，直接跳过
            for(int j = 0;j < n;j++)
            {
                C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
            }
        }
    }
}

2、利用二进制优化

这个方法需要用到模板类bitset，我们先来学习一下待会要用到的语法。

  • bitset是将数据用二进制储存下来并且进行编辑。在定义的时候，我们需要规定其长度，例如bitset<10>a，就是创建一个十位的二进制空间a。
bitset<10>a;//创建，a的初始状态为：0000000000

a.set(3);//将第3位变成1，此时a的状态为：0000000100

a.reset();//将所有位清零，此时a的状态为：0000000000

a.set(4),a.set(5),a.set(2);

a.count();//统计a中的1的数量，此时a的状态为：0000011010，a.count() = 3

在掌握以上语法之后，我们就可以来看看它的优化方法了：

首先我们创建两个数组：

bitset<1000>bt[805][3],ct[805][3];

在这里插入图片描述

然后，我们来看看它的运算过程：

在这里插入图片描述

如此如此便可得到结果，我们来看看最终的代码，这时候的时间复杂度为O(n2)：

#include<iostream>
#include<bitset>
using namespace std;
const int N = 805;
bitset<1000>bt[N][3], ct[N][3];

int main()
{
	int n, c;
	cin >> n;
	for (int i = 1; i <= n; i++)
		for (int j = 0; j < 3; j++)
			bt[i][j].reset(), ct[i][j].reset();

	for(int i = 1;i <= n;i++)
		for (int j = 1; j <= n; j++)
		{
			cin >> c;
			bt[i][c % 3].set(j);
		}

	for(int i = 1;i <= n;i++)
		for (int j = 1; j <= n; i++)
		{
			cin >> c;
			ct[j][c % 3].set(i);
		}

	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			int c = ((bt[i][1] & ct[j][1]).count()
				+ 2 * (bt[i][1] & ct[j][2]).count()
				+ 2 * (bt[i][2] & ct[j][1]).count()
				+ 2 * (bt[i][2] & ct[j][2]).count()) % 3;
			cout << c << ' ';
		}
		cout << endl;
	}
	return 0;

}

四、矩阵快速幂加速DP

1、Clarke and digits

这道题可是重量级。它运用的是矩阵快速幂和数位DP（数位DP将在之后吉吉会详细地讲解）。

首先 ，我们要考虑的是如何表示状态，根据题干，我们可以提取三个信息：长度、模数、数位。

我们创建一个三位数组，每一维度分别表示一条信息。如==a[i] [j] [k]==表示的长度为i、模数为j且尾数为k的数的数量。

下面来看看它的表示原理和状态更新的方式：

在这里插入图片描述

但是矩阵应该由二维数组表示，那我们开看看，能不能通过某种方式来减去一个维度。答案当然是可以的，先来看看压缩之后的结果：

a[i] [j] [k] = a[i] [j * 10 + k]

原理何在？我们知道，j是对7取模之后的结果，所以它的范围为[0,7]；k是尾数，尾数即最后一位数，所以它的范围为[0,9]。所以数组a的两维可以表示为一个7 * 10的矩阵，而上图的压缩方法是将其变成了一个1 * 70的矩阵，0 ~ 9表示原来的第一行的状态，10 ~ 19表示原来的第二行的状态，以此类推。那么，我们就得到一个二维数组，即我们需要的矩阵了。

学到这里，结果就明了了，构建完状态矩阵和初始矩阵之后，其他的套模板即可：

状态矩阵：

void init(int k)
{
    memset(bace.a,0,sizeof bace.a);
    for(int i = 0;i < 70;i++)
    {
        int x = i / 10,y = i % 10;//余数和尾数
        for(int j = 0;j < 70;j++)
        {
            int a = j / 10,b = j % 10;//上一状态的余数和尾数
            if((a * 10 + y) == x && y + b != k)
                bace.a[j][i] = 1;
        }
    }
}

初始矩阵：

for(int i = 1;i < 10;i++)once.a[0][(i % 7) * 10 + i] = 1;

完整代码如下：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 71, MOD = 1e9 + 7;

struct Matrix
{
    int a[N][N];
    Matrix()
    {
        memset(a, 0, sizeof a);
    }

    void init()
    {
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                a[i][j] = (i == j);
    }

    Matrix operator*(const Matrix& B)const
    {
        Matrix C;
        for (int i = 0; i < N; i++)
            for (int k = 0; k < N; k++)
                for (int j = 0; j < N; j++)
                    C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
        return C;
    }

    Matrix operator^(const int& t)const
    {
        Matrix A = (*this), res;
        res.init();
        int p = t;
        while (p)
        {
            if (p & 1)res = res * A;
            A = A * A;
            p >>= 1;
        }
        return res;
    }
};

Matrix bace, once;

void init(int k)
{
    memset(bace.a, 0, sizeof bace.a);
    for (int i = 0; i < 70; i++)
    {
        int x = i / 10, y = i % 10;
        for (int j = 0; j < 70; j++)
        {
            int a = j / 10, b = j % 10;
            if ((a * 10 + y) % 7 == x && y + b != k)
                bace.a[j][i] = 1;
        }
    }

    for (int i = 0; i < 10; i++)bace.a[i][70] = 1;
    bace.a[70][70] = 1;
}

int main()
{
    memset(once.a, 0, sizeof once.a);
    for (int i = 1; i < 10; i++)once.a[0][(i % 7) * 10 + i] = 1;
    int t;
    cin >> t;
    while (t--)
    {
        int l, r, k;
        cin >> l >> r >> k;
        init(k);
        Matrix ans1 = once * (bace ^ (l - 1));
        Matrix ans2 = once * (bace ^ r);
        ll ans = ans2.a[0][70] - ans1.a[0][70];
        ans = (ans % MOD + MOD) % MOD;
        cout << ans << endl;
    }
    return 0;
}

2、我也不知道这题叫啥，因为找不到

这里没有什么变编程上的语法知识，我们直接看看推导过程即可：

在这里插入图片描述

以下是费马小定理的衍生：

在这里插入图片描述

我们在求Kn时便可对p-1取模，防止数据过大

代码如下：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 3;
int n, a, b, c, p;

struct Matrix
{
	int a[N][N];

	Matrix()
	{
		memset(a, 0, sizeof a);
	}

	void init()
	{
		for (int i = 0; i < N; i++)
			a[i][i] = 1;
	}

	Matrix operator*(const Matrix& B)const
	{
		Matrix C;
		for (int i = 0; i < N; i++)
			for (int k = 0; k < N; k++)
				for (int j = 0; j < N; j++)
					C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % p;
		return C;
	}

	Matrix operator^(const int& t)const
	{
		Matrix A = (*this), res;
		int p = t;
		res.init();
		while (p)
		{
			if (p & 1)res = res * A;
			A = A * A;
			p >>= 1;
		}
		return res;
	}
};

ll pow(int a, int b)
{
	int ans = 1;
	while (b)
	{
		if (b & 1)ans = (ans * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return 1;
}

Matrix bace, once;

void init()
{
	memset(bace.a, 0, sizeof bace.a);
	memset(once.a, 0, sizeof once.a);
	bace.a[0][0] = c, bace.a[0][1] = 1, bace.a[0][2] = 1;
	bace.a[1][0] = 1, bace.a[2][2] = 1;
	once.a[0][0] = b, once.a[1][0] = 0, once.a[2][0] = b;
}

int main()
{
	cin >> n >> a >> b >> c >> p;
	if (n == 1)puts("1");
	else if (n == 2)cout << pow(a, b) << endl;
	else
	{
		p--;
		Matrix ans = bace ^ (n - 2);
		p++;
		ans = once * ans;
		cout << ans.a[0][0] << endl;
	}
	return 0;
}

在这里插入图片描述

Kicamon
关注
  • 2
    点赞
  • 踩
  • 2
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
专栏目录
【数论】矩阵快速幂，递推优化，模板
IH_LZH的博客
08-12 259
矩阵快速幂推导及模板
DP 优化 - 矩阵快速幂优化
Evan_song1234的博客
12-17 248
以矩阵快速幂方式优化 DP。
C++矩阵优化算法
曹营
11-30 3065
写在前面的话：此矩阵优化算法并非原创，笔者只是转述了一种智慧。  在C++中，定义一个矩阵通常是这样的:    class   MyMatrix    {            ........            public:            ........            float   data[50000];    }             
P4159 [SCOI2009] 迷路（矩阵快速幂，两点路径为k的方案数）
tomjobs的博客
01-19 263
思路： 具体参考https://www.luogu.com.cn/blog/qq-2056188203/mi-lu-scoi2009-ti-xie 简而言之，就是如果权值为1，要求两点之间经过kkk条边的路径方案数，只要将邻接矩阵进行kkk次方就好了。 本题权重为1~9，我们将每个点拆成10个点，两个点边权就通过拆成的点建边来表示，这样就成了权值为1的邻接矩阵形式了。 #include <cstdio> #include <cstring> #include <algorit.
常系数线性递推——多项式取模优化矩阵快速幂
deadpool66的博客
08-03 1040
题目：bzoj_4161 题目描述： &amp;amp;amp;nbsp;&amp;amp;amp;nbsp;&amp;amp;amp;nbsp;&amp;amp;amp;nbsp;&amp;amp;amp;nbsp;给定数列 {hn}前k项，其后每一项满足hn = a1*h(n-1) + a2*h(n-2) + … + ak*h(n-k)其中 a1,a2…ak 为给定数列。请计算 h(n)，并将结果对 1000000007 取模输出。 题解：
华东交通大学大学物理历年试卷以及上课ppt.rar
11-26
这是华东交通大学历年大学物理相关资源，包括试卷及相关答案，以及上课老师讲课相关的ppt,以供大学复习的时候使用，这些都是精选自各个学长之手，希望对大家有帮助。
华东交通大学829数据结构.pdf
08-13
华东交通大学829数据结构历年真题
华东交通大学-课程设计
05-21
2018年华东交通大学编译原理课程设计，里面有自己的课程设计报告。还要刘力月老师的课程设计题目，要求，报告模板。
华东交通大学2013-2019年数据结构.pdf
10-25
华东交通大学2013-2019年数据结构.pdf
华东交通大学2020年数据结构.pdf
10-25
华东交通大学2020年数据结构.pdf
矩阵快速幂优化递推式
宋宇超的博客
01-16 200
参考 练习
洛谷 P1939 【模板】矩阵加速（数列）：优化递推式的方法——矩阵快速幂
Karshilov的博客
12-15 1692
在大多数情况下，O(n)的效率都是值得骄傲的，然而，有时候并不是，比如如何在一秒钟内算出一个递推式的第1e9项，很明显O(n)不行了。 然而常数级又不太现实，除非你的数学非常好，这题又比较简单，你推了一个特征方程的通项公式…… 所以考虑log的做法：矩阵快速幂 如果你还不知道矩阵快速幂是什么，请走这边：传送门 对于这道题，嗯，模板嘛，已经告诉你了式子，就只需要考虑矩阵了，对于整个过程，我们
HDU 5318 (dp+矩阵快速幂优化)
努力
07-29 2204
The Goddess Of The Moon Time Limit: 6000/3000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others) Total Submission(s): 605    Accepted Submission(s): 259 Problem Description Chan
矩阵快速幂优化的动态规划
Leo_h
04-12 2078
为了防止写错
行列式运算法则 矩阵的运算及其运算规则:
热门推荐
ZJQ的博客
10-29 8万+
1、三角形行列式的值，等于对角线元素的乘积。计算时，一般需要多次运算来把行列式转换为上三角型或下三角型 2、交换行列式中的两行（列），行列式变号（交换） 3、行列式中某行（列）的公因子，可以提出放到行列式之外。（倍乘）（注：矩阵是全部元素都乘，都提取） 4、行列式的某行乘以a，加到另外一行，行列式不变，常用于消去某些元素。（倍加） 5、若行列式中，两行（列）完全一样，则行列式为0；可以推论，如果两...
矩阵乘法的分块优化
z-k的博客
03-04 1万+
普通矩阵乘法一般采用3重循环完成。 void dgemm (int n, double* A, double* B, double* C) { for (int i = 0; i &lt; n; ++i) for (int j = 0; j &lt; n; ++j) { double cij = C[i+j*n]; /* cij = C[i][j] *...
矩阵的特性和运算法则
weixin_42394591的博客
06-19 7363
基本运算法则：http://www2.edu-edu.com.cn/lesson_crs78/self/j_0022/soft/ch0605.html行列式的值的计算方法：点击进入逆矩阵：点击进入伴随矩阵：点击进入矩阵求导：点击进入...
hdu3802 矩阵快速幂&特征方程&降幂
bigtiao097的博客
08-12 594
矩阵快速幂，特征方程，降幂公式
指数型生成函数的范例——状压dp、矩阵快速幂优化及对指数型生成函数的理解
hans774882968的博客
08-06 496
从《组合数学》中的经典例题开始。 在n个排成一排的瓷砖上染色，有4种颜色，有2种要求出现次数为2的倍数，其他的无要求。求方案数，模998244353。 组合计数问题先看能否dp。显然可以考虑“前缀定义dp”，于是我们还要考虑当前这个颜色对每种有要求的颜色的奇偶性的影响。那么定义dp[i][S]为涂了前i个瓷砖，两种颜色出现次数各自模2的状态位为S的答案。则答案为dp[n][0]。转移分2种情况： 当前涂了无要求的颜色，贡献为dp[i-1][S] 当前涂了有要求的颜色，设是第j种有要求的颜色，j=0~1
华东交通大学数据结构829bbs
最新发布
10-07
论坛的名字"829bbs"代表的是华东交通大学数据结构课程的课程编号，显示了论坛的主题和定位。在论坛的主页上，用户可以看到各种不同的板块，如数据结构基础、算法分析与设计、实践项目等。每个板块上都有相应的帖子，...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值