この記事では素数判定アルゴリズムである”ミラーラビン法”をpythonで実装する方法について解説します。ミラーラビン法は”ミラーラビンテスト”と言われることもあります。
ミラーラビン法は、素数かどうかを判定したい数字pと素数判定の精度に関係するパラメータtを入力することで、「この数は素数である」もしくは「この数は素数ではない」という判定を出力します。判定の精度は100％ではなく、tの値によって変動します。
「素数判定の精度に関係するパラメータt」は1以上の整数を入力します。大きければ大きいほど判定の精度が上がります。
ミラーラビンテストが誤判定を出力する確率(素数ではない数を素数であると判定してしまう確率)は、1/4tである事が知られています。tを1大きくするごとに、誤判定の確立が1/4になります。
たとえばt=50とすれば、誤判定の確立は1/450以下となります。これは大体10-30くらいです。(1兆×1兆×100万)分の1くらいです。めちゃくちゃ小さい確率ですね。このくらい小さければ実用上は0とみなしてしまえます。なので「誤判定率は0じゃないの!?」と心配された方も、もう安心です!
しかも、1/4tというのは理論的な最悪値です。実際はそれよりもはるかに小さいことがほとんどであるといわれています。
細かい理屈はこの記事では書きません。実装の解説に注力しようと思います。
まずはどういったアルゴリズムでミラーラビンテストが動くのかを見てみましょう。
目次
ミラーラビン法のアルゴリズム
input: 判定したい数p,確率パラメータt
[step1] pが1の場合は「False」、2の場合は「True」と出力し終了する。
[Step2] p – 1 = 2skとなるようなs,kを求める(ただしkは奇数)。
[Step3] i = 1 とし、 以下を繰り返す。
[Step3-1] 1<a<pとなる整数aをランダムに生成する。
[Step3-2] b = ak mod p とする。
[Step3-3] j = 0としてj < sかつb ≢ ±1(mod ｐ)である限り、以下の処理を繰り返す。
[Step3-3-1] b = b2mod pを計算し、j +=1とする。
[Step3-4] j=s または (j>0 かつ b=1)なら「False」と出力して終了する。
[Step3-5] i = tなら「True」と出力して終了する。
[Step3-6] i += 1とする。
Step1に必要性について
[Step1]の処理は、pの入力を「3以上の整数」と限定すれば必要ありません。
このアルゴリズムだとpが1か2の時にバグります。[Step3-1]で乱数を生成する際に、例えばp=2だと「1<a<1のaを生成する」というプログラムとなってしまいます。これでは生成できるaが存在しませんね。
pが1か2だと生成できる乱数が存在しなくなってしまいます。こういったバグを発生させないために[Step1]の部分を書いておくと良いと思います。
step2を実現するための方法
[Step2]は次のような記述でしたね。
[Step2] p – 1 = 2skとなるようなs,kを求める(ただしkは奇数)。
この条件を満たすs,kのペアは一つとは限りません。条件さえ合っていれば何でもいいです。
これを実現する方法はいろいろありそうですが、一つは「とりあえずsは無視して考えて、kを決定してしまう」という方法があります。無視するとは、まずはs=1と置いてみる事です。kが決まったらそれに合うsを求めます。
まずk = (p-1)/2として、この時点でkが奇数ならそのままの値でkを決定します。もし偶数である場合は、奇数になるまで2で割り続けます。
kが決まればsは簡単に求められます。
2s = (p-1)/k ➡ s = log2｛(p-1)/k｝
という感じに変形して求めます。
pythonによる”ミラーラビン法”の実装
それでは上記のアルゴリズムをpython3によるコードで実装します。
先にコード全体を示す以下のようになります。
import math
import random
def MRPT(p,t):
"""Step1"""
if p == 1:
return False
if p == 2:
return True
"""Step2"""
k = (p-1)//2
while k%2 == 0:
k //= 2
s = int(math.log(((p-1)/k),2))
"""Step3"""
i = 1
while True:
"""Step3-1"""
a = random.randint(2,p-1)
"""Step3-2"""
b = pow(a, int(k), p)
"""Step3-3"""
j = 0
while ( j < s and b != 1 and b != p-1):
"""Step3-3-1"""
b = pow(b, 2 ,p)
j += 1
"""Step3-4"""
if (j == s or(j > 0 and b ==1)):
return False
"""Step3-5"""
if (i == t):
return True
"""Step3-6"""
i += 1
「ＭＲＰT」とはミラーラビン法を英語で書いた時の頭文字(Miller–Rabin primality test)です。
解説と注意点
いくつか解説や注意点を書いていきます。
インポートするモジュール
pow()関数とlog()関数を使うためにmathモジュールをインポートします。
乱数生成も行うので、randomモジュールをインポートします。
[Step2]
このステップはsとkを求めるプロセスです。先に解説したように、まずはkを求めています。
とりあえず「k = (p-1)/2」とかいてこれが奇数かどうかを調べます。pythonの場合while文のループに入る前にwhile分の条件あっているかどうかを調べられるため、while文に入る前にif文でチェックする必要はありません。
14行目でひたすらkを2で割り続け、kが奇数になるとwhileループから抜け出します。あとは、「s = int(math.log(((p-1)/k),2))」で「s = log2｛(p-1)/k｝」の計算をしてsの値を決定します。
[Step3]
「while True:」で無限ループとしています。ループの最後に「i+=1」として、「i=t」になれば終了するので、これで問題はありません。
[Step3-2]
pow関数は引数すべてが整数でなければバグります。この時点で実はkがfloat型となっているため、kをint型にしてやらなければなりません。
[Step２]の時点で「k = (p-1)/2」の部分(12行目)を「k =int((p-1)/2)」としてやってもよいです。もしくは除算記号ではなくシフト演算をつかって記述してやってもよいです。
シフト演算のほうがスマートで高速ですが、わかりやすさ重視で除算記号を使いました。除算命令はまあまあ重いので使用を避けると速くなります。
[Step3-3]
[Step3-3]で行うb ≢ ±1(mod p)の判定については「b != 1 and b != p-1」で表現しています。
[Step3-2],[Step3-3-1]より、b = ak mod pもしくは b = b2mod p ですから、b ≧ p-1 です(bはある値をpで割った余りなので、pを超えることはありえません)。したがって b ≡ ±1(mod p)となるのは、b = 1 か b = p-1 の時だけです。
したがって、「b != 1 and b != p-1」と書けばOKとなります。
ミラーラビン法のいいところ
ミラーラビン法のいいところは、tの値を選ぶことで誤判定率を把握して素数判定を行えることです。
誤判定率ゼロのアルゴリズムも存在しますが(有名なのはAKSアルゴリズム)、ミラーラビン法に比べると低速です。なのでこちらが使われることが多いです。
参考になりそうな他サイトへのリンク
この記事を書くにあたっていろいろ調べました。その中で見つけたWebサイトなどの中で、よさそうなページをメモしておきます。
・深川久『Miller-Rabin による素数の確率的判定法』
大学の資料か何かかな？かなりガッツリめに事が書いてあります。ミラーラビン法を支える定理とその証明など。理解するために必要な数学的な基礎知識についても書いてあります。
しっかり読めば勉強になりそう。今度読む。
・素数判定いろいろ – フェルマーテスト・ミラーラビン素数判定法
こちらはやさしい解説が書いてあります。pythonによるコード例もあります。私のこの記事でちょろっと述べた「kをシフト演算で求める」方法で書かれています。ミラーラビン法以外の素数判定法についての解説やpythonコードもあります。
私もフェルマーテストについての記事を書きたいね。
当ブログの関連記事
ここで紹介したプログラムを活用して素数を生成するプログラムを作ったよー、という記事です。ミラーラビン法の活用事例ですね。
こちらは既存の関数を使って素数判定をする方法の解説記事です。
お詫びと訂正(2020/07/25)
記事の中に誤りを見つけたので訂正しました。
この記事を作成した際に、なぜか2は素数ではないと勘違いしていたようで、そのせいでいくつか誤りがありました。
訂正前のプログラムだと、2を素数判定したときにFalseと返されます。実際はTrueと返されるのが正しいです。
訂正前に読んでくれた方にはご迷惑をおかけしました。申し訳ありませんでした。