【算法设计-分治】快速幂与龟速乘、矩阵乘与矩阵快速幂

目录

• 1. 快速幂
• 2. 龟速乘
• 3. 快速幂取模
• 4. 龟速乘取模
• 5. 快速幂取模优化
• 6. 矩阵乘
• 7. 矩阵快速幂
• 8. 矩阵快速幂用于递推式
  • 斐波那契数列
  • 递推式的一些套路
  • 参考资料

1. 快速幂

算法原理：

• 计算$3^{11}$：
  • $3^{11} = (3^5)^2 \times 3$
  • $3^5 = (3^2)^2 \times 3$
  • $3^2 = 3 \times 3$
  • 仅需计算 3 次，而非 11 次
• 计算 $3^{10}$：
  • $3^{10} = (3^5)^2$
  • $3^5 = (3^2)^2 \times 3$
  • $3^2 = 3 \times 3$
  • 仅需计算 3 次，而非 10 次

算法思路：

• 若指数是偶数，则将底数平方，指数除以 2。
• 若指数是奇数，则将底数平方，指数除以 2，再乘上底数。

算法代码：

typedef unsigned long long uLL;

// 快速幂 a^b
uLL power (uLL a, uLL b){
	uLL r = 1;
	while (b != 0){
		if (b & 1) 	// (b % 2 == 1)
			r = r * a;
		b = b >> 1; // (b = b / 2)
		a = a * a;
	}
	return r;
}

举例：

• 初始值：a = 3，b = 11
• 第 1 轮：（11 % 2 == 1）r=1x3=3，b=5，a=$3^2$=9
• 第 2 轮：（5 % 2 == 1）r=3x$3^2$=$3^3$=27，b=2，a=$(3^2)^2$=$3^4$=81
• 第 3 轮：（2 % 2 == 0）r 不变，b=1，a=$(3^4)^2$=$3^8$
• 第 4 轮：（1 % 2 == 1）r=$3^3$x$3^8$=$3^{11}$，b=0，a=$(3^8)^2$=$3^{16}$
• 得到 r = $3^3 \times 3^8$ = $3^{11}$

2. 龟速乘

算法原理：将其中一个乘数分解成 2 的幂次相加。

$12 \times a = 2^3 \times a + 2^1 \times a$

算法代码：

typedef unsigned long long uLL;

// 龟速乘 a*b
uLL mul (uLL a, uLL b){
	uLL r = 0;
	while (b != 0){
		if (b & 1) 	// (b % 2 == 1)
			r = r + a;
		b = b >> 1; // (b = b / 2)
		a = a + a;
	}
	return r;
}

3. 快速幂取模

初等数论中有如下公式：

(a × b) % m = ((a % m) × (b % m)) % m

推广：

(a × b × c...) % m = ((a % m) × (b % m) × (c % m) × ... ) % m

(a^b) % m = (a × a × a...) % m = ((a % m) × (a % m) × (a % m) × ... ) % m

算法代码：

typedef unsigned long long uLL;

// 快速幂取模 (a^b) % p
uLL powerMod (uLL a, uLL b, uLL p){
	uLL r = 1;
	while (b != 0){
		if (b & 1) 	// (b % 2 == 1)
			r = (r * a) % p;
		b = b >> 1; // (b = b / 2)
		a = (a * a) % p;
	}
	return r;
}

4. 龟速乘取模

算法原理：(a × b) % m = ((a % m) × (b % m)) % m

算法代码：

// 龟速乘取模 (a*b) % p
uLL mulMod (uLL a, uLL b, uLL p){
	uLL r = 0;
	while (b != 0){
		if (b & 1) 	// (b % 2 == 1)
			r = (r + a) % p;
		b = b >> 1; // (b = b / 2)
		a = (a + a) % p;
	}
	return r;
}

5. 快速幂取模优化

算法原理：注意到快速幂取模算法中的相乘操作可能会超出数据范围，因此可以将相乘操作转化为龟速乘取模。

原理依然是此公式：(a × b) % m = ((a % m) × (b % m)) % m，其中((a % m) × (b % m))即为龟速乘取模。

算法思路：快速幂 + 龟速乘结合。

// 快速幂取模防止爆炸 (a^b) % p
uLL powerModBig (uLL a, uLL b, uLL p){
	uLL r = 1;
	while (b != 0){
		if (b & 1) 	// (b % 2 == 1)
			r = mulMod(a, b, p) % p;
		b = b >> 1; // (b = b / 2)
		a = mulMod(a, a, p) % p;
	}
	return r;
}

6. 矩阵乘

#include <cstdio>
#include <iostream>
using namespace std;

#define MAX 100

struct Martix{
	int row;	// 行 
	int col;	// 列 
	int martix[MAX][MAX];
	Martix (int r, int c): row(r), col(c) {} // 构造函数，此时不能使用typedef！ 
};

Martix Multiply (Martix x, Martix y){
	Martix z(x.row, y.col);
	
	for (int i = 0; i < z.row; i++){
		for (int j = 0; j < z.col; j++){
			z.martix[i][j] = 0;
			for (int k = 0; k < x.col; k++)
				z.martix[i][j] += x.martix[i][k] * y.martix[k][j];
		}
	}
			
	return z;
}

int main(){
	int r, c;
	
	printf("请输入第一个矩阵的行和列：\n");
	scanf("%d %d", &r, &c);
	Martix x(r, c);
	printf("请输入%d行%d列的矩阵：\n", r, c);
	for (int i = 0; i < x.row; i++)
		for (int j = 0; j < x.col; j++)
			scanf("%d", &x.martix[i][j]);
	
	printf("请输入第二个矩阵的行和列：\n");
	scanf("%d %d", &r, &c);
	Martix y(r, c);
	printf("请输入%d行%d列的矩阵：\n", r, c);
	for (int i = 0; i < y.row; i++)
		for (int j = 0; j < y.col; j++)
			scanf("%d", &y.martix[i][j]);
	
	Martix result = Multiply(x, y);
	
	printf("结果是：\n");
	for (int i = 0; i < result.row; i++){
		for (int j = 0; j < result.col; j++)
			printf("%d ", result.martix[i][j]);
		printf("\n");
	}
	
	return 0;
}

7. 矩阵快速幂

#include <cstdio>
#include <iostream>
using namespace std;

#define MAX 100

struct Martix{
	int row;	// 行 
	int col;	// 列 
	int martix[MAX][MAX];
	Martix (int r, int c): row(r), col(c) {} // 构造函数 
};

Martix Multiply (Martix x, Martix y){
	Martix z(x.row, y.col);
	
	for (int i = 0; i < z.row; i++){
		for (int j = 0; j < z.col; j++){
			z.martix[i][j] = 0;
			for (int k = 0; k < x.col; k++)
				z.martix[i][j] += x.martix[i][k] * y.martix[k][j];
		}
	}
			
	return z;
}

Martix Power (Martix x, int k){
	Martix r(x.row, x.col);
	
	// 单位矩阵构建 
	for (int i = 0; i < r.row; i++){
		for (int j = 0; j < r.col; j++){
			if (i == j)
				r.martix[i][j] = 1;
			else
				r.martix[i][j] = 0;
		}
	}
	
	// 矩阵快速幂 
	while (k != 0){
		if (k & 1)
			r = Multiply(x, r);
		k = k >> 1;
		x = Multiply(x, x);
	}
	
	return r;
}

int main(){
	int r, k;
	
	printf("请输入行（列）：\n");
	scanf("%d", &r);
	Martix x(r, r);
	printf("请输入%d行%d列的矩阵：\n", r, r);
	for (int i = 0; i < x.row; i++)
		for (int j = 0; j < x.col; j++)
			scanf("%d", &x.martix[i][j]);
			
	printf("请输入指数k：\n");
	scanf("%d", &k);
	
	Martix result = Power(x, k);
	
	printf("结果是：\n");
	for (int i = 0; i < result.row; i++){
		for (int j = 0; j < result.col; j++)
			printf("%d ", result.martix[i][j]);
		printf("\n");
	}
	
	return 0;
}

8. 矩阵快速幂用于递推式

斐波那契数列

数列递推式：$f_{n} = f_{n-1} + f_{n-2}$

矩阵递推式：$\begin{bmatrix} f_n \\ f_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} f_{n-1} \\ f_{n-2} \end{bmatrix}$

矩阵公式：$\begin{bmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^n$

#include <cstdio>
#include <iostream>
using namespace std;

#define MAX 100

struct Martix{
	int row;	// 行 
	int col;	// 列 
	int martix[MAX][MAX];
	Martix (int r, int c): row(r), col(c) {} // 构造函数 
};

Martix Multiply (Martix x, Martix y){
	Martix z(x.row, y.col);
	
	for (int i = 0; i < z.row; i++){
		for (int j = 0; j < z.col; j++){
			z.martix[i][j] = 0;
			for (int k = 0; k < x.col; k++)
				z.martix[i][j] += x.martix[i][k] * y.martix[k][j];
		}
	}
			
	return z;
}

Martix Power (Martix x, int k){
	Martix r(x.row, x.col);
	
	// 单位矩阵构建 
	for (int i = 0; i < r.row; i++){
		for (int j = 0; j < r.col; j++){
			if (i == j)
				r.martix[i][j] = 1;
			else
				r.martix[i][j] = 0;
		}
	}
	
	// 矩阵快速幂 
	while (k != 0){
		if (k & 1)
			r = Multiply(x, r);
		k = k >> 1;
		x = Multiply(x, x);
	}
	
	return r;
}

int main(){
	int k;
		
	while (scanf("%d", &k) != EOF){	
		// 斐波那契数列的递推矩阵构建 
		// [1, 1]
		// [1, 0]
		Martix e(2, 2);
		e.martix[0][0] = 1; e.martix[0][1] = 1;
		e.martix[1][0] = 1; e.martix[1][1] = 0;
		
		Martix result = Power(e, k);
		
		for (int i = 0; i < result.row; i++){
			for (int j = 0; j < result.col; j++)
				printf("%d ", result.martix[i][j]);
			printf("\n");
		}
	}
		
	return 0;
}

输入和输出：

1
1 1
1 0
2
2 1
1 1
3
3 2
2 1
4
5 3
3 2
5
8 5
5 3
6
13 8
8 5
7
21 13
13 8
8
34 21
21 13
9
55 34
34 21
10
89 55
55 34

相关题目：（POJ3070）斐波那契数列f(n)，输入n，求f(n) % 10000，n <= 1e9。

解题思路：仅需在快速幂函数中添加% 10000即可。

递推式的一些套路

（1）数列递推式：$f_{n} = a \times f_{n-1} + b \times f_{n-2} + c$

矩阵递推式：$\begin{bmatrix} f_n \\ f_{n-1} \\ c \end{bmatrix} = \begin{bmatrix} a & b & 1 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} f_{n-1} \\ f_{n-2} \\ c \end{bmatrix}$

（2）数列递推式：$f_{n} = a \times f_{n-1} + b \times f_{n-2} + c^n$

矩阵递推式：$\begin{bmatrix} f_n \\ f_{n-1} \\ c^n \end{bmatrix} = \begin{bmatrix} a & b & c\\ 1 & 0 & 0 \\ 0 & 0 & c \end{bmatrix} \times \begin{bmatrix} f_{n-1} \\ f_{n-2} \\ c^{n-1} \end{bmatrix}$

（3）数列递推式：$f_{n} = c^n - f_{n-1}$

矩阵递推式：$\begin{bmatrix} f_n \\ c^n \end{bmatrix} = \begin{bmatrix} -1 & c \\ 0 & c \end{bmatrix} \times \begin{bmatrix} f_{n-1} \\ c^{n-1} \end{bmatrix}$

参考资料

基础算法—快速幂详解

矩阵快速幂（原理+模板）

KY21 递推数列