🤗

二次篩法 (QS) (2)

2022/05/15に公開

math

メモ

tech

二次篩法(QS) (1)の続きです。

選択

まずは「使える」関係式 (factor base で分解できた $f(x)$ ) を洗い出しましょう。

struct Smooth {
  Integer x;
  vector<int> exponents;
};

// Smooth な関係式だけ取り出す。
vector<Smooth> smooths;
for (Integer x : Range(a, b)) {
  if (cached_f[i] > 1) continue;
  smooths.push_back(Smooth {
    x,
    exponent[x]
  });
}
// これ以降 exponent と cached_f は捨てて大丈夫。

そしてその中からどれを使って分解に使う $(x,y)$ を構成するか見つけます。
まずは概要で述べたように指数の $\bmod{2}$ で考えるので変換します。つまり指数の

x	-1	2	3	11	19
4	1	0	2	0	1
5	1	1	4	0	0
11	1	1	1	1	0
13	1	1	2	0	0
14	0	0	2	0	0
15	0	1	0	0	1

が $\bmod{2}$ で

x	-1	2	3	11	19
4	1	0	0	0	1
5	1	1	0	0	0
11	1	1	1	1	0
13	1	1	0	0	0
14	0	0	0	0	0
15	0	1	0	0	1

となります。これにどの数式を使うか、という情報を付加させたいので $k$ 番目の式に $k$ 番目のビットを割り当てます。そしてこのビット情報があればどの式、つまりどの $x$ を利用するか分かるので $x$ の情報も削ります。

-1	2	3	11	19	式番号bit
1	0	0	0	1	000001
1	1	0	0	0	000010
1	1	1	1	0	000100
1	1	0	0	0	001000
0	0	0	0	0	010000
0	1	0	0	1	100000

また、ビット情報は疎であることが予想されるため、指数と式番号のビット情報は番号の set で持ちます。

struct Relation {
  set<int> indexs;  // indecies にするとパターン化がメンドイ
  set<int> exponents;
};

// 式情報をビット圧縮情報に変換
for (int i = 0; i < smooths.size(); ++i) {
  Relation r;
  r.indexs.insert(i);
  for (auto ip : smooths[i].exponents) {
    if (!r.exponents.insert(ip).second) {
      r.exponents.remove(ip);
    }
  }
}

後はこれ全体を大きなビット行列だと考えて依存関係を掃き出します。

vector<Relation> squares;
for (int i = relations.size() - 1; i >= 0; --i) {
  auto& ri = relations[i];
  if (ri.exps.size() == 0) {
    // 指数部分が全て偶数＝従属式
    squares.push_back(ri);
    continue;
  }
  int pid = *ri.exps.begin();
  for (int j = 0; j < i; ++j) {
    auto& rj = relations[j];
    if (rj.exps.contains(pid)) {
      rj ^= rj;  // Relation::operator^= は別途定義するとする
    }
  }
}

vector<vector<int>> candidates;
for (auto& s : squares) {
  candidates.push(vector<int>(s.indexs.begin(), s.indexs.end()));
}

これで candidates には平方数を作る式の組み合わせリストが入っていることになりますので、あとは本当に分解できているかそれぞれ確認するだけです。

for (auto& square : squares) {
  Integer x = 1;
  for (auto& eid : square) {
    x *= smooths[eid].x;
  }

  Integer y = 1;
  set<int> odd_exps;
  for (auto& eid : square) {
    for (auto& pid : smooths[eid].exponents) {
      if (!odd_exps.insert(pid).second) {
        odd_exps.remove(pid);
	y *= factor_base[pid];
      }
    }
  }
  
  x = gcd(x + y, n);
  if (x != 1 && x != n) {
    return x;
  }
  // Fail to factor N.
}

実験

とりあえず一通りできたので実行してみた。実験対象は

11629360743077306442685712558623 \\= 5195515532454019 \times 2238345871633717

でパラメータは

(a, b) = (\sqrt{N}-6000000, \sqrt{N}+6000000)\\ P = 20000

とすると

Factor Base: 1108 primes (0.004 sec)
# of smooth: 1471 (16.506 sec)
# of relations: 403 (78.266 sec)
11629360743077306442685712558623 = 5195515532454019 * 2238345871633717 (0.004 sec)

という感じになりました。少し smooth になる関係式が多いので探索範囲を絞っても良さそうですが、依存関係の算出に 80 秒弱もかかっているのでまずはこれの改善を目指します。→ 実験編 1

選択

実験

Discussion