素因数分解アルゴリズム(特にSQUFOF)のこと

主要な素因数分解アルゴリズム
SQUFOFについて

主要な素因数分解アルゴリズム
- 素因数の性質に依存するアルゴリズム(ρ法、p-1法、楕円曲線法)
- 素因数の性質に依存しないアルゴリズム(SQUFOF、二次ふるい法、一般数体ふるい法)
SQUFOFについて

主要な素因数分解 アルゴリズム

ここでいう素因数分解アルゴリズムは完全に素因数分解をするアルゴリズムではなく因数(約数)をひとつ見つけ出すアルゴリズムなので、完全な分解が必要なら再帰的に実行したり試し割り法と組み合わせる必要がある。

素因数の性質に依存するアルゴリズム(ρ法、p-1法、楕円曲線法)

都合のよい性質を持った素因数を含んでいる場合に成功するアルゴリズム。都合のよい素因数を含んでいる場合は、汎用のアルゴリズムよりも早く素因数を見つけ出したり大きな素因数を取り出したりすることができる。

ρ法
- アルゴリズムは単純。
- 小さい素因数を持っていると(分解したい数そのものは大きくても)早く終了する。
- もっと複雑な素因数分解アルゴリズムを使う場合でも、前処理(小さな素因数を事前に取り除く)や後処理(見つけた因数をさらに素因数分解する)に使われる。
- ρ法の説明→ポラードのρ法の考え方
p-1法
- 素因数 $p$ が「 $p-1$ は小さな素因数しか持っていない」という性質を持っている場合にうまくいく。「小さな素因数」というのがどれくらいの範囲かは補助的に与える数によって決まる(見つけられる素因数の範囲を広くするためには補助的に与える数を大きくしておかないといけない)。
楕円曲線法
- p-1法が整数の剰余乗算(とフェルマーの小定理)を使うのに対し、それに代えて楕円曲線上の点の加法(とラグランジュの定理)を使う。
- 整数係数の楕円曲線をひとつ取り、それをで考えたとき(つまり有限体上で考えとき)の曲線上の点の総数をとする(無限遠点も含めて数える)。このとき「は小さな素因数しか持っていない」という性質がある場合にうまくいく。
  - 楕円曲線 $C$ の取り方によって点の総数 $\sharp C$ が( $p+1-2\sqrt{p} \leq \sharp C \leq p+1+2\sqrt{p}$ の範囲で)変わるので、ある楕円曲線で分解がうまくいかなくても別の楕円曲線を試すことができる。
- 高速化のためのいろいろな手法がある(伊豆哲也、木田祐司「素因数分解と楕円曲線」(『数学のたのしみ 2005春』))。

素因数の性質に依存しないアルゴリズム(SQUFOF、二次ふるい法、数体ふるい法)

SQUFOF (SQUare FOrm Factorization)
- アルゴリズムは(ρ法ほどではないが)単純。
- 分解したい数の桁数が小さい範囲では優秀なアルゴリズム(Jérôme Milan Factoring Small to Medium Size Integers: An Experimental Comparison)。
- ρ法と同様、他のアルゴリズムの後処理などでも使われる。
- GNU Coreutilsのfactorコマンドはこのアルゴリズムをメインに使っている。
二次ふるい法(QS、Quadratic Sieve)
- アルゴリズムは複雑
- 高速。Msieve(複数多項式二次ふるい法MPQSを実装している)では60桁の数(30桁の素数×30桁の素数)の分解で数秒、80桁(40桁×40桁)で数分だった。
一般数体ふるい法(GNFS、General Number Field Sieve)
- 計算量のオーダー評価では現在のところ最速のアルゴリズム。
- 非常に複雑なので大きい数でないと二次ふるい法より速くならない(90桁〜100桁ぐらいから?)。

SQUFOFについて

SQUFOF(Square Form Factorization)について書かれた日本語の文章がほとんどないようなので説明を書いておく(SQUFOFの日本語名も特に存在しないみたい。平方形式分解?)。

連分数をもちいた素因数分解

SQUFOFはM. A. Morrison, J. Brillhart A Method of Factoring and the Factorization of F₇(1975)で使われたアルゴリズム(連分数を利用した素因数分解アルゴリズムの一種)から派生したものなので、まずはそちらを少し説明する。
ρ法は疑似乱数列の中から「 $x\equiv y \bmod p$ 」となるものを(素因数 $p$ を知らずに)見つけるという手法だった。
モリソンとブリルハートは「 $x^2 \equiv y^2 \bmod N$ 」となるもの見つけようとする。もし自明でない $x,y$ がうまく見つかると $\gcd(x\pm y, N)$ によって $N$ の因数が得られる(これは二次ふるい法や数体ふるい法でも使われる考え方)。
それを見つけるために連分数展開を利用する。
因数分解したい数 $N$ について $\sqrt{N}$ を連分数展開していく。

すると次の漸化式が成り立っている。

また連分数展開を途中で打ちきった近似連分数

を考えると、 $A_i,B_i$ も漸化式で表すことができる(この漸化式はSQUFOFでは使わないので省略する。連分数についての多くの解説に書いてある)。
このとき

が成り立つ。 ${}\bmod N$ で考えると

となり、等式の左辺が二乗の形になっている。そこでいくつかの $i$ について $(-1)^i Q_i \bmod N$ を集めてそれをかけ合わせて全ての素因数の指数が偶数になるようにできれば右辺も二乗の形になり「 $x^2 \equiv y^2 \bmod N$ 」となる $x,y$ が得られる。
これがモリソンとブリルハートのアルゴリズムの概要。

シャンクスの改良

モリソンとブリルハートのアルゴリズムは次のようして単純化することができる。

$\sqrt{N}$ の連分数展開を続けていって $i$ が偶数のときで $Q_i$ が完全平方になるものが現れるまで待てば、 $Q_i$ をいくつも集めずに平方となる解を得ることができる。

ただしこの場合、ある程度の長さ展開を続けていかないといけない。連分数の近似分数 $A_i/B_i$ の分子分母は展開が進むにつれどんどん増大していくので、 $A_i$ の増大を抑えるために ${}\bmod N$ で余りを取り続けないといけない。
SQUFOFを考案したシャンクス(D. Shanks)は連分数の展開過程を二次形式の変換に対応させることで、 $A_i$ を全く使うことなく因数を計算する方法を考案した(その理屈は省略する)。

アルゴリズムの説明

(おおよそJ. E. Gower, Samuel S. Wagstaff Jr. Square Form Factorizationの§3.1、3.3、3.5、5.2に基づく)
前提条件: 分解したい数 $N>0$ は、偶数でない、素数でない、完全平方でないこと(また $N$ が小さいと連分数展開の周期が短いので失敗する可能性が大きくなる)。
※ [Gower&Wagstaff]では理論的な取扱いのため $N \equiv 1 \pmod 4$ ときは $N$ を $2N$ に取り替えるという操作をしているが、この操作をしなくてもアルゴリズムは動作し多くの実装でもこの操作をしていないとある(§3.1)。
初期値と漸化式:

前進過程: 連分数展開の漸化式を使って展開を続けていき、が偶数でが完全平方になるものを見つける。
- ただし自明な解になってしまいそうな場合はそのは選ばず展開を続行する。次のように判定する。
  - 見つけた $\sqrt{Q_i}$ がそれ以前の展開のどこかに現れた $Q_i$ ( $Q_i$ が偶数なら $Q_i/2$ )と等しい場合は選択せず展開に戻る(展開に出てくる $Q_i$ はある値より小さいものだけを記憶しておけばいい。後述)。
- 展開のステップ数がある上限 $B$ を越えた場合は失敗とする。また記憶した値の数がある程度増加したら(50個ぐらい?)失敗にする。
後退過程: 前進過程の最後に得た項をに置き換えて、さらに展開を続ける。
- の漸化式は二つ前の値も使っているので、初期値の設定にはの関係を使う。前進過程の最後に得た項をとして次のようになる。
  - $Q'_0 = \sqrt{Q_i}$
  - $P'_1 = P_i + \sqrt{Q_i} \left\lfloor \frac{\large{q_0 - P_i}}{\large{\sqrt{Q_i}}} \right\rfloor$
  - $Q'_1 = \frac{\large{N - {P'_1}^2} }{ \large{\sqrt{Q_i}} }$
- 再び展開をおこない、と次の項が等しくなったところで展開を止める(後退過程での展開のステップ数は前進過程のおよそ半分ぐらい)。　
  - $Q_i$ が偶数なら $\frac{Q_i}{2}$ を見つかった因数とする。
  - $Q_i$ が奇数なら $Q_i$ を見つかった因数とする。

前進過程で $Q_i$ を記憶する部分は次のようにする。

を閾値とする。
- $Q_i$ が偶数ならば、 $\frac{Q_i}{2} \leq L$ の場合に $\frac{Q_i}{2}$ を記録する。
- $Q_i$ が奇数ならば、 $Q_i \leq L$ の場合に $Q_i$ を記録する。
- (小さい $Q_i$ は非常に少ないので、まず $Q_i \leq 2L$ かどうかをチェックして他の判定は後にするのが良い)
完全平方の $Q_i$ が見つかったとき、 $\sqrt{Q_i}$ がそれまでに記憶したどの値とも違っていれば後退過程に進む。
展開ステップの上限回数 $B$ は $B=4L$ 。

例

N=63375401385616362433としてみる。L= $\left\lfloor \sqrt{2\sqrt{N}} \right\rfloor$ =126181

Q₀=1、P₁=7960866874、Q₁=65830557
(記録) i=21127の時、Q_i=223456≦2Lなので、223456/2=111728を記録。
(完全平方) i=42462(偶数)の時、Q_i=12483145984は完全平方。√Q_i=111728は記録されているので展開を続ける。
(記録) i=45890の時、Q_i=114549≦Lなので、114549を記録。
(記録) i=58124の時、Q_i=1803≦Lなので、1803を記録。
(完全平方) i=59398(偶数)の時、Q_i=634082761は完全平方。√Q_i=25181は記録されていないので後退過程に進む。(このときP_i=7554780076)
展開の初期値はQ₀=25181、P₁=7960848882、Q₁=11376179489になる。
i=29807でP_i=P_i+1=7692565499。このときQ_i=7692565499は奇数なので、7692565499を答えとする。

63375401385616362433 / 7692565499 = 8238526067

アルゴリズムの改良

$N$ を展開する代わりに小さな奇数 $m$ を取って $mN$ を展開する: ある $m$ で探索が失敗しても別の $m$ を試すことができる。 $m$ 次第で少ない展開数で $Q_i$ を見つけられるかもしれない。ただし自明な因数 $m$ を除くようにアルゴリズムを少し変更する。

変更点1: $Q_i$ を記録するところ: $\frac{Q_i}{\gcd(Q_i, 2m)} \leq L$ ならば、 $\frac{Q_i}{\gcd(Q_i, 2m)}$ を記録する。
変更点2: $P_i=P_{i+1}$ となったところ: $\frac{Q_i}{\gcd(Q_i, 2m)}$ を見つかった因数とする。
性能の良い $m$ は3,5,7,11およびそれらからいくつかを選んで(全部もあり)かけ合わせたもの。

複数の $m$ で $mN$ を展開していく: 例えば「何百ステップか展開するごとに別の $m$ での展開に移る」というようにして複数の $mN$ を同時展開していくと、どれかの系列で早く $Q_i$ を見つけられてトータルで実行時間を減らせる(かもしれない)。
完全平方な $Q_i$ の選択条件の変更: $Q_i$ の記憶と比較の処理を上で説明したものよりも複雑なものにすると、適切な $Q_i$ の発見までの展開数を減らすことができる(上のやり方だと適切な $Q_i$ を飛ばす可能性がある)。ただし処理時間の点で(どれくらい)有利かは不明。
計算の省力化(アルゴリズム自体の改良ではない): 計算に時間がかかる部分は $q_i$ を求めるときの除算と $Q_i$ が完全平方かを判定するところ。 $q_i$ は $q_i=1$ となることが多い(連分数展開しているので $q_i=0$ になることはない)ので、はじめに $q_i=1$ かどうかを除算なしで判定するというやり方が考えられる。また完全平方かどうかはルートを取ってみないでも判定できる場合がある。

二乗したものを二進数で見ると末尾が必ず $001\overbrace{0\cdots\cdots0}^{2n}$ となる。そのためまずこの形かどうかを調べると判定の効率がよくなるかもしれない。

Schemeによるプログラム例

;; Gauche以外では以下の定義を必要に応じて追加。
;; SRFI-8
(define-syntax receive
  (syntax-rules ()
    ((receive formals expression body ...)
     (call-with-values (lambda () expression)
                       (lambda formals body ...)))))
;; R6RS
(define (exact-integer-sqrt x)
  (let ((sqrt-x (inexact->exact (floor (sqrt x)))))
    (values sqrt-x (- x (* sqrt-x sqrt-x)))))

;; Gauche固有
(define (floor->exact x)
  (inexact->exact (floor x)))

;; Gauche固有
(define-syntax let1
  (syntax-rules ()
    ((_ var expr expr0 ...)
     (let ((var expr)) expr0 ...))))

;; 入力値のチェックや記録したQiの個数が多くなったときの処理は省略した
;; n>0 nは偶数、素数、完全平方ではないこと
;; 乗数multは小さな奇数 1, 3, 5, 7, 11, 3*5, 3*7, etc.
(define (squfof n mult)
  (let* ((nm (* n mult))
         (L (floor->exact (sqrt (* 2 (sqrt nm)))))
         (m*2 (* 2 mult))
         (L*m*2 (* m*2 L))
         (max-steps (* 4 L)))
    (receive (q0 nm-q0*q0) (exact-integer-sqrt nm)
      ;; 以上で初期パラメータ計算終了

      ;; PiとQiの漸化式
      (define (next-terms Qi-prev Pi Qi)
        (let* ((qi (quotient (+ q0 Pi) Qi))
               (Pi-next (- (* qi Qi) Pi))
               (Qi-next (+ Qi-prev (* qi (- Pi Pi-next)))))
          (values Pi-next Qi-next)))

      ;; Qi/gcd(Qi, 2m) <= Lの時、Qi/gcd(Qi, 2m)を記録する
      (define (update-saved-Qis saved-Qis Qi)
        (if (not (<= Qi L*m*2))
          saved-Qis
          (let1 tmp (/ Qi (gcd Qi m*2))
            (if (<= tmp L)
              (cons tmp saved-Qis)
              saved-Qis))))

      ;; 前進過程
      (define (forward-cycle i Qi-prev Pi Qi saved-Qis)
        (if (> i max-steps)
          #f
          (let1 saved-Qis (update-saved-Qis saved-Qis Qi)
            (receive (Pi-next Qi-next) (next-terms Qi-prev Pi Qi)
              (if (even? i)  ; Q_{i+1}を調べるので偶奇判定が逆になる
                (forward-cycle (+ i 1) Qi Pi-next Qi-next saved-Qis)
                (receive (sqrtQ rem) (exact-integer-sqrt Qi-next)
                  (if (or (not (zero? rem))
                          (memv sqrtQ saved-Qis))
                    (forward-cycle
                      (+ i 1) Qi Pi-next Qi-next saved-Qis)
                    (backward-cycle Pi-next sqrtQ))))))))

      ;; 後退過程
      (define (backward-cycle P-last sqrtQ)
        (let1 Pi-tmp (+ P-last
                        (* sqrtQ (floor->exact (/ (- q0 P-last) sqrtQ))))
          (let loop ((Qi-prev sqrtQ)
                     (Pi Pi-tmp)
                     (Qi (/ (- nm (* Pi-tmp Pi-tmp)) sqrtQ)))
            (receive (Pi-next Qi-next) (next-terms Qi-prev Pi Qi)
              (if (not (= Pi Pi-next))
                (loop Qi Pi-next Qi-next)
                (/ Qi (gcd Qi m*2)))))))

      (forward-cycle 1 1 q0 nm-q0*q0 '()))))