幂模m与逐次平方法（快速幂算法）

1 计算幂模$m$

  如何计算：
$$5^{100000000000000}(mod12830603)$$
呢？这么大的数字显然是不能硬算的，如果$12830603$是素数，我们可以使用费马小定理简化计算，即便不是素数，也可以使用欧拉公式。但即便简化了计算，剩余的计算量依旧大的离谱：
$$5^{6546640}(mod12830603)$$

  进行这样的大数计算是很实用的，它可以用来做加密，而且产生的密码使用目前最先进的破译技术也不能破译。用来计算$a^k(modm)$的一个巧妙想法叫做逐次平方法。在叙述一般方法之前，先用一个例子体现它的思想。

$$7^{327}(mod853)$$

  第一步创建表格，给出$7,7^2,7^4,7^8,7^{16},\cdots (mod853)$的值。注意，要得到表中的每个项，我们仅仅需要平方前面的数。由于在平方前总是模$853$进行简化，所以从未使数大于$853^2$，下面是$7^{2^k}(mod853)$的列表

$$\begin{array}{rl}& 7^1\equiv 7\equiv 7(mod853)\\ & 7^2\equiv (7^1{)}^2\equiv 7^2\equiv 49\equiv 49(mod853)\\ & 7^4\equiv (7^2{)}^2\equiv 49^2\equiv 2401\equiv 695(mod853)\\ & 7^8\equiv (7^4{)}^2\equiv 695^2\equiv 483025\equiv 227(mod853)\\ & 7^{16}\equiv (7^8{)}^2\equiv 227^2\equiv 51529\equiv 349(mod853)\\ & 7^{32}\equiv (7^{16}{)}^2\equiv 349^2\equiv 121801\equiv 675(mod853)\\ & 7^{64}\equiv (7^{32}{)}^2\equiv 675^2\equiv 455625\equiv 123(mod853)\\ & 7^{128}\equiv (7^{64}{)}^2\equiv 123^2\equiv 15129\equiv 628(mod853)\\ & 7^{256}\equiv (7^{128}{)}^2\equiv 628^2\equiv 394384\equiv 298(mod853)\end{array}$$

  下一步，将指数$327$表成$2$的幂次和：

$$327=256+71=256+64+7=256+64+4+3=256+64+4+2+1$$
现在使用$327$的二进制展开式计算
$$7^{327}=7^{256+64+4+2+1}=298\ast 123\ast 695\ast 49\ast 7(mod853)$$
要完成$7^{327}(mod853)$的计算，我们只需要程这$5$个数，然后用模$853$进行简化。如果所有$5$个数的乘积对我们来说尝试太大，则可以乘前两个数简化，再与第三个数相乘简化，以这种方式，我们就不需要对任何大于$853^2$的数进行计算，因此：
$$\begin{array}{rl}298\ast 123\ast 695\ast 49\ast 7(mod853)& =828\ast 695\ast 49\ast 7(mod853)\\ & =538\ast 49\ast 7(mod853)\\ & =772\ast 7(mod853)\equiv 286(mod853)\end{array}$$

2 逐次平方法

  算法（逐次平方计算$a^k(modm)$）：用下述步骤计算$a^k(modm)$的值：

1. 将$k$表成$2$的幂次和:
   $$k=u_0+u_1\ast 2+u_2\ast 2^2+u_3\ast 2^3+\cdots +u_r\ast 2^r$$
   其中$u_i$为$0$或$1$
2. 使用逐次平方法制作模$m$的$a$的幂次表
   $$\begin{array}{rl}& a^1\equiv A_0(modm)\\ & a^2\equiv (a^1{)}^2\equiv A_0^2\equiv A_1(modm)\\ & a^4\equiv (a^2{)}^2\equiv A_1^2\equiv A_2(modm)\\ & a^8\equiv (a^4{)}^2\equiv A_2^2\equiv A_3(modm)\\ & ⋮\\ & a^{2^r}\equiv (a^{2^{r-1}}{)}^2\equiv A_{r-1}^2\equiv A_r(modm)\end{array}$$
   注意要计算表的每一行，仅需要取前一行最末的数，平方它然后用模$m$简化。也注意到表有$r+1$行，其中$r$是第一步中$k$的二进制展开式中$2$的最高指数。
3. 乘积：
   $$A_0^{u_0}\ast A_1^{u_1}\ast A_2^{u_2}\ast \cdots \ast A_r^{u_r}(modm)$$
   同余于$a^k(modm)$。注意到所有$u_i$是$0,1$。因此这个数实际上是$u_i=1$的那些$A_i$的乘积。

3 算法实现

  下述算法是进行逐次平方法的更有效方法。

算法流程：

• 置$b=1$
• 当$k\ge 1$时，循环：
  1.如果$k$是奇数则$b=a\ast b(modm)$
  2.置$a=a^2(modm)$
  3.置$k=k/2$（如果$k$是奇数，转下）
• 结束循环
• 转到$b$的值

下面来解释一下该算法，首先需要知道一个取模运算的运算法则，即：
$$(a\ast b)\%p=(a\%p\ast b\%p)\%p$$
对于一个大数次幂，可以通过降幂的方法来进行加速计算，降幂方法就是指数除以$2$，底数开平方，例如$a^{11}$，我们可以通过以下步骤来计算：

$$\begin{array}{rl}{a}^{11}& =(a^2{)}^5\ast a\\ & =(a^2{)}^4\ast a^2\ast a\\ & =(a^4{)}^2\ast a^2\ast a\\ & =(a^8{)}^1\ast a^2\ast a\end{array}$$
这个过程我们发现，$a^{11}$的结果就是在降幂过程中，幂为奇数时的各个底数的乘积，即$a^{11},(a^2{)}^5,(a^8{)}^1$的底数$a,a^2,a^8$。

def zcpf(a, k, m):
    b = 1
    a %= m
    while k >= 1:
        if k % 2 == 1:
            b = (a * b) % m
        a = a**2 % m
        k //= 2
    return b

print(zcpf(7, 327, 853))

输出：

286

4 参考资料

《数论概论》第四版 P74-P77