ρ法

最終更新日:2003.05.31

IDM | 素因数分解 | ρ法

紹介

現在実用されている主要な素因数分解アルゴリズムは、大きく分けて群論系とふるい系に分かれます。しかし、ρ法(Rho method)はその何れにも属さない少々特殊なアルゴリズムです。

名前も特殊で、他のアルゴリズムはその理論的特徴から命名されていますが、ρ法は処理の流れを図にしたときに、ギリシャ文字のρ(ロー)に似た形になることからそのように呼ばれています。

発明者J.M. Pollardの名を冠してPollard's Rho methodと呼ばれたり、乱数を用いた解法であるためモンテカルロ法(Monte-Carlo method)と呼ばれたりもします（註1）。

対象合成数Nの素因数pを発見するため計算量が、Nの大きさよりはpの大きさに依存するという特性があります。そのため、巨大な合成数Nを他のアルゴリズムで素因数分解するための前処理として比較的小さな素因子を取り除いておく目的にも使われます。

一番上に戻る

アルゴリズムの基本

{x_n}を、x_i+1 = f(x_i)によって定まり、x_n mod Nが乱数に近い振る舞いをするような整数列とします。また、{y_i}をy_i = x_2iで定まる数列とします。

i=1,2,3,...に対して順に、g_i = gcd(y_i-x_i, N)を計算します。

このとき、当然の事ながらg_iはNの因数です。それは自明な因数1,Nかもしれませんが、そうでないこともあります。自明でない因数g_iが発見されたとき、Nがg_i×(N/g_i)に分解されたことになります。

これの繰り返しをある程度継続することにより、いくつかの自明でない因数が発見され、結果としてNが分解されていくことが期待できます。

実装に於いては疑似乱数列を生成する関数f(x)として、f(x) = (x²+c) mod N　(ただしcは0,-2以外の整数)を用いるのが普通です（註2）。これは乱数列ではありませんが、かなり乱数に近い振る舞いをすることが経験的に知られています。

次節に述べるように、「mod Nでの振る舞いが乱数に近いこと、また一方ではそれが確定的でもあること」、{x_i}がこの両者を同時に満たすものであることがρ法の要です。

一番上に戻る

説明

上記アルゴリズムはあまりにも不可解に見えることでしょう。このようなまわりくどい手法が意味を持ち得るのは、偏に「Nの大きさには殆ど関係なく、素因数pの大きさにのみ依存した時間でpを発見し得る」という特性に依ります。

基本原理

x_j ≡ x_k (mod p)が成り立つようなj,kがあったとします。

x_j≡x_k (mod p)とは即ちx_j-x_kがpの倍数ということですから、gcd(x_j-x_k, N)を計算することでp（もしくはpを含むNの因数）を得ることができます。

このようなj,kは実際に存在するでしょうか？　x_i mod pについて考えてみます。その値は高々p個の値しかとり得ません。従って（鳩の巣原理より）、x_j₀ ≡ x_k₀ (mod p)を満たすようなp+1以下の値j₀,k₀が存在することになります。

ただし、次項に述べるようにj₀,k₀そのものを求めようとするとかえって無駄が多くなるため、実際にはもう少し大きめのj,kでx_j≡x_k (mod p)なるものを用います。

x_iは完全な乱数列ではなく、x_i+1 = f(x_i)によって定まる数列であるため、 x_j₀ ≡ x_k₀ (mod p)が成立したならば、その後の項についても帰納的にx_j₀+1 ≡ x_k₀+1 (mod p), ...が成立します。即ち、数列{x_i}は遅くとも第p項に至るまでに、mod pで循環を始めます（この循環図の形がアルゴリズムの名前のもとになったものです）。このことにより、条件を満たすj,kは無数に存在しますので、その中で比較的計算しやすいものを探すわけです。

j,kの探し方

さて、x_j ≡x_k (mod p)なる(j,k)を探すといっても、実際にはpは未知ですので、合同式が成立するかどうかを判定することができません。合同式の成立を検出してからGCDを計算するのではなく、とにかく毎回GCDを計算してみる形を取ります。もしGCDが1だったときは合同式が不成立だったのだと考えて次の(j,k)の候補をあたることになります。

ただし、(1,2), (1,3), (2,3), (1,4), (2,4), (3,4), (1,5), ...と、添え字(j,k)のあらゆるパターンの中を探していては効率が悪いです。そこで次の方法でこの問題を解決します。

j₀,k₀を、x_j₀ ≡ x_k₀ (mod p)が成り立つ最小の自然数(ただしj₀<k₀)とし、 dをd≡-j (mod (k₀-j₀))を満たすような最小の自然数とします（k₀-j₀は数列{x_i mod p}の循環周期です）。

このとき、2(j₀+d)-(j₀+d) = j₀+d ≡j₀-j₀ = 0 mod (k₀-j₀) ですから、x_2(j₀+d) ≡ x_j₀+d (mod p)です。即ち、あらゆる(j,k)ではなく(j, 2j)の形をしたパターンのみを探せば、j=j₀+dとなった時点でうまくx_j ≡ x_2j (mod p)となることになります。

d<(k₀-j₀)であり、k₀≦p+1です。これよりx_j ≡ x_2j (mod p)なる最小のjは、j = j₀+d < j₀+(k₀-j₀) = k₀ ≦p+1であり、従ってj≦pです。このことにより、素因数pを発見するまでに必要な計算量がNの大きさよりはpの大きさに依存するというρ法の特性が導かれます。

無保証性

j₀,k₀がx_j₀≡x_k₀ (mod p)だけでなくx_j₀≡x_k₀ mod Nも満たしてしまう場合があります。

言い換えれば、x_i mod Nが循環する周期とx_i mod pが循環する周期が等しい場合です。

このとき、得られるGCDは全てNそのものとなってしまいますので、いつまで経ってもpを発見することはできません。

ρ法による分解では、間違った答えを出すことはありませんが、答えが得られるとは限りません。ただし、実用上はかなりの確率で有効な答えが導かれることが知られています。また、答えが得られなくても、最初に用いたf(x)=x²+cとは異なるcを用いて再試行することで、多くの場合、答えを得ることができます。

一番上に戻る

Rubyによる実装

x_iの値は途中でとても大きくなってしまい、計算コストが高くなります。

どのみち必要なのはNとのGCDであって値そのものではないため、新たなx_iの値を求める毎にNによる剰余だけを保存しておけば、値をN未満に抑えることができます。

GCD(N, x_i-x_j) = GCD(N, (x_i mod N)-(x_j mod N))であるため、x_iの代わりにx_i mod Nを使っても差し支えないわけです。


#!/usr/local/bin/ruby

# Euclidの互除法でGCDを求める関数。
def gcd(a,b)
  while (a%=b)!=0 do
    a,b = b,a
  end
  return b
end

# nのp因子を全て取り除きながら、その都度pを表示する関数
def div_while_divisible(n,p)
  while true do
    quot, rem = n.divmod(p)
    return n if rem!=0
    print "#{p}\n"
    n = quot
  end
  return n
end
  
n = ARGV.shift.to_i    # 第一引数から対象合成数を得る。

c = 1                  # 仮にcは1に固定する
y = x = rand(n)        # x_1をランダムに決定。
y = (x**2+c)%n         # y = (x^2+c) mod n

ubound = Math.sqrt(n).to_i+1      # ループ回数の上限を(√n)+1にする。

i=0
while (i+=1) <= ubound do
  g = gcd(y-x, n).abs     # g = |GCD(y-x, N)|
  if g==n then
    #数列がmod Nで循環を始めたので
    #これ以上続けても新たな素因数は発見できない
  elsif g!=1 then        # gが非自明な因数のとき、
    n = div_while_divisible(n,g)
    ubound = Math.sqrt(n).to_i + 1     # ループ回数の上限を再設定
  end
  
  x = (x**2+c)%n              # x,yの次の項を計算
  y = ((x**2+c)**2+c)%n
end

print "#{n}\n"                # 残った因子を出力

注釈

c

f(x) = (x²+c) mod nにおけるcは、仮に1に固定してしまいました。必要ならば、分解しきれなかったときにcを変えて再試行するような処理を書き加えるべきでしょう。

ubound

上で議論したように、現在の乱数生成関数f(x)によって素因数pを発見可能ならば、それはp+1回ループした時点までに発見されます。

また、素朴試し割り法のときに議論したように、√nまでの素因数が全て発見されれば、その時点で残っている因数もまた素因数です。従って、ループは(√n)+1回行なえば十分であると言えるでしょう。

（他の小さな素因数の発見に失敗して、(√n)以上の大きな素因数p_largeの分解にのみ、(√n)+1回目以降で成功するという可能性もないとは言いませんが、期待薄です）。

一番上に戻る

Brentによる改良

上記では、y_i=x_2iとしておいて、gcd(y_i-x_i, n)をチェックしました。つまり、{x_i}の項のうち、添え字が(i,2i)というパターンを満たすようなのペアに対してGCDをとっていたわけです。

しかし、ここには無駄があります。y₂=x₄を計算する過程で必然的にx₃も求めていますがこの計算結果は捨ててしまっており、gcd(y₃-x₃,n)を計算する段階に於いてx₃を再計算しています。

Richard P. Brentは、(i,2i)というパターンを他のものに変更することを提案しました。(2,1), (3,2), (4,2), (5,4), (6,4), (7,4), (8,4), (9,8), (10,8), ...というものです。一般にiに対しi未満の最大の2のべきが対応するというパターンです。

Brentの改良を取り入れれば、計算した値を無駄なく（しかも過剰でもなく）数列循環の検出に用いることができます。また、yの値は直前のxの値から流用できるため、独自に計算する必要はありません。


#!/usr/local/bin/ruby

# 関数gcd, 関数div_while_divisibleの定義は省略。
# 前節を見よ。

n = ARGV.shift.to_i
c = 1
y = x = rand(n)        # x_1をランダムに決定。 => y_1に代入
x = (x**2+c)%n         # x_2を計算

ubound = Math.sqrt(n).to_i + 2

next_2pow = 2          # 現在のi以上の最小の2の羃
i=1

while (i+=1) <= ubound
  g = gcd(y-x, n).abs
  if g==n then
    break
  elsif g!=1 then
    n = div_while_divisible(n,g)
    ubound = Math.sqrt(n).to_i + 2
  end

  if i==next_2pow then         # iが2の羃に達したら
    y = x                      # y_iを次の項に移す(直前のxの値を使える）
    next_2pow *= 2;              # 次の2の羃を設定
  end

  i += 1
  x = (x**2+c)%n
end

print "#{n}\n"

一番上に戻る

素朴試し割りとの組み合わせ

ρ法の「確実に素因数を発見できるか分からない」という欠点を補うためには、予め素朴試し割り法を用いて、小さな素因数を全て取り除いておくことが有効です。

試し割りの計算量とNの大きさを鑑みて、適当な上限naive_uboundを定め、それ以下の素因数は試し割りで処理することにします。これにより、次のような利点がもたらされます。

2,3,5,...あたりの極めて小さな素因数がNに多数含まれる場合、これらがあちこちのGCDに分散して混入することを防止できます。

GCDは素因子pを返すとは限らず、pを含むNの因数を返す可能性もあります。そのため、極めて小さな素因数sが残っていると、pではなくps, ps², ...などが返される可能性が高くなってしまいます。

p,qがある程度大きな素因数の場合であれば、たまたまpqがGCDとして返ることはあっても、pq²はおそらくNを超えてしまうでしょう。極端に大きな次数を持つGCDが返る心配はないと言えます。しかし、小さな素因数sの場合はあちこちのGCDにs², s³が紛れ込んで収拾がつかなくなり兼ねません。
素数判定が少し容易になります。

得られたGCDがなお合成数である疑いがあるならば、そのGCDを再帰的に素因数分解する必要があります（これまでの実装では行なっていませんが）。しかし、素数であるかの判定にはある程度のコストが掛かります。

naive_uboundまでの素因数が確実に取り除かれているならば、その後のρ法の適用で発見された(naive_ubound)の2乗以下の因数は全て素因数です。ある数Mが合成数ならば必ず√M以下の素因数を持つためです（素朴試し割り法のときの議論を参照してください）。大小比較のみならばコストは低いので、素数判定が少しだけ容易になります。

試し割りを併用した実装


#!/usr/local/bin/ruby

# 関数gcd, 関数div_while_divisibleの定義は省略。
# 前記のソースを見よ。

n = ARGV.shift.to_i

STDERR.print "naive method started\n"

#適当に上限を定めてみる
naive_ubound = (Math.sqrt(n)<0xffff)?Math.sqrt(n):0xffff


p=2; q=3
n = div_while_divisible(n,p)
n = div_while_divisible(n,q)
p=1; q=-1;
while q <= naive_ubound
  p += 6; q += 6;
  n = div_while_divisible(n,p)
  n = div_while_divisible(n,q)
end

naive_ubound_squared = naive_ubound**2;

# 試し割りだけで完全に分解されてしまった場合、そのまま終了
if n <= naive_ubound_squared
  print "#{n}\n" if n!=1
  exit(0)
end


# ρ法の部分は今までと大体同じ

print "rho method started\n"
c = 1
y = x = rand(n)
x = (x**2+c)%n

next_2pow = 2
i=1
rho_ubound = Math.sqrt(n).to_i + 2
while (i+=1)<= rho_ubound
  g = gcd(y-x, n).abs
  if g==n then
    break
  elsif g!=1 then
    n = div_while_divisible(n,g)
    break if n <= naive_ubound_squared    #残りが素数と確定したら、終了
    rho_ubound = Math.sqrt(n).to_i + 2
  end

  if i==next_2pow then
    y = x
    next_2pow *= 2
  end

  i += 1
  x = (x**2+c)%n
end
print "#{n}\n"

一番上に戻る

GCD回数の節約

前節において小さな素因数は取り除いてしまうことにしました。ρ法を適用する時点ではある程度大きな素因数しかないわけですから、合同式y_i≡x_i mod pが成立してGCDが1以外になるのは稀である（周期k₀-j₀が長い）と言えます。

そこで、毎回GCDを計算するのではなく何回かまとめてGCDをとることにします。何らかのrをとっておいて、どうせ大部分のy-xはNと互いに素であろうと考えて、次のようにするわけです。

g = gcd((y_i-x_i)(y_i+1-x_i+1)...(y_i+r-x_i+r), N)

わざわざ積を計算するのは無駄に感じられる知れませんが、積の計算はGCDの計算に比べてかなり速いため、結局は時間短縮になります。

y-xの積はそのまま保存しておく必要はありません。どうせGCDをとるのに使うわけですから、掛けあわせる度にNによる剰余をとってそれを保存しておけば十分です。

rの取り方

rを必要以上に大きくとると、GCDをとる前に複数の素因子が出現してしまうかもしれません。そうすると、それら素因子の積しか得られないわけで、困ってしまいます。適切なrをとることが大切です。

……といっても、具体的にどの程度にすべきかは難しい問題です。以下では、仮にr=10に固定したソースを示します。

GCDを節約した実装


#!/usr/local/bin/ruby

# 関数gcd, 関数div_while_divisibleの定義は省略。
# 前記のソースを見よ。

n = ARGV.shift.to_i

STDERR.print "naive method started\n"

naive_ubound = (Math.sqrt(n)<0xffff)?Math.sqrt(n):0xffff


p=2; q=3
n = div_while_divisible(n,p)
n = div_while_divisible(n,q)
p=1; q=-1;
while q <= naive_ubound
  p += 6; q += 6;
  n = div_while_divisible(n,p)
  n = div_while_divisible(n,q)
end

naive_ubound_squared = naive_ubound**2;

# 試し割りだけで完全に分解されてしまった場合、そのまま終了
if n <= naive_ubound_squared
  print "#{n}\n" if n!=1
  exit(0)
end


print "rho method started\n"
r = 10
c = 1
y = x = rand(n)
x = (x**2+c)%n

next_2pow = 2
i=1
rho_ubound = Math.sqrt(n).to_i + 2
product = 1
while (i+=1)<= rho_ubound
  product *= (y-x)
  product %= n

  if i%r == 0 then
    g = gcd(product, n).abs
    if g==n then
      break
    elsif g!=1 then
      n = div_while_divisible(n,g)
      break if n <= naive_ubound_squared    #残りが素数と確定したら、終了
      rho_ubound = Math.sqrt(n).to_i + 2
    end
    product = 1
  end

  if i==next_2pow then
    y = x
    next_2pow *= 2
  end

  i += 1
  x = (x**2+c)%n
end
print "#{n}\n"

一番上に戻る

実装

実装一覧
言語	所在	サイズ(bytes)	備考
ruby	impl/rho/brent_rho.rb	1,299	文中のものと同一。

一番上に戻る

参考文献

Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest: 精選トピックス　アルゴリズムイントロダクション第3巻 p.195-199, 1995, ISBN4-7659-0247-8
Donald E. Knuth, 準数値算法 The Art of Computer Programing (3), 1986, ISBN4-7819-0426-2
Donald E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms (3rd Edition), 1997, ISBN:0201896842

下記資料は直接参照してはいません。上記の参考文献の中で参考文献として挙げられているオリジナルの論文です。

J.M. Pollard: A Monte Carlo method for factorization. BIT, 15:331-334, 1975
Richard P. Brent: An improved Monte Carlo factorizationi algorithm. BIT, 20(2):176-184, 1980

一番上に戻る

註

モンテカルロ法

「確定的な問題を乱数を用いて解く」類の解法を一般に、賭博都市の名前に因んでモンテカルロ法(Monte Carlo method)と呼ぶそうです。

有名なモンテカルロ法の適用として円周率の計算があります。正方形[-1,1]×[-1,1]内の座標(x,y)をランダムに生成し、それが単位円に含まれるかどうかをチェックします。この試行を繰り返すにつれて、乱座標(x,y)が単位円に含まれた割合は（円の面積）/（正方形の面積） = π/4に近づくと期待できます。このことを利用して、πを求めようというものです。

f(x)

理論上f(x)は、合同関係を保つZ/nZからZ/nZへの写像で乱数的な振る舞いを期待できるものであれば何でも良いです。例えばx³+1でも問題ありません。

ただ、計算コストのことなどを考えるとやはりx²+cあたりが適当なのでしょう。

一番上に戻る

変更履歴

2002.09.05: 全体的に再構成した。
2002.09.04: 公開。

一番上に戻る

IDM | 素因数分解 | ρ法