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

# どのくらい速いの？
素数判定法といえば、 *エラトステネスの篩* を思い浮かべると思います。
与えられた数字が $N$ の時、 $\sqrt{N}$ の長さの配列を用意して、最も小さい素数から試しに割る方法なので、時間計算量も、空間計算量も、めちゃくちゃ食います。
素数の個数を扱うので、厳密にオーダーは計れませんが、割り算の回数を $f(N)$ とすると、
```math
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})$[^0000] あたりではないかと。
$N=10^6$ 辺りから不都合が出るのではないかと思います。（憶測ばかりですね）

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


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

競プロサイト AtCoder の[ある問題](https://atcoder.jp/contests/arc143/tasks/arc143_b "AtCoder Regular Contest 143 - B - Counting Grids")を解けず、解説を見たところ、後回しにしていたモジュロ演算（余剰演算）の実装がやっぱり必要だなと感じました。
競プロやっている限り、二項係数からは逃げられませんよね。

自作クラスの参考に、[AtCoder-Library](https://github.com/universato/ac-library-rb/blob/main/lib/modint.rb "Github(universato/ac-library-rb/lib/modint.rb)") を読んでいた[^0030]時、モジュロクラスの特異メソッド、 prime? メソッドで知らない挙動をするアルゴリズムを目にしました。

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

**コード内にある「2 7 61」をキーワードとして検索をかけてみれば、[ミラー-ラビン素数判定法](https://ja.wikipedia.org/wiki/%E3%83%9F%E3%83%A9%E3%83%BC%E2%80%93%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95 "ミラー-ラビン素数判定法 - Wikipedia")という単語にありついたわけです。**（不思議）

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

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

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

とりあえず、私自身がわかりやすいように書いたコードを載せておきます。
```ruby: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.判別条件

このアルゴリズムは、奇素数の特徴を利用して判別しています。よって、奇素数になり得ない数は直ちに判別してしまいましょう。
```Ruby
return false if n <= 1  # 1以下は非素数とする
return true if n == 2  # 奇素数のみの判定なので、2は先に素数判定とする
return false if n.even?  # 偶数は合成数
```
まず、１は素数ではありません。素数の定義は、約数に自身と１のみを持つ自然数とあります。
１が素数だと、素因数分解で $1^ \infty$ の項が出てきていまいます。
負の数を含め、１以下を非素数としましょう。

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

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

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

## 2.許容範囲

[Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%9F%E3%83%A9%E3%83%BC%E2%80%93%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95 "ミラー-ラビン素数判定法 - Wikipedia")には、以下の記述があります[^2000]。
> もし n < 4,759,123,141 なら、a = 2, 7, 61 について調べればよい。
> もし n < 341,550,071,728,321 なら、a = 2, 3, 5, 7, 11, 13, 17 について調べれば十分である。

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

```Ruby
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)$
この数も、なぜ必要になるかは後ほど。
```Ruby
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$ の項数と等しいです。

<details><summary>"(d & -d)"についてはこちら</summary>

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

とりあえず、 $d = 1011100000_{(2)}$ と例えてみましょう。
対して、 $-d$ を $d$ の基数の補数とするならば、
```math
\begin{eqnarray}
-d &=& 0100011111_{(2)}+1 \\
   &=& 0100100000_{(2)} \\
d &=& 1011100000_{(2)}
\end{eqnarray}
```
なんということでしょう。下６桁が綺麗に全く同じです。
１と０がひっくり返った上で、１を足したときに繰り上がりを繰り返す手順を行ったので、下の桁から０が連なる桁数を数えることができます。
仕上げに、下記を記す。
```math
d \land -d = 0000100000_{(2)}
```
（この`(d & -d)`表現も、AtCoder-Libraryのフェニック木で学んだ…）
（というか、どこかのフェニック木の記事を読ませた方がいい気がする）

---

</details>

もう、単純に奇数になるまで割り続けてもいい↓
```Ruby
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 を返す。

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

|          | 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.許容範囲](https://qiita.com/drafts/0310a959da872a7f890e/edit#2許容範囲)で用意した配列`[2,7,61]`の要素を $a$ にを試せば充分とされている。
また、[3.(N-1)を2で割り続けた数を求める](https://qiita.com/drafts/0310a959da872a7f890e/edit#3n-1を2で割り続けた数を求める)で用意した奇数 $d$ は、上記の説明の $d$ に当てはまる。

以上を踏まえて、
```Ruby
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](https://ja.wikipedia.org/wiki/%E3%83%9F%E3%83%A9%E3%83%BC%E2%80%93%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95 "ミラー-ラビン素数判定法 - Wikipedia") のコードも、 [Library](https://github.com/universato/ac-library-rb/blob/main/lib/modint.rb "Github(universato/ac-library-rb/lib/modint.rb)") のコードもどちらも`１`が出現したら、その時点で合成数かどうかを判断しています。
$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}$
1. (1.)の中でも、$a^{{2^r}d} \equiv -1 \pmod{n}$ （ $r$ は、$0$ 〜 $n-1$の整数）
1. または、$a^{d} \equiv 1 \pmod{n}$

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

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

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

# 作る上で 「二分累乗法」 は特に重要
実は、Ruby内では二分累乗法が組み込まれています。[^5000]
言語によっては、自分で実装しなければならないので、一応ここで短く説明しておきます。

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

```math
a^{233} = a^{1} \cdot a^{8} \cdot a^{32} \cdot a^{64} \cdot a^{128}
```
累乗というのは、ただ愚直に掛け続ければ、膨大に時間をかけてしまう計算です。
しかし、余剰演算の上では、上のようにたった数回掛けるだけで済んでしまいます。

# まとめ

1. [「3以上の奇素数」のみに限定する。](https://qiita.com/drafts/0310a959da872a7f890e/edit#1判別条件)
1. [計算可能範囲を意識する。](https://qiita.com/drafts/0310a959da872a7f890e/edit#2許容範囲)（32bitの正整数ならば、気にしなくていい。）
と同時に、サンプルのための配列を用意する。
1. [(n-1)を奇数になるまで割り続けた数 d を用意する。](https://qiita.com/drafts/0310a959da872a7f890e/edit#3n-1を2で割り続けた数を求める)
1. [サンプル全てで条件が通れば](https://qiita.com/drafts/0310a959da872a7f890e/edit#4リストとd乗その２乗)、それは素数である。

```ruby: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も公開できるように書き換えないとなぁ。）

[^0000]: 単調増加のため、$\Theta$記法にしています。
[^0010]: 四則演算等のオーダーは無視しています。実際には、掛け算のオーダーを加え、$O((\log(N))^3)$ほどになります。64bitを超える掛け算なんて、常識外ですからね。無視します。
[^0020]: $log(N)$は簡単にいえば、桁数を表しているようなものだと思ってください。１億なら、9桁ですよね。（ $log_{10}(100,000,000)=8$ ）
[^0030]: AtCoder-Library の存在を最近知ったばかりなんです。。
[^2000]: 論文とかではなく、ごめんなさい。でも、事実を確認できれば問題ないと思うので。
[^4000]: この表を書くためにこの記事書いた！
[^5000]: 意味わからん。セルフで実装した私がバカみたいじゃないか。