高速な累乗計算

累乗(x^n)を単純に計算すると、オーダーは O(n)となり効率が悪いです。そこで、nを2の累乗に分解して計算する高速化手法が一般に知られています。

たとえば、3 の 11 乗を計算する場合を考えましょう。11 は 1 + 2 + 8 に分解できます。

この累乗の系列では、ある値は一つ前の値を2乗することで計算できます。たとえば、こうです。

3^1 = 3
3^2 = 3 * 3 = 9
3^4 = (3^2)^2 = 9^2 = 81
3^8 = (3^4)^2 = 81^2 = 6561

よって、3^11 は以下のように計算できます。

3^11 = 3^(1+2+8)
     = 3^1 × 3^2 × 3^8
     = 3 × 9 × 6561
     = 177147

この方法のオーダーは、O(log2(n)) です。

遅延評価風に

RSA のために高速な累乗計算を Haskell で実装したことがありました。そのときのコードは、こんな感じです。

import Data.Bits
power x n = iter n seq 1
    where
      seq = iterate (\y -> y*y) x
      iter 0 _ ret = ret
      iter b (s:ss) ret = let next = if odd b
                                     then ret * s
                                     else ret
                          in iter (shiftR b 1) ss next

遅延評価を生かして、x を2乗、2乗と乗じていく無限数列(seq)を作っているところがポイントです。

このコードはちゃんと動きますが、右シフト(shiftR)を使っているので、ダサいなぁと思っていました。(いや、`mod` 2 を使ってもいいんですが、それでもダサいですね。)

美しいアルゴリズム

現在輪講で読んでいる SOE に、累乗を高速に計算でき、しかも美しいアルゴリズムが載っていました。それは、こうです。

power x n
    | n == 0    = 1
    | even n    = power (x*x) (n `div` 2)
    | otherwise = x * power x (n - 1)

分かりますか？

• 偶数のときは、結果は変わらず、x を 2 乗し、n を半分にする
  • x^6 = (x^2)^3
• 基数のときは、結果に x を掛け、ｘはそのままで、n を一つ減らす
  • x^7 = x × x^6

これが理解できれば、アルゴリズムが正しいのは明らかです。

オーダーは O(log2(n)) です。なぜなら、n が奇数のときは、1を引くので、次は必ず偶数に。n が偶数のときは、半分にして、次は偶数か奇数に。という訳で、最低でも2回に一回は、n のビット数が減るのです。

美しいなぁ。感動した！

The Haskell School of Expression: Learning Functional Programming through Multimedia