F# で素数判定 (Miller-Rabin)

Miller-Rabin 素数判定法は、与えられた奇数が素数か否かを確率的に判定する手法。
合成数を間違って素数と判定してしまう可能性はあるが、素数を見逃すことはない。
非常に大きい数が素数かどうかを単発で判定することを高速に行える(連続して判定したい場合は、他の方法が効果的なことも)。


こんなのを何に使うのかは全く見当も付かないが、まあそれはともかくとして F# で書いてみた。
ポイントは、「出来るだけ F# らしく書いてみた(つもり)」というところ。

let inline square n = n * n
let rec modpowI (a:bigint) n m = 
    if n = 1I then a
    elif n % 2I = 0I then (modpowI a (n / 2I) m |> square) % m
    else ((modpowI a (n - 1I) m) * a) % m
let rec fac2pow s d = if d % 2I = 0I then fac2pow (s + 1) (d / 2I) else (s, d)

let miller_rabin k n =
    let n1 = n - 1I
    let s, d = fac2pow 0 n1
    let mx = if n1 > (bigint System.Int32.MaxValue) then System.Int32.MaxValue else int n1
    let rand = new System.Random()
    let f _ =
        let rec check ar r = ar = n1 || (r <> 1 && check (ar * ar % n) (r - 1))
        let ad = modpowI (bigint(rand.Next(1, mx))) d n
        ad = 1I || check ad s
    Seq.forall f {1..k}

let isprime n = if n <= 2I then n = 2I else n % 2I <> 0I && miller_rabin 3 n


printfn "%A" (isprime 1150921I)  // prime
printfn "%A" (isprime 1150922I)  // not prime
  • 小さい数ではわざわざ Miller-Rabin を使う意味があまりないので、bigint 限定。
  • 乱数を使うが、System.Random は符号付き int の範囲の乱数を返すだけ。まじめに作ろうかと迷いつつも、結局 int の最大値の範囲までの乱数を用いるようにした。
  • 確率的手法ゆえ、試行回数が増えれば正しい結果が得られる確率は上がる。iprime の末尾あたりの "3" がその試行回数。