フロイドの循環検出法とポラード・ロー法による素数判定

今日の料理

自宅超会議の様子です。クローズドソース。あと汚ない。

概要

ある数の因数を見つけ出すためには、ポラード・ロー法という方法がある。このポラード・ロー法は、フロイドの循環検出法と呼ばれる、循環を見つける方法が関連してくる。今回は、簡単なリストを使い、この循環法を確認したのちに、ポラード・ロー法について説明をする。

はじめに

最近の日課は素数について調べることになっているのだけれど、yuguiさんのサイトにて『UBASICによるコンピュータ整数論』という本が紹介されていて、「へー」なんで思っていたりした。その中に、ロー法と呼ばれる奴が紹介されていたので実装してみるといいのかなと思って実装することにしてみた。以下はそのメモである。

循環する合同式

という式における任意の数pについて観察してみよう。例えば、53の変化の一覧を出したい場合、次のようにコードが書ける。

def modulo_test p, i
  a = 2
  result = [a % p]
  1.upto i do |i|
    a = (a ** 2) % p
    result << a % p
  end
  return result
end

このとき、modulo_test(53, 50)としたときの一覧を出力すると、下のようになる。

[2,  4, 16, 44, 28, 42, 15, 13, 10, 47, 36, 24,
46, 49, 16, 44, 28, 42, 15, 13, 10, 47, 36, 24,
46, 49, 16, 44, 28, 42, 15, 13, 10, 47, 36, 24,
46, 49, 16, 44, 28, 42, 15, 13, 10, 47, 36, 24,
46, 49, 16]

見ればわかるように、3番目から始まって、12番目で循環している。このように、任意の数を割った余りはどこかで循環するように出来ている。この数同士を周期的のペアと頭においておく。

そこで、同様に53の2乗である2809で変化を観察してみよう。

[  2,   4,   16,  256,  929,  678, 1817,  914, 1123, 2697, 1308,  183,
2590, 208, 1129, 2164,  293, 1579, 1658, 1762,  699, 2644, 1944, 1031,
1159, 579,  970, 2694, 1989, 1049, 2082,  437, 2766, 1849,  248, 2515,
2166, 526, 1394, 2217, 2148, 1526,   15,  225,   63, 1160,   89, 2303,
417, 2540, 2136]

ぱっと見、循環しているように見えないが、例えば1031 - 183 = 848 と2809の最大公約数を求めると、53という数が出てくる。また、同様に2515 - 1031 = 1484 と2809の最大公約数のペアにしても、同じように53が求められる。

理屈的には、任意の場所でとなるようなiとjが確率的に期待される場合、この差はpの倍数となる。それはそうで、ある数から余りを引けば、元の割った数の倍数になるからだ。

これを応用すれば、pの因数を見つけだしたいnに対し、というアプローチが考えられるわけだ。

とはいえ、理屈はそうだとしても、問題はこのようなペアを探すことだろう。ここで、作った配列の要素一つごとに、その差と調べたい数の最大公約数を見つけるのは、試し割り法、つまり数を一つずつ割って試すのと変わらない。むしろ直感的なだけ、後者のほうがマシである。

問題は循環していて、差が倍数になるような数を見つけ出す方法だ。このときフロイドの循環検出法が出てくる。

循環する連結リストの作成

ここで、先ほどの問題から離れて、単純になんらかの循環したデータについて、ループしているかどうかを判定するための方法について考察することにする。

そこで、循環したデータ構造をてっとりばやく作るなら、連結リストだろうということになる。実践的な連結リストはRakeの中にあるのだが、もうすこしシンプルにしよう。

ここで考える連結リストのかたちは、自分の値と、次への要素にリンクする単純なものだ。アスキーアートにするならば、次のような感じだ。

[ 1 |  ] --> [ 2 |  ] --> NULL

そこで、次のようなクラスを作成する。

class SimpleListNode
  def initialize(value)
    @value = value
    @next = nil
  end

  def cons(nextnode)
    @next = nextnode
    return self
  end

  def empty?
    @next.nil?
  end

  def next
    @next
  end

  def to_s
    "[value=#{@value} | ] --> " +
      if empty?
        "NULL"
      else
        @next.to_s
      end
  end
end

必要最小限に整えたものである。これを使う場合:

[1] pry(main)> require './linkedlist.rb'
[2] pry(main)> a = SimpleListNode.new "hoge"
[3] pry(main)> b = SimpleListNode.new "fuda"
[4] pry(main)> a = a.cons b
[5] pry(main)> a.to_s
=> "[value=hoge | ] --> [value=fuda | ] --> NULL"

といったように使うことができる。

さて、さっそく循環リストを作ってみる。

[6] pry(main)> b = b.cons a
[7] pry(main)> a.to_s
SystemStackError: stack level too deep
from /usr/lib/ruby/vendor_ruby/pry/pry_instance.rb:355

スタックレベルが深くなりすぎているということは、相互参照によって呼びだしが際限なく続いているからで、ひとまず成功だと言える。

フロイドの循環検出法の考え方

ここで、どのようにループを検出するかを考えるわけだ。単純に考えると、要素を配列に次々と入れていって、配列の中に存在している要素とでくわした場合、ループと判定したらいいと考えられる。これは力まかせの方法だが、もうすこし良い方法がある。

ところで、良くギリシアの古典的な問題に、アキレスと亀が使われる。いわば、アキレスが前を歩いている亀に追いつくためには、亀が通った点を通らなければ……、という奴である。

そこで、位置関係を逆にして、アキレスに先に走ってもらい、亀は後ろを走ってもらうことにしよう。亀よりもアキレスのほうが速いので、アキレスが前を走っているならば、亀は追いつくことができない。なので、もし走るコースが一直線になっているならば、アキレスはゴールまで一直線になる。

しかし、ここで、もしコースがどこかで周になっていたらどうだろうか。そうするとアキレスはいつのまにか亀を後ろから追い抜いてしまうだろう。

これがフロイドの循環検出法の考え方だ。これを、循環 した連結リストでチェックする場合、次のようになる。ついでにto_sのほうも、ループかどうかを調べられるようにしよう。

class SimpleLinkedNode

  # ...

  def loop?
    if empty? || @next.empty?
      return false
    end
    slow_p = @next
    fast_p = @next.next
    while slow_p || fast_p
      return true  if slow_p == fast_p
      return false if slow_p.empty? || fast_p.empty? || fast_p.next.empty?
      slow_p = slow_p.next
      fast_p = fast_p.next.next
    end
  end

  def to_s
    "[value=#{@value} | ] --> " +
      if loop?
        " .. INFINITY .."
      elsif empty?
        "NULL"
      else
        @next.to_s
      end
  end

  # ....

end

あわせて読みたい

• 循環小数の循環節を求める「ウサギとカメのアルゴリズム」 : 未知の破片

ポラード・ロー法の実装

このような考え方が、先のポラード・ロー法にも使えるということになる。

def rho n
  a = 1 + Random.rand(n - 3)
  x = y = Random.rand(n)
  d = 1
  g = ->(x) { return (x ** 2 + a) % n }
  while d == 1
    x = g[x]
    y = g[g[y]]
    d = (x - y).gcd(n)
  end
  return d
end

このコードは、ここを参考に若干の手直しを行なっている。

ここで唐突に乱数を使っているので、あれっとなるかもしれないけれど、実はこのポラード・ロー法、乱択アルゴリズムなので、一回だけでやると、合成数もこのテストをすりぬけてしまう。実際、一回だけやると次のような素数リストになってしまう。

[4, 5, 7, 8, 9, 10, 11, 13, 15, 16, 17, 19, 20,
21, 23, 27, 29, 31, 32, 33, 37, 41, 43, 47, 49,
53, 54, 56, 59, 61, 67, 71, 73, 79, 80, 83, 85,
89, 93, 95, 97]

素人目に見ても、明かに合成数だな、とわかるほどの乱雑さ。しかし、試行回数を重ねると、明かにふるい落とされているのがわかる。以下は100000のリストに対して、10回ほど重ねることで、判定をミスった数をふるい落としている様子だ。要素を載せると大変なので、素数候補を集めた配列のサイズを載せておく。

9989
9607
9597
9592
9592
9592
9592
9592
9592
9592

この結果からわかるように、4回目から結果が安定してきている。

ちなみに、Rubyには組み込みライブラリとしてPrimeがあり、それを利用したところ、9592という結果になったので、ほぼ確実に素数だけ残ったと言える。

おわりに

とはいえ、実際のポラード・ロー法の利用法としては、任意の大きな数の素因数を求めるのに使われるので、素数判定として使うのは若干変なのだけど、概要だけつまむのには丁度いいかな、と思ったので、この形式になった。

また、今回は省略したけれども、ポラード・ロー法のフロイド循環検出も余計であるとして、これを使わない方法もある。

参考文献

• lecture notes - Factoring