二次篩法 (QS) (高速化)
二次篩法 (QS) (2)の続きです。
ひとまず計算はできるようになりましたが、計算時間がかかり過ぎるのが難点です。ザッと時間配分を見ると依存関係式の算出部分に多くの時間を割いていることがわかります。パラメータを変更してその割合を減らすことも可能ですが、今回はプログラムの変更で高速化を目指しましょう。
依存関係式の算出部分で気になるのが利用する数式の情報の持ち方です。試しに出力すると403の関係式のうち393を利用して計算している（ケースがある）とのことなので set で 1 個ずつ持つのが非効率なのは明らかです。なのでここの高速化を試みてみましょう。
式番号管理の変更
なのでまずは64bitずつビット情報に圧縮して計算するようにしてみましょう。式番号
struct Relation {
- set<int> indexs; // indecies にするとパターン化がメンドイ
+ map<int, uint64> indexs;
set<int> exponents;
};
// 式情報をビット圧縮情報に変換
for (int i = 0; i < smooths.size(); ++i) {
Relation r;
- r.indexs.insert(i);
+ r.indexs[i / 64] = 1ULL << (i % 64);
for (auto ip : smooths[i].exponents) {
vector<vector<int>> candidates;
for (auto& s : squares) {
- candidates.push_back(vector<int>(s.indexs.begin(), s.indexs.end()));
+ candidates.push_back({});
+ auto& candidate = candidates.back();
+ for (p : s.indexs) {
+ for (int i = 0; i < 64; ++i) {
+ if ((p.second >> i) & 1) {
+ candidate.push_back(p.first * 64 + i);
+ }
+ }
+ }
}
では計算時間を見てみましょう。（1回勝負）
Factor Base: 1108 primes (0.003 sec)
# of smooth: 1471 (11.549 sec)
# of relations: 403 (24.889 sec)
11629360743077306442685712558623 = 5195515532454019 * 2238345871633717 (0.003 sec)
ということで依存性解析の部分が3倍速くなりました。
指数管理の変更
変更による高速化が思った以上に大きかったので調子に乗って指数部分も同様の管理にしてみました。
struct Relation {
map<int, uint64> indexs;
- set<int> exponents;
+ map<int, uint64> exponents;
};
// 式情報をビット圧縮情報に変換
for (int i = 0; i < smooths.size(); ++i) {
Relation r;
r.indexs[i / 64] = 1ULL << (i % 64);
for (auto ip : smooths[i].exponents) {
- if (!r.exponents.insert(ip).second) {
- r.exponents.remove(ip);
+ int index = ip / 64;
+ r.exponents[index] ^= 1ULL << (ip % 64);
+ if (r.exponents[index] == 0) {
+ r.exponents.remove(index);
}
}
}
auto pid = *ri.exponents.begin();
+ uint64 bit = pid.second & (~pid.second + 1)
for (int j = 0; j < i; ++j) {
auto& rj = relations[j];
- if (rj.exponents.contains(pid)) {
+ if (rj.exponents.contains(pid.first) && (rj.exponents[pid.first] & bit)) {
rj ^= rj; // Relation::operator^= は別途定義するとする
}
}
この形で実験した結果、
Factor Base: 1108 primes (0.003 sec)
# of smooth: 1471 (12.771 sec)
# of relations: 403 (6.110 sec)
11629360743077306442685712558623 = 5195515532454019 * 2238345871633717 (0.003 sec)
と依存性解析の部分が更に4倍速くなりました。元のと比べると13倍ですね。
式管理の密行列化
さて冒頭の導入部分で
試しに出力すると403の関係式のうち393を利用して計算している（ケースがある）
とありますが、これは素直に考えると map(疎行列が前提) で管理するよりも vector(密行列が前提) で管理する方が効率的な気がします。ただ依存性解析の序盤では明らかに疎行列なので段階的に密行列へ移行させることができるといい気がします。
なので今回の実装では
- 密行列は 0 番目から順番に適用範囲を伸ばしていく
- 密行列の次の 64bit に当たる情報を半分の式が保持していたら移す
という形で実装してみます。 LLL？知らないですね。
struct Relation {
+ vector<uint64> indexd; // 密行列部分
map<int, uint64> indexs;
map<int, uint64> exponents;
};
vector<Relation> squares;
for (int i = relations.size() - 1; i >= 0; --i) {
... // 指数のbitを消す処理
+ int nextd_index = ri.indexd.size(); // 次にindexdに移すかもしれないindex
+ int num_using = 0; // ↑を使ってるrelationの数
+ for (int j = 0; j < i; ++i) {
+ num_using += relations[j].indexs.contain(nextd_index);
+ }
+ if (num_using * 2 >= i) {
+ // indexs -> indexd の移行
+ for (int j = 0; j < i; ++i) {
+ auto& sparse = relations[j].indexs;
+ uint64 bits = sparse.contain(nextd_index) ? sparse[nextd_index] : 0;
+ relations[j].indexd.push_back(bits);
+ }
+ }
}
vector<vector<int>> candidates;
for (auto& s : squares) {
- candidates.push(vector<int>(s.indexs.begin(), s.indexs.end()));
+ candidates.push_back({});
+ auto& candidate = candidates.back();
+ const auto& dense = s.indexed;
+ for (int i = 0; i < dense.size() * 64; ++i) {
+ if ((dense[i / 64] >> (i % 64)) & 1) {
+ candidate.push_back(i);
+ }
+ }
}
まぁまたメインっぽい作業は operator^= の中なんですが。
とりあえず再度時間計測すると
Factor Base: 1108 primes (0.003 sec)
# of smooth: 1471 (13.779 sec)
# of relations: 403 (5.684 sec)
11629360743077306442685712558623 = 5195515532454019 * 2238345871633717 (0.004 sec)
ということでちょっぴり高速化しましたが、誤差程度、もしくはデータ依存の可能性が高そうですね。ということで次は複数データでの検証をできるようにしたいです。→実験編
Discussion