模重複平方演算法

在RSA演算法裡頭經常要用到“求x的n次方模m”這樣的過程，通常使用O(log(n))的模重複平方演算法來實現，提高效率。

其實這是大二學的《資訊安全數學基礎》裡面的內容，那時為了考試需要（手算+寫出很羅嗦的過程），還專門寫了代碼放在Blog空間裡考試的時候用—……

同樣O(log(n))的遞迴演算法其實很容易理解：

/* C */int square(int x) { return x * x; } /* 這裡可不能用宏喲 */int fast_mod(int x, int n, int m){    if (n == 0)        return 1;    else if (n % 2 == 0)        return square(fast_mod(x, n / 2, m)) % m;    else         return x * fast_mod(x, n - 1, m) % m;}#Pythonfast_mod = lambda x, n, m: 1 if n == 0 else fast_mod(x, n / 2, m) ** 2 % m if n % 2 == 0 else x * fast_mod(x, n - 1, m) % m;Scheme(define (even? n) (= (remainder n 2) 0))(define (mod-fast x n m)    (define (square x) (* x x))     (cond ((= n 0) 1)          ((even? n) (remainder (square (mod-fast x (/ n 2) m)) m))           (else (remainder (* x (mod-fast x (- n 1) m)) m))))(mod-fast 79 24 33);

但是SICP要求寫一個迭代版的。印象中我是記得可以把 n 寫成二進位（比如n=13， 1101），然後一位一位地推。

試推了一下從高位到低位，倒是挺簡單的，記B[i]為第 i 位的值，通過A[i] = A[i+1]^2 * x^B[i] 從高到低計算出A[0]，就得到結果了。但是問題是，為了從高到低計算，又得一次遞迴，似乎不能滿足要求。

於是只好反過來，從低位往高位推。這個過程其實也挺簡單的，舉例來說：

當n=13，即二進位的 1101 = 1 * 2^3 + 1 * 2^2 + 0 * 2^1 + 1 * 2^0時，最終結果 

ans = x^n % m
    = x^(1 * 2^3 + 1 * 2^2 + 0 * 2^1 + 1 * 2^0) % m
    = [x^(1 * 2^3)] * [x^(1 * 2^2)] * [(x ^ (0 * 2^1)] * [x ^ (1 * 2^0)] % m

也就是說，從低到高，在第 i 位的時候，將 x^(Bit[i] * 2^i) % m 乘到結果中即可。這裡可以稍微變換一下：僅當Bit[i] == 1的時候，將x^(2^i) % m乘進去即可。所以這裡可以用一個輔助的變數 b 來儲存 x^(2^i) % m，在每次迭代的過程中 b =
b^2 % m 。

於是實現就容易了：

/* C */int fast_mod_iter(int x, int n, int m){    int a = 1, b = x; //i=0的時候b = x^(2^0) = x    while (n)    {        if (n % 2 == 1)            a = a * b % m;        b = b * b % m;        n /= 2;    }    return a;}#Pythondef fast_mod(x, n, m):     a = 1     b = x     while True:        if n == 0:            return a        if n % 2 == 1:            a = a * b % m         b = b * b % m         n /= 2;Scheme(define (even? n) (= (remainder n 2) 0))(define (mod-fast-iter x n m)    (define (iter a b n)        (cond ((= n 0) a)              ((even? n)                (iter a (remainder (* b b) m) (/ n 2)))              (else                (iter (remainder (* a b) m) (* b b) (/ (- n 1) 2)))))    (iter 1 x n))(mod-fast-iter 79 24 33);

轉自大牛 http://www.felix021.com/blog/read.php?2112