素因数分解プログラムについてのメモ

http://blog.practical-scheme.net/shiro/20130216a-factorizeを読んで疑問に思ったので調べた。

敢えて実装済みのMonte Carlo Factorizationで挑戦してみよう。

gosh> (use math.prime)
#<undef>
gosh> (time (mc-factorize 280671392065546467397265294532969672241810318954163887187279320454220348884327))
;(time (mc-factorize 280671392065546467397265294532969672241810318954163 ...
; real 552.461
; user 552.310
; sys    0.110
(162425297 215940091 358456949 369941863 369941863 479871607 706170617 481362815814826159)

9分ちょいかかった。

試しにGauche 0.9.4 pre3を手元で動かしてみるともっと時間がかかった。

gosh> (time (mc-factorize 280671392065546467397265294532969672241810318954163887187279320454220348884327))
;(time (mc-factorize 280671392065546467397265294532969672241810318954163 ...
; real 884.309
; user 883.320
; sys    0.050

で本題は次の部分。

GNFSでこれが1秒で計算できるなら、 Monte Carlo FactorizationをスキップしてGNFS実装した方が良かったかなあ。

GNFS(一般数体ふるい法)は少し解説を読んだだけでもかなり実装がたいへんそうなアルゴリズムなので、実際にGNFSの実装をうたっているmsieveはともかくPerlやPythonのそれほど専門向けではなさそうなライブラリでGNFSを使っているのだろうかと疑問に感じた。それで調べてみるとやはりGNFSを使っているわけではないようだ。
どうやらPythonのSymPyはポラードのρ法とp-1法を使っていて、PerlのMath::Prime::Utilはそれに加えて(特にGNU Multi-Precision Library使用版で)いくつかのアルゴリズムを併用しているみたい。

他の実装との速度比較

上と同じ数について因数分解の時間を計測してみると次のようになった。

Python  factorint                           1.70(sec.)
Python  factorint (only rho method)         4.91
Python  factorint (only p-1 method)        13.32
Perl    factor (GMP)                        0.19
Perl    factor (Pure Perl)                236.28

こうしてみるとアルゴリズム上あまり違いがないと思われるPythonのρ法版と比べてGaucheの結果が非常に悪い。しかしGaucheがそんなに遅いとは思えない。(ABC予想についてのデータを作ったときも十分速かった印象がある)
ということはアルゴリズムの実装のどこかにバグがありそう。

ポラードのρ法

Gaucheはρ法(モンテカルロ法)のブレントによる変形版を実装しているのだけど、ポラートによるオリジナルのアルゴリズムがかなり単純なのでそっちをWikipediaを見ながらコードに移した。

;;; Gaucheの実装でmc-find-divisor-1と呼ばれている部分を置き換え。
(define (mc-find-divisor-1 n x0)
   (pollard-rho n x0))

(define (pollard-rho n x0)
  (define (f x) (remainder (+ (* x x) 2) n))

  (let loop ((U x0)
             (V x0))
    (let ((U (f U))
          (V (f (f V))))
      (let1 g (gcd (- U V) n)
        (cond ((= g 1) (loop U V))
              ((= g n) #f)
              (else g))))))

実行してみる。

gosh> (time (mc-factorize 280671392065546467397265294532969672241810318954163887187279320454220348884327))
;(time (mc-factorize 280671392065546467397265294532969672241810318954163 ...
; real   3.183
; user   3.160
; sys    0.000
(162425297 215940091 358456949 369941863 369941863 479871607 706170617 481362815814826159)

いっきに速くなってPythonと比べても遜色のない速度になった。やはり元の実装のどこかに問題があったみたい。ただどこに問題があるのかよく分からない。いかにも手続き的な感じの多重ループは(代入を避けると)Schemeのコードに移しづらいので、元のアルゴリズムとSchemeのコードをうまく比較できなかった。

ブレントのアルゴリズム

しょうがないのでブレントのアルゴリズムを半ば機械的にSchemeのコードに移して、同じ因数分解を実行してみた。

(time (mc-factorize 280671392065546467397265294532969672241810318954163887187279320454220348884327))
;(time (mc-factorize 280671392065546467397265294532969672241810318954163 ...
; real   0.729
; user   0.730
; sys    0.000
(162425297 215940091 358456949 369941863 369941863 479871607 706170617 481362815814826159)

さらに速くなった。何か重要な処理を飛ばしているために速くなっている可能性もある(変形版の方が遅い?みたいな文章もあった)。それにρ法の実行速度は入力した数の大きさではなく素因数の大きさで決まるようなのでこれだけ見ても性能はよくわからない。素因数の大きさに依存しないアルゴリズムを使っても80桁ぐらいの数で数時間くらいかかるみたい(Appendix 1. Factorization resultsなど)。
けれど、取りあえずよいことにする。終わり。
(追記: 『UBASICによるコンピュータ整数論』(電子版)P.163の

問題 13 rho1 ではGCDを毎回計算している。これは大変時間のかかる計算なので普段はA1-A2[上のプログラムでの(- U V)]の積を作っておいて100回に1回だけGCDを計算するようにプログラムを改造せよ。

に従ってポラードのρ法のプログラムを書き換えたらいくらか速くなった。

(define (mc-find-divisor-1 n x0)
  (pollard-rho2 n x0))

(define (pollard-rho2 n x0)
  (define (f x) (remainder (+ (* x x) 2) n))
  (define N-loop 100)

  (let loop ((U x0)
             (V x0))
    (receive (U V W) (let inner-loop ((i 0) (U U) (V V) (W 1))
                       (if (> i N-loop)
                         (values U V W)
                         (let ((U (f U))
                               (V (f (f V))))
                           (inner-loop (+ i 1) U V
                                       (* W (remainder (- U V) n))))))
      (let1 g (gcd W n)
        (cond ((= g 1) (loop U V))
              ((= g n) #f)
              (else g))))))

追記終わり)

追記: ポラードのρ法の考え方

こういう場合、彼は思いがけぬ偶然にたよったのである。合理的な連脈の糸をたどれぬこのような場合には、彼は冷静に、そして慎重に不合理の連脈に従った。(…)なにか一つでも手がかりがある場合なら、こんな方法は下の下だが、ぜんぜん手がかりがない場合には、これこそ最善の方策だ。
(G・K・チェスタトン「青い十字架」(中村保男訳))

合成数 $n$ の因数(素因数とは限らない)を見つけたいとする。
ポラードのρ法は簡単に言えば次のようなアルゴリズム。

ランダムな数の列 $x_1,x_2,\ldots,x_k,\ldots$ を作っていく。これは単純な擬似ランダム列でかまわないけれど、線形合同法は望ましくない。

例えば適当に選んだ初期値 $x_0$ と漸化式 $x_{k+1}=f(x_k)=(x_{k}^2+2) \bmod n$ を使う。

$\gcd(x_{k}-x_{2k},n)>1$ となるものを見つける。次のように計算する。

$x_{k}$ と $x_{2k}$ を記録するための変数 $U,V$ を用意し、初期値を $U=V=x_0$ とする。

$U\leftarrow f(U)$ と $V\leftarrow f(f(V))$ を計算し $\gcd(U-V,n)$ を求める。gcdの値が1以外になるまで $U,V$ を更新し続ける。

見つかった $\gcd(x_{k}-x_{2k},n)>1$ は $n$ の因数になる。

ただし見つけたものがnだった場合はアルゴリズム失敗。ランダム列の初期値を取りかえてやり直すか、あきらめる)。

このアルゴリズムは、次の性質に基づいている。

nの素因数のひとつをpとする。 $0 \leq a,b \leq n$ となる $a$ と $b$ について

$a \equiv b \bmod p$ ならば、 $g =\gcd(a-b,n)>1$ となる。このとき $g$ は $n$ の因数である(ただし $a$ と $b$ が等しい場合は $g=n$ となる)。

例えばn=299として、ランダムにaとbを選んで、a=208,b=52を選んだとする。
このとき偶然にも、nのある素因数pについて $a \equiv b \bmod p$ が成り立っている。そのため
$\gcd(208-52,299)=13\,>1$
となり、299の因数(この場合、素因数)13が得られる。
このような偶然がどれくらいの確率で起こるかを考えてみる。
$a$ と $b$ をデタラメに選んでみる。このとき $a \equiv b \bmod p$ となる確率は $1/p$ 。したがって $1/p$ とだいたい同じ程度の確率で $g =\gcd(a-b,n)>1$ となる。(ここで $p$ は最小の素因数)。
そして $g>1$ のとき因数 $g$ が見つかったことになる(ただしその見つかった因数が $g=n$ である可能性もある。確率は $p/n$ ぐらい)。

しかし $a$ と $b$ を毎回ランダムに選ぶと $a \equiv b \bmod p$ となるまでに平均 $p$ 回程度の繰り返しが必要になってしまい意味がない(それなら普通に小さい順に $n$ を試し割りしていけばいい)。
もっと少ない回数で $a \equiv b \bmod p$ となるものを見つけたい。
そのために次のようにする。

$\{0,\ldots,n-1\}$ の値を取り乱数っぽく振る舞う関数 $f(x)$ を用意する。この関数は $x\equiv y \bmod p$ ならば $f(x)\equiv f(y) \bmod p$ となる性質を持っている必要がある(この性質は $f(x)$ が「多項式(mod n)」の形のときは満たされている)。そして適当な初期値 $x_0$ から始めて

という疑似乱数列を作っていく。この系列の中から $x_{k} \equiv x_{l} \bmod p$ となるものを見つける。 $\bmod p$ の余りは $p$ 通りしかないので、いずれどこかで(系列の $r$ 番目と $s$ 番目)で

となるものが出てくる。このとき $f(x)$ の性質から、ここから先は

というように常に一致する。つまり ${}\bmod p$ で考えると、 $s$ 番目から系列の循環が始まる(ロー法という名前はこの循環にいたる様子を「ρ」の形に見立てている)。
そして系列が循環していることはいわゆるウサギとカメのアルゴリズム(フロイドの循環検出法 - Wikipedia)で検出できる。
ふたつの系列を計算していき

$g =\gcd(x_{2k}-x_{k},n)>1$ となるかを見ていれば循環が検出でき $n$ の因数が見つかることになる。

ρ法にかかる時間の見積り

ここで問題になるのは同じ値が登場するまでにどれくらいかかるか。最悪の場合( ${}\bmod p$ でのすべての値を網羅する場合) $p$ 回かかってしまう。
もし値をランダムに選んでいったとするといわゆる「誕生日問題」になり、同じ値が出るまでの平均回数はおおよそ $\sqrt{p}$ に比例した値になる。そのため $f(x)$ 系列をランダムだとみなせば循環までにかかる回数は $\sqrt{p}$ に比例すると見積もれる(ただし $f(x)$ が線形合同法だと循環の長さが $p$ に比例してしまうのでρ法と相性が悪い)。
実行にかかる時間は $n$ よりも素因数 $p$ の大きさに支配されている。もちろん $n$ で割った余りや $n$ との最大公約数を取る処理があるので $n$ にも依存しているけれど。そのため比較的小さな素因数を持っている場合は高速に実行できると予想される。うまくいかないのは次のとき。

素因数がどれも大きいために循環するまでにかかる時間が大きくなる場合(→他のアルゴリズムを使うかあきらめる)。
循環までにかかる時間が期待値よりもずっと大きくなる場合(単に運が悪かった? $f(x)$ や初期値 $x_0$ と $n$ の相性が悪かった?)。
見つかった因数 $g$ が元の数 $n$ と等しかった場合( $f(x)$ がランダムなら非常に低確率。 $n$ と $f(x)$ の相性が悪いと発生しやすいかもしれない?)。

ブレントのアルゴリズムのコード

使ったブレントのアルゴリズムのコードは以下のもの。

(define (mc-find-divisor-1 n x0)
  (brent n x0))

(define (brent n x0)
  (define m (ceiling->exact (log n 2)))
  (define (f x) (remainder (+ (* x x) 2) n))
  (define (f^ k x)
    (let loop ((k k) (x x))
      (if (= k 1)
        (f x)
        (loop (- k 1) (f x)))))

  (define (first-step)
    (outer-loop))

  (define (outer-loop)
    (let loop ((y x0) (r 1) (q 1))
      (let* ((x y)
             (y (f^ r y)))
        (receive (g y ys q) (middle-loop x y r q)
          (let1 r (* r 2)
            (if (> g 1)
              (second-step g x ys)
              (loop y r q)))))))

  (define (middle-loop x y r q)
    (let loop ((k 0) (y y))
      (let1 ys y
        (receive (y q) (inner-loop (min m (- r k)) x y q)
          (let ((g (gcd q n))
                (k (+ k m)))
            (if (or (>= k r) (> g 1))
              (values g y ys q)
              (loop k y)))))))

  (define (inner-loop ntimes x y q)
    (let loop ((i 1) (y y) (q q))
      (if (> i ntimes)
        (values y q)
        (let* ((y (f y))
               (q (remainder (* q (abs (- x y))) n))) ; absは不要のはず
          (loop (+ i 1) y q)))))

  (define (second-step g x ys)
    (cond ((not (= g n)) g)
          (let loop ((ys (f ys)))
            (let1 g (gcd (abs (- x ys)) n)  ; absは不要のはず
              (if (> g 1)
                (if (= g n) #f g)
                (loop (f ys)))))))

  (first-step))