前提
この記事はNim 1.6.14でC++(gcc)にトランスコンパイルすることを前提に書かれています。
modpowって？
aをn乗したものをmで割った余りを求めるものです。繰り返し二乗法というアルゴリズムを用いることで、$O(\log(n))$で計算ができます。
計算途中の結果が最大で$(m-1)^2$となるので、mが32bit整数を超えるような値の場合はオーバーフローに気を付けて実装する必要があります。
作り方
C++の128bit整数を呼ぶことで、簡単に作成することができます。
proc mul128(a, b, m: int): int {.importcpp: "(__int128)(#) * (#) % (#)", nodecl.}
proc powmod(a, n, m: int): int =
result = 1 mod m
var
a = a
n = n
while n > 0:
if n mod 2 != 0: result = mul128(result, a, m)
if n > 1: a = mul128(a, a, m)
n = n shr 1
return result
使用例
64bit整数について確定的に判定できるミラーラビン素数判定を作成するときに、内部で用いる。
proc mul128(a, b, m: int): int {.importcpp: "(__int128)(#) * (#) % (#)", nodecl.}
proc powmod(a, n, m: int): int =
result = 1
var
a = a
n = n
while n > 0:
if n mod 2 != 0: result = mul128(result, a, m)
if n > 1: a = mul128(a, a, m)
n = n shr 1
proc isprime(N: int): bool =
let bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
if N == 2:
return true
if N < 2 or (N and 1) == 0:
return false
let N1 = N-1
var d = N1
var s = 0
while (d and 1) == 0:
d = d shr 1
s += 1
for a in bases:
var t: int
if a mod N == 0:
continue
t = powmod(a, d, N)
if t == 1 or t == N1:
continue
block test:
for _ in 0..<(s-1):
t = powmod(t, 2, N)
if t == N1:
break test
return false
return true