今日のラーメン

f:id:nisemono_san:20160410115447j:plain

概要

素数を検査する方法には、いわゆる「ふるい法」と呼ばれる方法を使った決定的な方法と、「だいだいこのような式に該当するならば、素数だろう」といった推測に基づく、確率的な方法がある。確率的な方法の一つに、「フェルマーの小定理」を利用したテスト方法がある。しかし、これは「カールマイケル数」と呼ばれる合成数が引っかかる。

はじめに

暫らくの間、『はじめての数論』を読みながら、素数判定について調べていた。そして、素数チェックの方法として有名な「フェルマーテスト」が紹介されていた。そういえば、この「フェルマーテスト」について調べたことなかったな、と思ったので、改めて実装しながら理解したことをメモしたいと思う。

フェルマーの小定理

まず、「フェルマーテスト」の方法は、「フェルマーの小定理」から作られたものである。

「フェルマーの小定理」は次のようになる。素数になる任意の数(ここではpとしよう)があり、pの倍数ではない(あるいは割りきることがない)任意の整数aを定義する(式としてあらわすなら ${a \not\equiv 0 (mod p) }$ )。このとき、 ${a^{p-1} \equiv 1 (mod p)}$ になる。

逆に考えれば、 ${a^{p-1} \equiv 1 (mod p)}$ を満たすpは素数ではないか、と推測できる。「はじめに」で伸べたように、これは必ずしも素数ではないことがわかっている。しかし、 ${a^{p-1} \equiv 1 (mod p)}$ という式を満たす合成数のカールマイケル数をあらかじめ表として保持しておけば、実用に困ることはないたろう(SICPによれば、このような数は100000000以下においては255個だと言われている)。

しかし、実際に ${a^{p-1}}$ をそのまま計算するのは無茶である。例えば、50093は素数だが、 ${2^{50092}}$ はだいだい15000桁になるので、数を求めてから、その余りを計算するのは現実的ではない。

とはいえ、これは実際に数を出すときの話であって、今回の関心は「余り」にしかない。極論を言ってしまえば、[{tex:{2^{50092}}]がどんな数であろうとも、 ${1 (mod p)}$ を満たすかどうかにしか、興味がない。そこで、このような数に対して、どのように余りを見つけるかについて、考えてみることにしたい。

繰り返し自乗法(表ありバージョン)

さて、現在持っている情報としては、aのp - 1乗するということであり、そして推測としては、pが素数であるならば、余りが1になることであった。

まず最初に、 ${a^{p}}$ に対して余りを出す方法を考えてみよう。このとき、繰り返し自乗法という方法が使える。

この方法は、 ${a^{b+c} (mod p) \equiv a (mod p) * b (mod p)}$ という法則性に基づいている。例えば、 ${5^{3 + 2} (mod 7) \equiv 4 (mod 7) * 6 (mod 7) \equiv 3 (mod 7)}$ といったようにだ。この場合ならば小さい計算であるので、手でできるものだが、例えば ${7^{327} (mod 853) }$ だった場合はどうだろうか。

まず最初に、327が二進法で分解できることに着目し、それぞれの7の2k乗でどのような余りになるか調べてみよう。例えば、 ${ 7^{1} \equiv 7 (mod 853)}$ だし、 ${ 7^{2} \equiv 49 (mod 853) }$ だ。指数法則に基づき、 ${ 7^{4} \equiv (7^{2})^{2} \equiv 695 (mod 853) }$ といったように、順々に表ができる。

これをOCamlのコードで直したのが次のものになる:

let pow2modlist x n m =
  let rec pow2mod x n a m xs =
    let int_a = int_of_float a in
    if n < int_a then xs
    else
      let next_mod = int_of_float ((float_of_int x) ** 2.) mod m in
      pow2mod next_mod n (a *. 2.) m ([(int_a, x mod m)] @ xs)
  in pow2mod x n 1. m []

さて、このあと327を2のk乗の数で分解すると、256 + 64 + 4 + 2 + 1となる。前から順番に、余りを掛けていくといいということになる:

let mymodcalc_with_list x n m =
  let rec mycalc x n mlist =
    match mlist with
      [] -> x
    | (y, z) :: mlist ->
       if n > y then mycalc (x * z mod m) (n - y) mlist
       else mycalc x n mlist
  in if n == 0 then 1
     else if n == 1 && x < m then x
     else mycalc x n (pow2modlist x n m)

繰り返し自乗法(表無しバージョン)

さて、この方法は表を作らなければいけないために、非常に効率が悪い。そこで、表を省略したものが上記書に掲載されているので、それをRubyにしたものを掲載する。

def modcalc(a, k, m)
  b = 1
  while k >= 1
    b = (a * b)  % m if k.odd?
    a = (a ** 2) % m
    k = (k / 2).floor
  end
  return b
end

また、OCamlだと次のようになる:

let mymodcalc a k m =
  let is_odd x = (x mod 2) == 1 in
  let rec mycalc a k m b =
    let next_b = if is_odd k then a * b mod m else b in
    if k < 1 then b
    else mycalc (int_of_float (float_of_int a ** 2.) mod m) (k lsr 1) m next_b
  in
  mycalc a k m 1

さて、これによってフェルマーの小定理を利用した素数のテストができるようになる。

フェルマーの素数テスト

このエントリを書く前に、さらっと他のアルゴリズムを参考にしたのだが、いわゆる「フェルマーの小定理」のそのまま( ${a^{p-1} \equiv 1 (mod p)}$ )を利用する方法と、この式を少し変形したもの( ${a^{p} \equiv a (mod p)}$ )を利用するパターンに分かれている。後者の変形については、a < pの場合、 ${a^{p} \equiv a^{p-1}*a^{1} \equiv 1 * a \equiv a (mod p)}$ と変形できるためだろう。今回は、本書に従い、後者で実装する。

とはいえ、この素数テスト自体は非常に簡単で、ランダムにaに値する数を選びだし、上の関数を利用して、余りが同じになるかどうかをチェックするだけで良い。コードに直すと下のようになる。

let fermat_test x =
  let random_int = ref (Random.int (int_of_float (sqrt (float_of_int x))) + 2) in
  (mymodcalc !random_int x x) == !random_int

また、Rubyなら、下のように書ける:

def fermat_test x
  limit = Math.sqrt(x).ceil
  2.upto(limit) do |n|
    return false if (modcalc n, x, x) != n
  end
  true
end

テストを潜りぬける数としてのカールマイケル数

さて、この方法はSICPには「数学と工学の違い」と書いてあるように、このテストをすりぬける数というのが存在する。それが冒頭に觝れたカールマイケル数である。

このカールマイケル数の定義はとても単純で、 ${a^{p} \equiv a (mod p)}$ を満たす合成数という、そのままの定義が与えられている。このような数のリストを、いわば「ふるい法」と比較してみると、カールマイケル数がどのようなものであるかという一覧を取りだすことが可能である。以下がRubyで100000以下のカールマイケル数を出力した場合のリストである。

しかし、このように従来の「ふるい法」の素数リストと、「フェルマーテスト」で引っかかった数を比較するのは、あまりにも冗長である。カールマイケル数にはある一定の特徴があることが確認されており、それを利用して生成するのが一番いいだろう。とはいえ、今回のエントリは意外に長くなりすぎたので、後日にこれは課題として残しておく。

まとめ

今回は、素数判定テストとして、フェルマーの小定理を使った素数判定法について、自分が理解した範囲のことを書いてみた。これよりも、もう少し正確な方法としては「ラビン・ミラー法」というものが紹介されていたので、こちらも実装してみたい。また、カールマイケル数を出力する方法についても、また後日改めて調べなおしたい。