$a \times b \pmod{n}$の答えを高速に求めるアルゴリズム。$a,b,n$は整数、$0 \le a,b \lt n$。数字が大きくなっても平気。
上記の計算のどこに時間かかるって、コンピュータは、加算は得意、乗算もそこそこできるけど、除算は四則演算の中では苦手。モンゴメリ乗算は、mod nの部分を工夫して、除算を1回のみで行えるようにしてしまう。
1回のみだったら変わらないじゃん、となるのだが、$a \times b \times c \times d ... \pmod{n}$と複数回かけても1回の除算で済む。このような場合、正攻法では、数字が大きくなるとオーバーフローの問題や、桁数に比例して計算時間がかかるため、現実的には$((a \times b \bmod n) \times c \bmod n) \times d \bmod n ...$と、乗算毎に割ることになる。それが1回で済むようになるため、効果が大きい。
概要としては、うまく式変形することで、$\bmod n$の式から$n$による除算を取り除き、任意の定数$R$による除算と剰余（と、その他除算以外の演算）に置き換えてしまう。この$R$は任意に選べるので、コンピュータの得意な2の冪乗数にすることで、実質的にビット演算で行えるようになる。
我々も割り算は苦手でも、$123456789 \div 1000 = 123456 \ \text{あまり} \ 789$のように、割る数が10,100,1000…なら数字を区切るだけで求められるのと同じ。
詳細な説明はWikiとか他のサイトに譲るとして、pythonのコード。冪剰余も求められる。
class Montgomery:
def __init__(self, n):
self.n = n
self.nb = n.bit_length()
# Rを、Nより大きい最小の2の冪乗数とする
# R^2 mod n : この1回だけ除算が必要になる
self.r2 = (1 << (self.nb * 2)) % n
# Rを2の冪乗とすることで、mod Rをビットマスクで求められるようになる
self.mask = (1 << self.nb) - 1
# N * N' = -1 mod R となるN'の導出
# Rを2の冪乗とすることで加算とビットシフトで求められるようになる
self.nr = 0
t = 0
vi = 1
for _ in range(self.nb):
if t & 1 == 0:
t += n
self.nr += vi
t >>= 1
vi <<= 1
def reduction(self, t):
"""モンゴメリリダクション"""
c = t * self.nr
c &= self.mask
c *= self.n
c += t
c >>= self.nb
if c >= self.n:
c -= self.n
return c
def mul(self, a, b):
"""a * b mod n を計算"""
return self.reduction(self.reduction(a * b) * self.r2)
def exp(self, a, b):
"""a ^ b mod n を計算"""
p = self.reduction(a * self.r2)
x = self.reduction(self.r2)
y = b
while y:
if y & 1:
x = self.reduction(x * p)
p = self.reduction(p * p)
y >>= 1
return self.reduction(x)
pythonには、$a^b \pmod{n}$の計算ならpow(a,b,n)があるから、わざわざ自分で書く必要性は無いんだけどね。組み込みな分、圧倒的に速いし。
import timeit
mg = Montgomery(1000000007)
print(timeit.timeit('mg.exp(123456789, 987654321)', number=1000, globals=globals()))
# => 0.0361
print(timeit.timeit('pow(123456789, 987654321, 1000000007)', number=1000, globals=globals()))
# => 0.0037