快速幂算法实现

问题：如何快速计算

x^n

，

x^{-n}

？

提示：

x^{-n} = (\frac{1}{x})^n

快速幂，二进制取幂（Binary Exponentiation，也称平方法）是一个在 $O(\log n)$ 时间内计算 $x^n$ 的技巧性算法。快速幂算法的本质是分治算法。

递归版本

这个版本能直观理解计算过程。

// 计算x的n次方，x可能为负数
Pow(x, n):
	return n >= 0 ? QUICK-MUL(n,x) : 1.0 / QUICK-MUL(-n,x)

QUICK-MUL(N,x)
	if N == 0:
		return 1.0
	y = QUICK-MUL(N / 2 ,x) // 这里是整除，存在向下取整的过程
	return N % 2 == 0 ? y * y : y * y * x

复杂度分析：

时间复杂度： $O(\log n)$ ，即为递归的层数。
空间复杂度： $O(\log n)$ ，即为递归的层数。这是由于递归的函数调用会使用栈空间。

迭代版本

假设要计算 $x^9$ ：

形式化来讲：

\begin{aligned} n&=2^{i_{0}}+2^{i_{1}}+\dots+2^{i_{k}} \\ x^n&=x^{2^{i_{0}}}\times x^{2^{i_{1}}}\times\dots \times x^{2^{i_{k}}} \end{aligned}

因此算法可以这样写，把 n 看成一个二进制数，从 n 的低位开始依次检查二进制位，如果为 1，则把当前 x 乘到结果中。每次检查下一位时，x 自身翻一番。

def quickMul(N):
	ans = 1.0
	# 贡献的初始值为 x
	x_contribute = x
	# 在对 N 进行二进制拆分的同时计算答案
	while N > 0:
		if N % 2 == 1:
			# 如果 N 二进制表示的最低位为 1，那么需要计入贡献
			ans *= x_contribute
		# 将贡献不断地平方
		x_contribute *= x_contribute
		# 舍弃 N 二进制表示的最低位，这样我们每次只要判断最低位即可
		N //= 2
	return ans

复杂度分析：

时间复杂度： $O(\log n)$ ，即为对 n 进行二进制拆分的时间复杂度。
空间复杂度： $O(1)$ 。

快速幂的应用

模指数运算（模意义下的取幂）

推论 TY-6（分配律）：m 是正整数且 a 和 b 是整数，那么有

$(a + b) \bmod m = ((a \bmod m) + (b \bmod m)) \bmod m$
$a\cdot b \bmod m = ((a \bmod m)(b \bmod m)) \bmod m$

这个定理其实也解释了为什么在一些算法题中，当题目要求给出最终计算结果取模时，即便我们在过程中取模也不会出错的原因。

假设我们要计算 $b^n \bmod m$ ，其中 n 是一个二进制表示的数：

\begin{aligned} n&=(a_{k-1}a_{k-1}\dots a_{1}a_{0})_{2} \\ b^n &= b^{a_{k-1}\cdot 2^{k-1}+\dots+a_{1}\cdot 2^{1}+a_{0}} = b^{a_{k-1}\cdot 2^{k-1}}\dots b^{a_{0}} \end{aligned}

需求出 $b,b^2,\dots,b^{2^{k-1}}$ ，再把 $a_{j}=1$ 的那些项相乘。
为了提高效率并减少空间需求，每乘一项后，都做一次模 m 运算以缩小结果值。

MODULAR-EXPONENTIATION(b: int, n: 二进制数 ,m: 正整数)
	x=1
	power = b % m
	for i=0 to k-1
		if n[i]==1
			x = (x*power) % m
		power = (power*power) % m
	return x

此外，我们还可以使用费马小定理加速算法过程。

费马小定理

如果 p 是素数，对每个整数 a 都有 $a^p ≡ a ~(\bmod p)$ ；
如果 p 是素数， $p \not\mid a$ ，那么 $a^{p-1} ≡ 1 ~(\bmod p)$ 。

例子：利用费马小定理化简运算

求 $7^{121}\bmod 13$

解题步骤：

根据费马小定理， $7^{12} ≡ 1 (\mod 13)$
那么根据分配律有 $(7^{12})^k ≡ 1 (\mod 13)$ ，k 是任意整数
$7^{121} = 7^{12 \cdot 10 +1} = (7^{12})^{10} \cdot 7 ≡ 1^{10} \cdot 7 ≡ 7 (\mod 13)$

再捋顺下思路，有： $7^{121} ≡ 7^{121 \bmod 12} (\mod 13)$ 。

即有 $x^n ≡ x^{n\bmod (p-1)} (\bmod p)$ 。

如果指数超过 64 位整数的范围，可尝试欧拉降幂。

矩阵快速幂

🟨 P3390 【模板】矩阵快速幂 - 洛谷

import java.util.*;

public class Main {
    private static int MOD = 1000000007;

    public static void main(String[] args) throws Exception {
        Scanner in = new Scanner(System.in);
        
        // 输入构造
        int n = in.nextInt();
        long k = in.nextLong();
        long[][] A = new long[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                A[i][j] = in.nextInt();
            }
        }

        // 单位矩阵
        long[][] ans = getUnitMatrix(n);

        // 快速幂
        ans = fast_expo(n, k, A, ans);

        // 输出构造
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                System.out.print(ans[i][j] + " ");
            }
            System.out.println();
        }
    }

    private static long[][] fast_expo(int n, long k, long[][] A, long[][] ans) {
        while (k > 0) {
            if ((k & 1) == 1) {
                ans = mult(ans, A, n);
            }
            k >>= 1;
            A = mult(A, A, n);
        }
        return ans;
    }

    private static long[][] getUnitMatrix(int n) {
        long[][] ans = new long[n][n];
        for (int i = 0; i < n; i++) {
            ans[i][i] = 1L;
        }
        return ans;
    }

    public static long[][] mult(long[][] A, long[][] B, int n) {
        long[][] ans = new long[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                for (int k = 0; k < n; k++) {
                    ans[i][j] += (A[i][k] * B[k][j]) % MOD;
                }
                ans[i][j] %= MOD;
            }
        }
        return ans;
    }
}