LoginSignup
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@Shirona

めちゃ速い！ミラー-ラビン素数判定法

Last updated at Posted at 2022-07-17

どのくらい速いの？

素数判定法といえば、 エラトステネスの篩 を思い浮かべると思います。
与えられた数字が $N$ の時、 $\sqrt{N}$ の長さの配列を用意して、最も小さい素数から試しに割る方法なので、時間計算量も、空間計算量も、めちゃくちゃ食います。
素数の個数を扱うので、厳密にオーダーは計れませんが、割り算の回数を $f(N)$ とすると、

N \cdot \sum_{a=1} ^ {\log{\left( \sqrt{N} \right)}} \left(\frac{1}{2^a} \right) < f(N) <
N \cdot \sum_{a=1} ^ {\sqrt{N}} \left( \frac{1}{a} \right)　

だと思います。
なので、割り算自体のオーダーを無視すれば、 $\Theta(N^{3/2})$1 あたりではないかと。
$N=10^6$ 辺りから不都合が出るのではないかと思います。（憶測ばかりですね）

それに比べ、ミラー-ラビン素数判定法のオーダーは、$O(\log(N))$です23。
とても早いです。初めて見た時は、あまりの計算の少なさに疑いを持ち、この判定法の名前にたどり着いたくらいですから。

アルゴリズムを見つけた動機

競プロサイト AtCoder のある問題を解けず、解説を見たところ、後回しにしていたモジュロ演算（余剰演算）の実装がやっぱり必要だなと感じました。
競プロやっている限り、二項係数からは逃げられませんよね。

自作クラスの参考に、AtCoder-Library を読んでいた4時、モジュロクラスの特異メソッド、 prime? メソッドで知らない挙動をするアルゴリズムを目にしました。

Rubyist でかつ、ミラー-ラビン素数判定法を知らない方は１分だけでも眺めてみてください。
（他の言語でも AtCoder-Library に適用されているかもしれません。）

コード内にある「2 7 61」をキーワードとして検索をかけてみれば、ミラー-ラビン素数判定法という単語にありついたわけです。（不思議）

Rubyには、Primeクラスがありますが、一つの数が素数かどうか判別するという点に関しては、より高速になります。少なくとも、私が作ったコードでは。

というわけで、こういった手順を踏まえれば、「自分でもメソッドが作れるよ！」と誰でも思えるように、記事を書いていきたいと思います!

ミラー-ラビン素数判定法の手順

とりあえず、私自身がわかりやすいように書いたコードを載せておきます。

prime?.rb
def prime?(n)
  return false if n <= 1  # 1以下は非素数とする
  return true if n == 2  # 奇素数のみの判定なので、2は先に素数判定とする
  return false if n.even?  # 偶数は合成数
  
  list = n < 4759123141 ? [2, 7, 61].select{|x| x<n} :
         n < 341550071728321 ? [2, 3, 5, 7, 11, 13, 17] :
         [2,3,5,7,11,13,17,19,23,29,31].inject((3..1000).to_a){|r,i|r.select{|x|x% i!=0||x==i}}
  d = n-1  # # XXX1 => XXX0
  zero_len = (d & -d).bit_length-1
  d >>= zero_len  # XX100..0 => XX1
  list.each do |a|
    y = a.pow(d, n)  # a ** d % n (二分累乗法)
    y = n-1 if y==1  # a**d(奇数乗)でy==1の時、y==n-1に統合する
    (zero_len-1).times do  # d..(n-1)/4
      break if y == n-1
      y = y ** 2 % n  # a**((2**r)*d)
    end
    return false if y != n-1
  end
  return true
end

これで本当に走るのか？と思ったら、Primeクラスの結果と比べてみましょう。
（想定桁以内で例外が出たら、すみません。）

1.判別条件

このアルゴリズムは、奇素数の特徴を利用して判別しています。よって、奇素数になり得ない数は直ちに判別してしまいましょう。

return false if n <= 1  # 1以下は非素数とする
return true if n == 2  # 奇素数のみの判定なので、2は先に素数判定とする
return false if n.even?  # 偶数は合成数

まず、１は素数ではありません。素数の定義は、約数に自身と１のみを持つ自然数とあります。
１が素数だと、素因数分解で $1^ \infty$ の項が出てきていまいます。
負の数を含め、１以下を非素数としましょう。

次に、2は最も小さい素数です。trueを返しましょう。

残りの偶数は、すべて合成数です。falseを返しましょう。

あとは、３以上の奇数のみになりました。

2.許容範囲

Wikipediaには、以下の記述があります5。

もし n < 4,759,123,141 なら、a = 2, 7, 61 について調べればよい。
もし n < 341,550,071,728,321 なら、a = 2, 3, 5, 7, 11, 13, 17 について調べれば十分である。

というわけで、こういった配列を用意します。
何に使う配列かは、後ほど。

list = n < 4759123141 ? [2, 7, 61].select{|x| x<n} :
       n < 341550071728321 ? [2, 3, 5, 7, 11, 13, 17] :
       [2,3,5,7,11,13,17,19,23,29,31].inject((3..1000).to_a){|r,i|r.select{|x|x% i!=0||x==i}}

実は、ミラー-ラビン素数判定法は、確率的アルゴリズムに属するもので、この配列の情報が少なすぎると、"合成数であるのに、素数と判定してしまった"というケースが発生します。※かなり大きい数字に対しての濫用は禁物です。
実際に、4759123141は48781と97561を約数に持つ合成数です。にもかかわらず、ここでlistに[2, 7, 61]を渡すと、trueを返してしまいます。
しかし、正整数のみ扱う32bitの最大数でも 4,294,967,296 。 4,759,123,141 以下です。滅多に問題には出現しないので、[2, 7, 61]のみの準備で構わないでしょう。

それに、忘れるところでした。一行目の.select{|x| x<n}は重要です。用意する要素は判別する要素より小さくなければなりません。

３行目は、ロマンですね。限度がいくつかわかりません。数学者の仕事です。

3.(N-1)を2で割り続けた数を求める

後に必要になるのは、$d,2d,4d,8d, \dots (n-1)/8,(n-1)/4,(n-1)/2$ になります。
上記は、 $d$ のみ奇数で、残りは偶数です。$d \cdot 2^r$ $(0 \leq r \leq zero\_len-1)$
この数も、なぜ必要になるかは後ほど。

d = n-1  # # XXX1 => XXX0
zero_len = (d & -d).bit_length-1  # d:XXXX100..0 => (d & -d):100..0
d >>= zero_len  # XX100..0 => XX1

計算数１っぽい表現が好きなので、少し読みにくいですね。２進数表記で下桁から１が出現するまでの０の桁数をzero_lenに格納しています。
zero_lenは、$d,2d,4d,8d, \dots (n-1)/8,(n-1)/4,(n-1)/2$ の項数と等しいです。

"(d & -d)"についてはこちら

(d & -d)がどのような数になるかを説明するためには、ビット論理積と基数の補数についての前提知識が必要になります。それぞれ一つずつの記事になってしまうので、割愛。

とりあえず、 $d = 1011100000_{(2)}$ と例えてみましょう。
対して、 $-d$ を $d$ の基数の補数とするならば、

\begin{eqnarray}
-d &=& 0100011111_{(2)}+1 \\
   &=& 0100100000_{(2)} \\
d &=& 1011100000_{(2)}
\end{eqnarray}

なんということでしょう。下６桁が綺麗に全く同じです。
１と０がひっくり返った上で、１を足したときに繰り上がりを繰り返す手順を行ったので、下の桁から０が連なる桁数を数えることができます。
仕上げに、下記を記す。

d \land -d = 0000100000_{(2)}

（この(d & -d)表現も、AtCoder-Libraryのフェニック木で学んだ…）
（というか、どこかのフェニック木の記事を読ませた方がいい気がする）

もう、単純に奇数になるまで割り続けてもいい↓

d = n-1  # # XXX1 => XXX0
d >>= 1 while d[0]==0

4.リストとd乗、その２乗

資料を読んでいて、最も躓くのがこの工程。
証明に目を取られて、どの部分が実用的なのかも見つかりづらいと思いました。

$a^d \ne 1\ \bmod\ n$ で、かつ $[0, s - 1\ ]$ の範囲の全ての $r$ について $a^{{2^r}d}\ \bmod\ n \ne\ -1$ なら、composite を返す。

対偶の条件を視覚的にわかるようにすれば、こんな感じ。6

1 (mod n) ... n-1 (mod n)
$a^d$ ○ ○
$a^{2d}$ ○
$a^{4d}$ ○
: : : :
$a^{\frac{n-1}{4}}$ ○
$a^{\frac{n-1}{2}}$ ○
$a^{n-1}$

用意した $a$ 全てに対して、それぞれ○に当てはまる数が少なくとも一つあれば、次の $a$ を走査する手順となる。
$1 \leq a \leq (n-1)$ を満たす全ての正整数 $a$ において上記に当てはまるならば、それは素数と言えるとなっている。
しかし、$n < 4759123141$ の範囲ならば、2.許容範囲で用意した配列[2,7,61]の要素を $a$ にを試せば充分とされている。
また、3.(N-1)を2で割り続けた数を求めるで用意した奇数 $d$ は、上記の説明の $d$ に当てはまる。

以上を踏まえて、

list.each do |a|
  y = a.pow(d, n)  # a ** d % n (二分累乗法)
  y = n-1 if y==1  # a**d(奇数乗)でy==1の時、y==n-1に統合する
  (zero_len-1).times do  # d..(n-1)/4
    break if y == n-1
    y = y ** 2 % n  # a**((2**r)*d)
  end
  return false if y != n-1
end
return true

と、私は書きました。

しかし、 Wikipedia のコードも、 Library のコードもどちらも１が出現したら、その時点で合成数かどうかを判断しています。
$a^{n-1} \equiv 1 \pmod{n}$ かどうか調べるだけではいけないのかと思ったりもしました。結局、証明にまで手を回さねばならなくなる。。

頭が混乱したので、一度定理を整理する。

  • $n$が素数であり、かつ $a \not\equiv 0 \pmod{n}$ $\Rightarrow$ $a^{n-1} \equiv 1 \pmod{n}$ （フェルマーの小定理より）
  • $n$が素数であり、かつ $a^2 \equiv 1 \pmod{n}$ $\Rightarrow$ $\pmod{n}$ において、 $a \equiv 1$ または、 $a \equiv -1$
  • $a^2 \equiv 1 \pmod{n}$ かつ $a\not\equiv 1,-1\pmod{n}$ $\Rightarrow$ $n$ は合成数

残るイレギュラーの可能性:

  1. 合成数でも、$a^{n-1} \equiv 1 \pmod{n}$
  2. (1.)の中でも、$a^{{2^r}d} \equiv -1 \pmod{n}$ （ $r$ は、$0$ 〜 $n-1$の整数）
  3. または、$a^{d} \equiv 1 \pmod{n}$

ちなみに、下記は $p$ を奇素数として表現してややこしい流れになりますが、

二乗して $1\pmod{p}$となるのは、$1$ および $-1$ だけということになる。

らしいので、計算途中-1が出現する前に1が出現すれば、判断する数は合成数ということになります。
よって、参照先のコードは特に問題ないことになります。

作る上で 「二分累乗法」 は特に重要

実は、Ruby内では二分累乗法が組み込まれています。7
言語によっては、自分で実装しなければならないので、一応ここで短く説明しておきます。

以下の式を見れば、実装はすぐにできるでしょう。

a^{233} = a^{1} \cdot a^{8} \cdot a^{32} \cdot a^{64} \cdot a^{128}

累乗というのは、ただ愚直に掛け続ければ、膨大に時間をかけてしまう計算です。
しかし、余剰演算の上では、上のようにたった数回掛けるだけで済んでしまいます。

まとめ

  1. 「3以上の奇素数」のみに限定する。
  2. 計算可能範囲を意識する。（32bitの正整数ならば、気にしなくていい。）
    と同時に、サンプルのための配列を用意する。
  3. (n-1)を奇数になるまで割り続けた数 d を用意する。
  4. サンプル全てで条件が通れば、それは素数である。
prime?.rb
def prime?(n)
  # 1. 「3以上の奇素数」のみに限定する。
  return false if n <= 1  # 1以下は非素数とする
  return true if n == 2  # 奇素数のみの判定なので、2は先に素数判定とする
  return false if n.even?  # 偶数は合成数
  
  # 2. 計算可能範囲を意識する。
  #    と同時に、サンプルのための配列を用意する。
  list = n < 4759123141 ? [2, 7, 61].select{|x| x<n} :
         n < 341550071728321 ? [2, 3, 5, 7, 11, 13, 17] :
         [2,3,5,7,11,13,17,19,23,29,31].inject((3..1000).to_a){|r,i|r.select{|x|x% i!=0||x==i}}
  
  # 3. (n-1)を奇数になるまで割り続けた数 d を用意する。
  d = n-1  # # XXX1 => XXX0
  zero_len = (d & -d).bit_length-1
  d >>= zero_len  # XX100..0 => XX1
  
  # 4. サンプル全てで条件が通れば、それは素数である。
  list.each do |a|
    y = a.pow(d, n)  # a ** d % n (二分累乗法)
    y = n-1 if y==1  # a**d(奇数乗)でy==1の時、y==n-1に統合する
    (zero_len-1).times do  # d..(n-1)/4
      break if y == n-1
      y = y ** 2 % n  # a**((2**r)*d)
    end
    return false if y != n-1
  end
  return true
end

最後までご閲覧いただき、ありがとうございます！
（Git-Hubも公開できるように書き換えないとなぁ。）

  1. 単調増加のため、$\Theta$記法にしています。 ↩

  2. 四則演算等のオーダーは無視しています。実際には、掛け算のオーダーを加え、$O((\log(N))^3)$ほどになります。64bitを超える掛け算なんて、常識外ですからね。無視します。 ↩

  3. $log(N)$は簡単にいえば、桁数を表しているようなものだと思ってください。１億なら、9桁ですよね。（ $log_{10}(100,000,000)=8$ ） ↩

  4. AtCoder-Library の存在を最近知ったばかりなんです。。 ↩

  5. 論文とかではなく、ごめんなさい。でも、事実を確認できれば問題ないと思うので。 ↩

  6. この表を書くためにこの記事書いた！ ↩

  7. 意味わからん。セルフで実装した私がバカみたいじゃないか。 ↩

1
0
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?