模取幂（求一个数的幂对另一个数的模运算）

即求a^b mod n的值，a/b为非负整数，n是一个正整数

设<bk, bk-1, ..., b1, b0>是b的二进制表示，随着c的值从0到b成倍增长，下面过程最终计算出a ^ b mod n

MODULAR-EXPONENTIATION(a, b, n)

    c = 0

    d = 1                    //初始化

    let <bk, bk-1, ..., b1, ...b0> be the binary representation of b

    for i = k downto 0

        c = 2c

        d = (d * d) mod n

        if bi == 1

            c = c + 1

            d = (d * a) mod n

    return d

下面证明算法的两个性质【循环不变式】

1.    c的值与b的二进制表示的前缀<bk, bk-1, ..., bi+1>相同

2.    d = a^c mod n

初始化：最初，i = k，前缀<bk, bk-1, ..., bi+1>为空，对应于c=0，且d = 1 = a^0 mod n

保持：

    令c'和d'为for循环的一次迭代的结束处c和d的值，则它们是下一次迭代前的值；

    每次更新迭代，c' = 2c(bi=0)或c' = 2c + 1(bi=1)，使得c满足条件1

    【以下过程建立在基础上：[(a mod n) * (b mod n)] mod n = (a * b) mod n】

    如果bi = 0，d' = d^2 mod n = (a^c mod n)^2 mod n =  (a^c)^2 mod n = a^2c mod n = a^c' mod n

    如果bi = 1，d' = (d^2 mod n) * a mod n = (d^2 * a) mod n = (a^c)^2 * a mod n = a^(2c+1) mod n = a^c' mod n

    从而条件二成立

终止：

    i = -1，c = b；从而d = a^c mod n = a^b mod n；即d为结果

时间复杂度：a,b,n 都是 k位的数，需要算数运算的总次数是O(k)，而需要的位操作的总次数是O(k^3)，因为每一轮迭代中有d*d