二次篩法 (QS) (1)

改めて GNFS の勉強をし直してみたくなったので脳内で分かってるつもりになっている事を確認しつつ実装していく予定で書き始めてみました。まずは理論も実装も簡単な QS からスタートです。

篩法

概要

まずは今後やっていく篩系の素因数分解アルゴリズムに共通する概要を説明します。各篩法では最終的に因数分解したい合成数 $N$ に対して

$$x^2 \equiv y^2 \pmod{N}$$

となる非自明な整数組 $(x,y)$ を見つけることになります。
すると $x^2-y^2=(x+y)(x-y)\equiv 0 \bmod{N}$ から $x\pm y$ が $N$ の素因数をバラバラに持っていることを期待して $x\pm y$ と $N$ の最大公約数を求めます。その結果が $1$ や $N$ になった場合は残念ながら分解失敗となりますが、他の値が出ればそれは $N$ の(素)因数ということで成功となります。

ではこの $(x,y)$ の組み合わせはどのように求めるのか、という点がそれぞれのアルゴリズムで異なります。$N=667$ を例として以下の 5 式を見てみましょう。

$$\begin{align*} 31^2 & \equiv 294 = 2 \cdot 3 \cdot 7^2\\ 32^2 & \equiv 357 = 3 \cdot 7 \cdot 17\\ 37^2 & \equiv 35 = 5 \cdot 7\\ 38^2 & \equiv 110 = 2 \cdot 5 \cdot 11\\ 39^2 & \equiv 187 = 11 \cdot 17\\ \end{align*}$$

この両辺をかけると

$$\begin{align*} (31\cdot32\cdot37\cdot38\cdot39)^2 &\equiv (2\cdot3\cdot5\cdot7^2\cdot11\cdot17)^2\\ 144^2 &\equiv 86^2 \end{align*}$$

となるので $(x,y)=(144,86)$ として ${\rm GCD}(N,x+y)=23$ と分解できます。

この一連の流れが「篩法」に分類される素因数分解アルゴリズムの概要で、それぞれのアルゴリズムでは主に最初にあげたような関係式を見つける方法が異なります。

二次篩法

ここではこの分野で最も基礎的な二次篩法 (QS; Quadratic Sieve) を紹介します。歴史的にはこれ以前にも連分数法という同様の手法がありましたが、篩らしい篩法ではないので割愛します。

概要

二次篩法は以下のように大きく4つのステップから構成されます。

1. パラメータの選択
2. 篩
3. 有効式の選択
4. 分解

最初は後の作業で探索する範囲を決めるパラメータを設定します。具体的には $x$ の探索範囲 $[a,b)$ と使う素数の上限 $P$ です。なお $[a,b)$ については方法としては必須ではないのですが計算時間を短くするという目的から $\sqrt{N}$ を含む範囲で設定することを仮定します。
また、探索幅 $b-a$ や $P$ は小さければ探索時間が短くなりますが、必要なデータが十分に揃う確率が低くなります。逆に大きければ必要データは揃う確率は上がりますが探索に時間を要するようになります。
なるべく短い時間でほぼ確実に分解するためには、少し小さめの桁数の分解で実験を繰り返した経験からこれらのパラメータに最適と思われる値を設定することになります。

ステップ 2. の篩作業では 1. で決めた範囲 $[a, b)$ にある整数 $x_k$ について

$$f(x) = x^2 - N$$

という式の値の素因数分解を考えます。つまり $x^2 \bmod N$ 一般ではなく $f(x)$ に限定します。すると探索範囲として $\sqrt{N}$ を含むあまり広くない範囲を設定しているはずなので $f(x_k)$ は $\sqrt{N}$ 程度の大きさになりそこそこ小さな素数だけで分解できることが期待されます。

例として $N=187$ について $[a,b)=[4,16)$、$P=20$ と設定して計算した結果を見てみましょう。$P$ 以下の素数だけで分解できる ($P$-smooth という) 値だけ抜き出しました。

$$\begin{align*} 4^2 & \equiv -171 = (-1) \cdot 3^2 \cdot 19\\ 5^2 & \equiv -162 = (-1) \cdot 2 \cdot 3^4\\ 11^2 & \equiv -66 = (-1) \cdot 2 \cdot 3 \cdot 11\\ 13^2 & \equiv -18 = (-1) \cdot 2 \cdot 3^2\\ 14^2 & \equiv 9 = 3^2\\ 15^2 & \equiv 38 = 2 \cdot 19 \end{align*}$$

次にステップ 3. の有効式の選択では前のステップで得られた $f(x_k)$ の情報から、かけ合わせて平方数になる組み合わせを選択します。上記の例では $\bm{x}=\{4,5,15\}, \{5,13\}, \{5,13,14\}$ といった組み合わせが考えられます。

このステップは見方を少し変えると各素数の指数を $\bmod{2}$ で考えたときに従属関係になる式の組み合わせを求める操作になりますので、一般的にはそのビット列を元に計算する方式が取られています。

そして 4. の分解ステップで既に篩法の概要で述べたように ${\rm GCD}(N, x\pm y)$ の計算をすることで非自明な因数を見つけられることが期待されます。ステップ 3 で複数の従属組み合わせが見つかっていれば、そのいずれかから非自明な解を得られる可能性が高くなります。

アルゴリズム

ここまでは数式の上で原理を理解するための概要説明であり、具体的な手順として不明なところや素直に実行しては効率的ではないところがあったかと思いますので順に説明しましょう。

篩

篩ステップでは $f(x_k)$ の素因数分解を行うことになります。が、これを直接計算すると結局試し割りの連続となるので効率がよくありません。

for (Integer x = a; x < b; ++x) {
  Integer fx = f(x);
  for (int p : factor_base) {
    while (fx % p == 0) {
      exponent[x][p]++;
      fx /= p;
    }
  }
  cached_f[x] = fx;  // == 1 なら FB-smooth
}

なので複数の $f(x_k)$ の因数分解を纏めて効率的に行う事を考えます。
素数 $p$ がある $x=x_k$ について $f(x)$ を割り切れたとしましょう。

$$p | x_k^2-N$$

すると $x=x_k+p$ についても $p$ で割り切ることができるはずです。

$$f(x_k+p) = x_k^2 + 2px_k + p^2 - N = f(x_k) + p(2x_k+p) \\ \therefore p | f(x_k+p)$$

なので $f(x)$ を $p$ で割り切れる $x$ を 1 つ見つける事ができればその前後 $p$ 毎に $p$ で割り切れることがわかります。そしてそのような $x$ は方程式

$$x^2-N \equiv 0 \pmod{p}$$

を満たすので、$\bmod{p}$ における $N$ の平方根を探ることになります。詳細は別記事に纏めましたのでご覧ください。

for (Integer x : Range(a, b)) {
  cached_f[x] = f(x);
}
for (int p : factor_base) {
  for (int r : ModSqrts(n % p, p)) {
    for (int pk = p;; pk *= p) {
      Integer a0 = GetStart(pk, p, n % p, r);
      if (a0 >= b) break;
      for (Integer x = a0; x < b; x += p) {
        exponent[x][p]++;
        cached_f[x] /= p;
      }
    }
  }
}

ModSqrts(a, p) は $\bmod{p}$ での $a$ の平方根2つを返す関数で、GetStart()は $r^2\equiv n \pmod{p^k}$ を満たす $r$ のうち、探索範囲内の最小のものを求める関数です。
さて、この方法の最内ループを見てみるとエラトステネスの篩と同じような作業をしていることが分かるでしょう。これが「篩系」アルゴリズムと言われる理由です。

長くなりそうなので　(2) に続く