算法设计-分治快速幂与龟速乘

Posted 2023-03-07 Mount256

文章目录

• • 1. 快速幂
  • 2. 龟速乘
  • 3. 快速幂取模
  • 4. 龟速乘取模
  • 5. 快速幂取模优化

1. 快速幂

算法原理：

• 计算 3¹¹：
  • 3¹¹ = (3⁵)² x 3
  • 3⁵ = (3²)² x 3
  • 3² = 3 x 3
  • 仅需计算 3 次，而非 11 次
• 计算 3¹⁰：
  • 3¹⁰ = (3⁵)²
  • 3⁵ = (3²)² x 3
  • 3² = 3 x 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²=9
• 第 2 轮：（5 % 2 == 1）r=3x3²=3³=27，b=2，a=(3²)²=3⁴=81
• 第 3 轮：（2 % 2 == 0）r 不变，b=1，a=(3⁴)²=3⁸
• 第 4 轮：（1 % 2 == 1）r=3³x3⁸=3¹¹，b=0，a=(3⁸)²=3¹⁶
• 得到 r = 3³x3⁸ = 3¹¹

2. 龟速乘

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

12 x a = 23 x a + 21 x 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

(ab) % 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;