Julia の素数周りの実装

## 概要

Julia には，素数周りの関数が標準で組み込まれています．その実装にどんなアルゴリズムが採用されているのか調べてみました．

いきなりまとめですが，Julia の素数周りのアルゴリズムでは，

* 素数生成: <del>Atkin の篩</del> Wheel factorization を使った Eratosthenes の篩 [更新]
* 素数判定: Miller-Rabin 素数判定法
* 素因数分解: 試し割り法 + Pollard-Brent Rho 法

が採用されています．以下，それぞれについて軽く紹介します．ちなみに，この記事では，アルゴリズムの詳細や，計算量の説明はしません．誰かの Julia のコードを調査するとかっかかりにでもなったら幸いです．

この記事は，2015-07-13 時点での [コード](https://github.com/JuliaLang/julia/blob/b2528a6e66/base/primes.jl) をもとに書きました．動作確認にも，`b2528a6e66 ` の Nightly ビルドを使っています．

## 素数生成 = <del>Atkin の篩</del> Wheel factorization を使った Segment 版の Eratosthenes の篩 [更新]

素数生成アルゴリズムとして一番有名なのは，[Eratosthenes の篩](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes)というアルゴリズムでしょう．このアルゴリズムは単純なわりに動作速度が実用的で，1000 万くらいまでの素数ならばお手軽に求められます．Project Euler の問題を解くときなども Erathosthenes の篩を自前実装する人が多いようです．


<del>Julia で使われているのは，Eratosthenes とは若干趣の異なる [Atkin の篩](https://en.wikipedia.org/wiki/Sieve_of_Atkin) です．Atkin の篩は，計算量が Eratosthenes の篩よりも少ないことが知られています．コードは[このあたり](https://github.com/JuliaLang/julia/blob/d05a5790b6b06c05562295c155f69f736b2469eb/base/primes.jl#L9)です．
</del>

[2015-08-10 に更新．素数生成で使われているアルゴリズムが変更になりました]

Julia で使われているのは，Eratosthenes の篩を [Wheel factorization](https://en.wikipedia.org/wiki/Wheel_factorization) を用いて高速化したものになっています．もともとは，Atkin の篩というアルゴリズムが使われていましたが，Eratosthenes の篩を高速化したもののほうが実験ベースで実行速度が高速なことが知られている（[primesive](http://primesieve.org/)などの実装による）ため，Eratosthenes の篩に舵を切ったのだと思います．

さて，この素数生成の [PR](https://github.com/JuliaLang/julia/pull/12025) でも触れられているように，segmented sieve という手法を用いて CPU キャッシュサイズに最適化した実装をすることにより，より高速化することが可能なことが知られています．こちらも実装されるかもしれません．


試しに 10000000までの素数を求めてみます．

```
julia> gc; @time length(primes(10000000))
  0.056016 seconds (10 allocations: 10.543 MB, 6.05% gc time)
664579
```

僕の場合 10000000 あれば十分なことが多く，この速度は十分実用的です．


## 素数判定 = Miller Rabin 素数判定法

ある整数 n が素数かどうかは，n を n の平方根まで割ってみれば良いですね．これを「試し割り法」といいます．この方法はそれなりに大きい整数になってくると実用的ではありません．

ある整数 n が素数かどうか，素因数をがんばって計算しなくても「素数かどうか」という判定だけを高速に行う方法がないかな，と考えた人がいて，いろいろアルゴリズムが発見されています．[Fermat 法](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A7%E3%83%AB%E3%83%9E%E3%83%BC%E3%81%AE%E5%B0%8F%E5%AE%9A%E7%90%86#.E3.83.95.E3.82.A7.E3.83.AB.E3.83.9E.E3.83.BC.E3.83.86.E3.82.B9.E3.83.88)と呼ばれるものが比較的素朴な，確率的素数判定法です．これは Fermat の小定理という定理に基づいたアルゴリズムです．Fermat 法は，合成数を素数だよと誤判定してしまうことがあります，Carmichael 数はそのような素数のひとつです．

Fermat 法より良い確率的素数判定アルゴリズムに [Miller-Rabin 法](https://ja.wikipedia.org/wiki/%E3%83%9F%E3%83%A9%E3%83%BC-%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95)があります．こちらは Carmichael 数も合成数と正しく判定できます．この方法では，十分な試行を重ねることによって，誤判定の確率を無視できるほど小さくできることが知られています（たぶん）

多くの素数ライブラリに，Miller-Rabin 法が利用されていて，Julia もそうです．

使ってみましょう．

```
julia> gc; @time isprime(280671392065546467397265294532969672241810318954163887187279320454220348884327)
  0.005004 seconds (1.28 k allocations: 65.417 KB)
false

julia> gc; @time isprime(67280421310721)
  0.086759 seconds (43.43 k allocations: 1.921 MB)
true
```

こちらも，十分な速度で動作します．

確率的素数判定法には，ほかには Baillie–PSW というのもあります．こちらは Miller-Rabin よりも精度が高いがちょっと遅い，みたいなアルゴリズムで，Scheme 実装の Gauche などではこちらが用いられているようです．

Fermat 法も Miller-Rabin 法も，確率的素数判定法なのですが，AKS 素数判定法という，「多項式時間で動作」し，かつ「決定的」に素数判定を行うアルゴリズムもあります．しかし，多項式時間といえどあまり実用的な速度では動かないので，ライブラリなどで使われているのをみたことはありません．

## 素因数分解 = 試し割り法 + Pollard-Brent Rho 法

素数周りのアルゴリズムで最も盛んに実装が行われているものといえば，素因数分解アルゴリズムかと思います．古来より地道な研究が行われ，数多のアルゴリズムが作られています．小さめの（1時間で分解が完了する程度の）整数を分解するための実装も含めた詳しい話が[こちらのスライド](http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf) にありますので，参考にしてみてください．

Julia の素因数分解アルゴリズムは，小さい整数に対しては「試し割り法」，それなりに大きい整数に関しては，[Pollard-Brent Rho 法というものが使われています](https://github.com/JuliaLang/julia/blob/d05a5790b6b06c05562295c155f69f736b2469eb/base/primes.jl#L105)．ちなみに，なぜ Rho 法という名前なのかというと，このアルゴリズムを可視化したときに見える図が，[ρの形に見えるから](https://www.google.co.jp/search?q=algorithm+introduction+pollard+rho&espv=2&biw=1307&bih=1105&source=lnms&tbm=isch&sa=X&ei=53ujVZDII4rN0ATKvrOwAQ&ved=0CAgQ_AUoAw#tbm=isch&q=pollard+rho)，ということみたいです．

ちなみに，もともと Pollard Rho 法というのがあって，それを高速化したのが Pollard-Brent Rho 法です．

Pollard-Brent Rho 法は，素因数分解対象の合成数が比較的小さい素因数を持っている場合に威力を発揮します．なので，小さな素因数を持たない場合は別の，例えば SQUFOF とか MPQSを使ってみるとかが考えられますが，Pollard-Brent Rho 法に比べると，実装は若干面倒くさくなります．

まぁ細かい話は抜きにして，使えればいい！ というわけで，使ってみましょう．

```
julia> gc; @time factor(18446743979220271189)
  1.027358 seconds (4.24 M allocations: 80.591 MB, 14.16% gc time)
Dict{Int128,Int64} with 2 entries:
  4294967279 => 1
  4294967291 => 1
  
julia> gc; @time factor(280671392065546467397265294532969672241810318954163887187279320454220348884327)
  1.289839 seconds (5.41 M allocations: 111.762 MB, 14.55% gc time)
Dict{Base.GMP.BigInt,Int64} with 7 entries:
  369941863          => 2
  481362815814826159 => 1
  162425297          => 1
  479871607          => 1
  215940091          => 1
  358456949          => 1
  706170617          => 1
```

思ったより速かった．というのも，ここでは因数が 10 進 9 桁のものを使っているからですが…．

## 感想

Julia の素数周りの実装は，手堅い構成になっていると思います（うわ，偉そうだ）．言語組み込みの素数周りの実装としては十分といえるのではないでしょうか．といっても，もちろん Msieve とか，専用の素因数分解ユーティリティには太刀打ちできなことは確実なので，Julia 組み込みで無理そうだったら，やはり専用のユーティリティに頼るべきでしょうね．

いつしか数体篩法を Julia で実装し，4 桁 bit の素因数分解を達成する猛者が現れることを夢見つつ，終わりたいと思います．