Subset sumの論文を読んだ感想

実用的に高速

コード

SubsetSumNewton.cpp \( O(n + t\log t) \) と SubsetSumRecursive.cpp \(O(n + t \log^2 t)\)
SubsetSumNewton.cpp · GitHub

概要と意訳

Ce JinさんとHongxun WuさんのA Simple Near-Linear Pseudopolynomial Time Randomized Algorithm for Subset Sumを読んだhttps://arxiv.org/abs/1807.11597。著者を調べるとcodeforcesでLGMだった。
参考文献も読んだhttps://arxiv.org/abs/1004.3412

Subset Sum問題

正整数列\( s_1 \dots s_n \)と正整数\( t \)が与えられる。\( k=1 \dots t \)それぞれについて、\( s \)の部分列で、和が\( k\)になるものの個数を1000000007で割った余りを求めよ。
ただし、同じ部分列でもインデックスが異なれば別のものとして数える。

入力 s = { 1, 1, 2, 3 }, t = 7
出力 2 2 3 3 2 2 1
k=3になるものは {1, 2}, {1, 2}, {3}の3通り

Near-Linear Pseudopolynomial Time Randomized Algorithm

アルゴリズムは時間計算量\( \tilde{O}(n+t) \)。Õはsoft-Oと読み、\( \tilde{O}(n)\)は\(O(n \log^k n)\)の略記らしい。論文中にあるものは\(O(n + t \log^2 t)\)と\(O(n + t\log t)\)。
クラシックな\(O(nt)\)の動的計画法と比べると速い。
時間計算量は入力の長さだけでなく入力の値にも依存する。

「和がちょうどtになる部分列が存在する」ことを判定するにはmodをランダムな素数で試す必要がある。論文で述べられているのは上記の問題を解く決定的アルゴリズムなので、競技プログラマーには場合の数問題のほうがインパクトが大きい気がする。

アルゴリズム

詳しいことは論文を読むとして。解の母関数は\(A(x) =\prod_{i=1}^{n}(1+x^{s_i})\)で愚直に計算すれば動的計画法と同じ。積は対数を取ると和になり、指数で戻せばいいらしい。

計測

f:id:natsugiri:20180817010041p:plain
実行時間
縦軸は秒。横軸xは入力データがn=t=xであることを意味する。
青線newtonは指数関数を参考文献のニュートン法で実装した\(O(n + t\log t)\)アルゴリズム
赤線recursiveは指数関数を論文で紹介された方法で実装した\(O(n + t\log^2 t)\)アルゴリズム
黄線dpは動的計画法\(O(nt)\)アルゴリズム

ニュートン法はt+1を2の冪に切り上げるため階段状になった。
recursiveは分割統治なのでtが2の冪に近いと比較的ロスが少ない。2の冪であるときにニュートン法に追い付いているが、これはニュートン法の計算量は定数倍がやや大きいためだと思われる。

その他

  • 二つのアルゴリズムの違いは指数関数の実装だけだが、どちらも畳み込みで実装されている。こちらのmath314さんのNTTを利用した。

math314.hateblo.jp

  • modが998244353のときはnewtonとrecursiveは3倍速くなる。
  • グラフのrecursiveは細かく測定すれば滑らかになるはず?
  • dpはatcoderのclangだと2~3倍高速。
  • 存在判定問題の場合bitsetを使った\(O(nt/w)\)または\(O(n+t^2/w)\)のアルゴリズムも、現実的な入力サイズに対しては高速(wは32や64など)。
  • mod > t の必要性が言われているが理由がわからなかった。

LISとFenwick Tree

概要

  • Longest increasing subsequence 問題はFenwick treeを用いてO(n log m)で解ける(nは列の長さ、mは要素の最大値)
  • 実行時間としては二分探索を用いる方法と大きな差はない

動機

今朝のCodeforces Round #468の他人のコードを見るとLISをFenwick treeで解いている。それも考察の結果で偶然FenwickになったわけではなくLISはFenwickで解くという強い意志を感じる。実際それで1位を取っている。調べるとCodeforcesの記事GeeksforGeeksの解説がヒットする。つまりこの手法は常識であり、これを知らないままプログラミングコンテストをやっていたのは舐めプだ。

LIS

長さnの数列a[0], ..., a[n-1]が与えられる。数列から0個以上n個以下の要素を消去してできた数列a'[0], ..., a'[k-1]がa'[i] < a'[i+1]を満たすとき、a'を増加部分列と呼ぶ。最長の増加部分列の長さを出力せよ。
1 <= n <= 10^5
1 <= a[i] <= 10^5
a[i]は整数

二分探索による解法

プログラミングコンテストチャレンジブックやWikipediaに載っている。

int buf[MAX_N];
void lis() {
    const int INF = 1<<30;
    for (int i=0; i<N; i++) buf[i] = INF;
    for (int i=0; i<N; i++) *lower_bound(buf, buf+N, A[i]) = A[i];
    int len = lower_bound(buf, buf+N, INF) - buf;
    printf("%d\n", len);
}

INFはa[i]より大きい値であればよい。INFを定義したくない、または初期化をしたくないなら次のコードもよい。

int buf[MAX_N];
void lis() {
    int len = 0;
    for (int i=0; i<N; i++) {
	int k = lower_bound(buf, buf+len, A[i]) - buf;
	buf[k] = A[i];
	if (len == k) len++;
    }
    printf("%d\n", len);
}

Fenwick treeによる解法

Fenwick treeは最大値を出力するクエリの場合

  • 更新はincremental。つまり、より大きい値の更新が来た場合だけ書き換える
  • 出力質問の範囲はprefixのみ

2つの制約が付くがこれでよい。

const int SIZE = MAX_A + 1;
int buf[SIZE];

void update(int i, int x) {
    for (; i<SIZE; i|=i+1) buf[i] = max(buf[i], x);
}

int get_max(int r) {
    int ret = 0;
    for (; r; r&=r-1) ret = max(ret, buf[r-1]);
    return ret;
}

void lis() {
    for (int i=0; i<SIZE; i++) buf[i] = 0;
    for (int i=0; i<N; i++) {
	int k = get_max(A[i]);
	update(A[i], k+1);
    }
    int len = get_max(SIZE);
    printf("%d\n", len);
}

比較

よくある問題設定ならば実行時間はどれも大差ない。しかしa[i]の最大値制約が強い場合(100程度)はFenwick treeの解法が2倍以上速い。
逆に、a[i]の制約が緩い(10^9程度)場合や、a[i]がint型以外の場合は前処理でa[i]を0~n-1に写像しなければならず、やや面倒であり、オンライン問題には使いづらい。

長方形から森構造を作る平面走査アルゴリズム

問題

二次平面上に軸に平行な長方形がn個与えられる。
入力として与えられる長方形について、どの二つの長方形も、

  • 共通部分を持たない
  • 一つの長方形がもう一つの長方形の真に内側にある

のどちらかを満たすことを保証する。

これらの長方形の包含関係から根付き有向森を作りたい。
つまり、「長方形uが長方形vを真に含む」と「森の頂点uは頂点vの先祖である」が同値となるような森。

解法

二次元の場合はx_lowでソートしてxが小さいほうから大きいほうへ平面走査する。y座標を平衡二分木で管理して親頂点を見つける。

typedef vector<int> VI;
typedef long long ValType;

struct Rect {
    ValType x_low, x_high, y_low, y_high;
};

const int MAXN = 100000;
int N;
VI G[MAXN];
Rect R[MAXN];

struct OrderIdByX {
    const Rect *p;
    OrderIdByX(const Rect *p_) { p = p_; }
    bool operator()(const int i, const int j) const {
	return p[i].x_low < p[j].x_low;
    }
};

void forest_from_rects(VI *g, const Rect *rects, int n) {
    if (n == 0) return;

    VI idx(n);
    REP (i, n) idx[i] = i;
    OrderIdByX cmp(rects);
    sort(idx.begin(), idx.end(), cmp);
    set<pair<ValType, int> > se;

    REP (i_, n) {
	int v = idx[i_], p;
	auto it = se.lower_bound(make_pair(rects[v].y_low, -1));
	while (it != se.end() && it->second != -1 && rects[it->second].x_high <= rects[v].x_low)
	    se.erase(it++);
	if (it == se.end() || it->second == -1) p = -1; // v is a root;
	else g[p = it->second].push_back(v); // rects[v] is in rects[p];
	    
	se.emplace_hint(it, rects[v].y_low, p);
	se.emplace_hint(it, rects[v].y_high, v);
    }
}

int main() {
    forest_from_rects(G, R, N);
}

std::set::lower_boundでkeyと異なる型の値を渡す

less<>を使えばよい。

問題

(名, 姓)のpairがたくさんあってset< pair< string,string> >に格納されている。

(Alice, Acer)
(Anna, Cork)
(Bella, Bread)
(Bella, Field)
(Cherry, Card)

名が一つ与えられたときに姓が辞書順最小の人を見つけたい。
例えばBellaが2人いるが、Bellaの姓の辞書順最小のBreadを見つけたい。
setのほかの機能も使いたいのでこれらをsetで行いたい。

解法1

Bellaの一人目を見つけたいのでset::lower_boundが使えそう。
stringは偶然にも最小値を持つデータ型である。
空文字列が最小値なのでmake_pair("Bella", "")でlower_boundすれば見つかる。

解法2

空文字列を使わない方法がある
リファレンスのlower_bound(set::lower_bound - cpprefjp C++日本語リファレンス)を見るとc++14からsetの要素の型とは異なる型で使えるらしい。

template <class K>
iterator lower_bound(const K& x)

しかしこれはオーバーロードするだけではコンパイルできない。

// ERROR
namespace std {
    bool operator<(const pair<string, string> &x, const string &y) {
        return x.first < y;
    }
};

int main() {
    set<pair<string, string> > se;
    se.lower_bound(string("Bella"));
}

リファレンスの備考を見るとset::findを読めとある。findを読むとis_transparentが必要だとある。これはless< void>には定義されている。

つまりエラーの原因は比較関数がless< pair< string, string> >であること。

// GOOD
namespace std {
    bool operator<(const pair<string, string> &x, const string &y) {
        return x.first < y;
    }
};

int main() {
    set<pair<string, string>, less<> > se;
    se.lower_bound(string("Bella"));
}

解法3

自分で構造体にis_transparentを定義してもよい。

// GOOD
struct Cmp {
    bool operator()(const pair<string, string> &x, const pair<string, string> &y) const {
        return x < y;
    }

    bool operator()(const pair<string, string> &x, const string &y) const {
        return x.first < y;
    }

    typedef void is_transparent;
};

int main() {
    set<pair<string, string>, Cmp> se;
    se.lower_bound(string("Bella"));
}

感想

set< string>と文字列リテラルを比較するための仕組みらしいが、このように使ってもいいものかどうか。

SPCLN, SnackDown 2017 Online Elimination Round 解法案

SnackDown 2017 Online Elimination Round
Cleaning the Space (SPCLN)

https://www.codechef.com/SNCKEL17/problems/SPCLN

この問題は
https://www.codechef.com/problems/RINと同値。問題設定や変数名を変えただけ。Editorialもある。

解法

最小カット問題として解く。
とりあえず各space debrisの最大エナジーをとり、他のボタンで破壊するならその分のペナルティを受けると考える。
依存関係を保つために流量INFを張ってペナルティを最小カットとして求める。
グラフは「燃やす埋める問題」を多層にした感じ。RINのEditorialが詳しい。

PREFIXOR, SnackDown 2017 Online Elimination Round 解法案

SnackDown 2017 Online Elimination Round
Prefix XOR (PREFIXOR)

https://www.codechef.com/SNCKEL17/problems/PREFIXOR
https://www.codechef.com/viewsolution/14319736

問題

非負整数列A[1]..A[n]が与えられる。次のオンラインクエリに答えよ
(l, r):次の条件を全て満たすペア(i, j)の個数を出力せよ。

  1. l<=i<=j<=r
  2. 任意のk (i<=k<=j-1)について、(A[i]^A[i+1]^...^A[k]) <= (A[i]^A[i+1]^...^A[k+1])。^はxor。

n, q <= 400000

続きを読む

BLACKCOM, SnackDown 2017 Online Elimination Round

SnackDown 2017 Online Elimination Round 解法案
Black Nodes in Subgraphs (BLACKCOM)

https://www.codechef.com/problems/BLACKCOM/
https://www.codechef.com/viewsolution/14319524

問題

入力でN頂点の木が与えられる。各頂点は白か黒のどちらかの色が塗られている。Q個のクエリに答えよ。
クエリ(b, s) : 頂点数がちょうどs、そのうち黒頂点の数がちょうどbの部分連結グラフは存在するか。yes/noで答えよ。
N <= 5000
Q <= 100000

解法

動的計画法/木DP
全ての頂点について、部分木のdpテーブルを4つ作る。

dp1_min[v][s] := 部分木vに対して、「vを含むサイズsの連結部分グラフ」が含む黒頂点の最小の数
dp1_max[v][s] := 同上の最大の数
dp2_min[v][s] := 部分木vに対して、「vを含まないサイズsの連結部分グラフ」が含む黒頂点の最小の数
dp2_max[v][s] := 同上の最大の数

v=0を根とするとクエリ(s, b)に対する答えは

min(dp1_min[0][s], dp2_min[0][s]) <= b && b <= max(dp1_max[0][s], dp2_max[0][s])

である。
ボトムアップにdpテーブルを作る。(頂点v, vの子w)があるとき、wのテーブルを見てvのテーブルを更新。

sz[v] = 1;
EACH (e, G[v]) {
    for (int i=sz[v]; i>=1; i--) for (int j=sz[*e]; j>=1; j--) {
	amin(dp1_min[v][i+j], dp1_min[v][i] + dp1_min[*e][j]);
	amax(dp1_max[v][i+j], dp1_max[v][i] + dp1_max[*e][j]);
    }
 
    REP (j, sz[*e]+1) {
	amin(dp2_min[v][j], dp1_min[*e][j]);
	amin(dp2_min[v][j], dp2_min[*e][j]);
	amax(dp2_max[v][j], dp1_max[*e][j]);
	amax(dp2_max[v][j], dp2_max[*e][j]);
    }
 
    sz[v] += sz[*e];
}

dp1の更新が二重ループになっているが木全体の計算量はO(N^2)になる。

部分木wのサイズをR、部分木vのすでに見た頂点の数をLとしてN=L+Rのステップ数をもとめる(つまりコード中で、sz[v]==L, sz[*e]==Rのとき)。T(N)上記のdpの計算ステップ数とする。
T(N) >= T(L) + T(R) + L*R
N^2 = (L+R)^2 >= L^2 + R^2 + L*R
よってT(N) = O(N^2)