【算法设计-分治】快速幂与龟速乘、矩阵乘与矩阵快速幂
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}\)
