Pythonでロシア農民乗算法を使って複素数のべき乗を高速に計算する方法

この記事では、「ロシア農民乗算法（Russian Peasant Multiplication）」の考え方を応用して、複素数のべき乗を効率よく計算するPythonプログラムを解説します。

問題設定は次のとおりです。4つの整数 p、q、r、k が与えられたとき、複素数 (p + qi) の r 乗を求めます。結果は a + bi の形で表せるので、a mod k と b mod k のペアを答えとして返します。

たとえば p = 3、q = 0、r = 8、k = 10000 が入力された場合を考えてみましょう。q = 0 なので計算は単純な 3 の 8 乗、つまり 6561 になります。したがって出力は (6561, 0) です。

アルゴリズムのポイント

ロシア農民乗算法は、指数を半分ずつ減らしながら計算を進める手法で、いわゆる「繰り返し二乗法（バイナリ法）」と同じ発想です。複素数の場合も、次の性質を利用すれば同様に高速化できます。

(p + qi)² = (p² − q²) + (2pq)i

この式のおかげで、指数 r が偶数なら底を2乗して指数を半分にでき、奇数なら1回だけ掛けて偶数に帰着させられます。その結果、単純に r 回掛け算を繰り返す O(r) の方法に対し、O(log r) で計算が完了します。

解法の手順

• r = 0 のとき： 1 を返します。
• r = 1 のとき： (p mod k, q mod k) を返します。
• r が偶数のとき： 実部を (p*p − q*q) mod k、虚部を 2*p*q mod k として、solve(p, q, r/2, k) を再帰的に呼び出します。
• r が奇数のとき： まず solve(p, q, r−1, k) の結果を (pr, qr) として受け取り、((p*pr − q*qr) mod k, (p*qr + q*pr) mod k) を返します。これは複素数の掛け算 (p + qi)(pr + qr・i) を展開した形に相当します。

Pythonでの実装例

それでは、実際のコードを見てみましょう。

def solve(p, q, r, k):
   if r == 0:
      return 1
   elif r == 1:
      return (p % k, q % k)
   elif r % 2 == 0:
      return solve((p*p - q*q) % k, 2*p*q % k, r//2, k)
   else:
      (pr, qr) = solve(p, q, r-1, k)
      return ((p * pr - q * qr) % k, (p * qr + q * pr) % k)

print(solve(3, 0, 8, 10000))

入力

3, 0, 8, 10000

出力

(6561, 0)

出力の (6561, 0) は、3⁸ = 6561 であり虚部が 0 であることを示しています。期待どおりの結果ですね。

まとめ

ロシア農民乗算法の「指数を半分にしながら進める」という発想は、複素数のべき乗にもそのまま応用できます。また、途中で剰余を取りながら計算することで巨大な数による桁あふれを防げるため、競技プログラミングなどでも役立つ実用的なテクニックです。