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)

東京工業大学でTarjan教授を拝見した

講義に出席した。顔を拝むのが目的。
f:id:natsugiri:20170518155835j:plain:w500

Tarjan先生

我々が毎日書いているUnion-Findや強連結成分分解(lowlinkのやつ)の証明やらアルゴリズムやら考えたすごい人。
1時間の講義でZip Treeをマスターしよう。

Zip Tree

二分探索木の一つ。詳細な解説は講義に出た人の特権なので避けるがTreapにとても似ている。正直、名前聞いたことなかった。
Treapと比較して2つの違いがある

ランク

Treapはランダムな値を優先度としてノードに割り当てるがzip treeではそれをランクと呼ぶ。Treapでは経験上ノード1つに対してO(log n)ビットの追加メモリが必要だが、Zip TreeはO(log log n)ビット、多分。。。
実装ではこんな感じでいいだろうか?Treapのランクが32ビットなのに対し、Zip Treeは31以下(5ビット)で済んでいる。Treapで乱数の下位5ビットにすると壊滅する。

mt19937 engine;

int zip_tree_rnk() {
    unsigned x = engine() | (1u<<31);
    return __builtin_ctz(x);
}

int treap_rnk() {
    return engine();
}

insert / delete

名にZipを冠するがinsert / deleteの操作から来ているようだ。トップダウンでできるのでループ処理で実装するのが簡単なはず。実装は一番下。

実装・実行

Insert 200000

Treap
Height 42; Ave. depth 21.255840;

Zip Tree
Height 52; Ave. depth 22.656550;

実装量はほぼ同じ。実装安を目指したので値の重複を許してinsertされる再帰関数にしてしまったためZipTreeがやや遅いかもしれない。ランクは結局どちらもintで持たせたので省メモリの恩恵を受けていない。
上記の実装上の2つの違いは独立で、Treapの実装で1つだけ変更を加えてやるとかもできる。

実装は上半分がTreap、下半分がZip Tree。より良い実装はそのうち考える。

#include<climits>
#include<cassert>
#include<set>
#include<random>
#include<stdio.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<string.h>
using namespace std;

typedef long long LL;

//////////////////////////////////////////////////
// Common;
//////////////////////////////////////////////////

mt19937 engine(19480430);

struct Node {
    int val;
    int rnk;
    Node *l, *r;
    Node(int val_, int rnk_) : val(val_), rnk(rnk_), l(0), r(0) {}
};

typedef Node* Tree;

Tree find(Tree t, int val) {
    while (t) {
	if (val < t->val) {
	    t = t->l;
	} else if (t->val < val) {
	    t = t->r;
	} else {
	    return t;
	}
    }
    return 0;
}


//////////////////////////////////////////////////
// Treap;
//////////////////////////////////////////////////

// Subfunctions;
Tree rotate_right(Tree x) {
    Tree y = x->l;
    x->l = y->r;
    y->r = x;
    return y;
}

Tree rotate_left(Tree x) {
    Tree y = x->r;
    x->r = y->l;
    y->l = x;
    return y;
}

Tree treap_erase_minimum(Tree t, Tree y) {
    if (!t->l) {
	y->val = t->val;
	y->rnk = t->rnk;
	Tree x = t->r;
	delete t;
	return x;
    } else {
	t->l = treap_erase_minimum(t->l, y);
	return t;
    }
}

Tree treap_balance(Tree t) {
    if (t->l && t->rnk <= t->l->rnk && (!(t->r) || t->l->rnk >= t->r->rnk)) {
	t = rotate_right(t);
	t->r = treap_balance(t->r);
    } else if (t->r && t->rnk < t->r->rnk) {
	t = rotate_left(t);
	t->l = treap_balance(t->l);
    }
    return t;
}

// Main functions;
int treap_rnk() {
    return engine();
}

Tree treap_insert(Tree t, int val, int rnk) {
    if (!t) {
	return new Node(val, rnk);
    } else {
	if (val < t->val) {
	    t->l = treap_insert(t->l, val, rnk);
	    return (t->l->rnk >= t->rnk? rotate_right(t): t);
	} else {
	    t->r = treap_insert(t->r, val, rnk);
	    return (t->r->rnk > t->rnk? rotate_left(t): t);
	}
    }
}

Tree treap_erase(Tree t, int val) {
    if (!t) {
	return 0;
    } else {
	if (val < t->val) {
	    t->l = treap_erase(t->l, val);
	    return t;
	} else if (t->val < val) {
	    t->r = treap_erase(t->r, val);
	    return t;
	} else if (!t->l || !t->r) {
	    Tree x = (t->l?: t->r);
	    delete t;
	    return x;
	} else {
	    t->r = treap_erase_minimum(t->r, t);
	    return treap_balance(t);
	}
    }
}

//////////////////////////////////////////////////
// Zip Tree;
//////////////////////////////////////////////////

// Subfunctions;
Tree unzip_insert_root(Tree t, int val, int rnk) {
    Tree x = new Node(val, rnk);
    Tree low = x, high = x;
    while (t) {
	if (val < t->val) {
	    high = high->l = t;
	    t = t->l;
	} else {
	    low = low->r = t;
	    t = t->r;
	}
    }
    low->r = high->l = 0;
    swap(x->l, x->r);
    return x; 
}

Tree zip_erase_root(Tree t) {
    bool b = false;
    Tree x = t;
    Tree low = t->l, high = t->r;
    while (low && high) {
	if (low->rnk < high->rnk) {
	    (b? x->r: x->l) = high;
	    x = high;
	    b = false;
	    high = high->l;
	} else {
	    (b? x->r: x->l) = low;
	    x = low;
	    b = true;
	    low = low->r;
	}
    }
    (b? x->r: x->l) = (low?: high);
    x = t->l;
    delete t;
    return x;
}

// Main functions;
int zip_tree_rnk() {
    unsigned x = engine() | (1u<<31);
    return __builtin_ctz(x);
}

Tree unzip_insert(Tree t, int val, int rnk) {
    if (!t) {
	return new Node(val, rnk);
    } else if (val < t->val) {
	if (t->rnk <= rnk) {
	    return unzip_insert_root(t, val, rnk);
	} else {
	    t->l = unzip_insert(t->l, val, rnk);
	    return t;
	}
    } else {
	if (t->rnk < rnk) {
	    return unzip_insert_root(t, val, rnk);
	} else {
	    t->r = unzip_insert(t->r, val, rnk);
	    return t;
	}
    }
}

Tree zip_erase(Tree t, int val) {
    if (!t) {
	return 0;
    } else if (val < t->val) {
	t->l = zip_erase(t->l, val);
	return t;
    } else if (t->val < val) {
	t->r = zip_erase(t->r, val);
	return t;
    } else {
	return zip_erase_root(t);
    }
}

//////////////////////////////////////////////////
// Debug;
//////////////////////////////////////////////////
double depth_sum;
int height(Tree t, int depth=0) {
    if (t) {
	depth_sum += depth;
	return max(height(t->l, depth+1), height(t->r, depth+1)) + 1;
    } else {
	return -1;
    }
}

void show(Tree t, int depth=0) {
    if (t) {
	show(t->l, depth+1);
	for (int i=0; i<depth; i++) printf("..");
	printf("%d\n", t->val);
	show(t->r, depth+1);
    }
}

int last;
void verify(Tree t) {
    if (t && t->l) {
	assert(t->l->val <= t->val);
	assert(t->l->rnk < t->rnk);
	verify(t->l);
    }

    assert(last <= t->val);
    last = t->val;

    if (t && t->r) {
	assert(t->val <= t->r->val);
	assert(t->rnk >= t->r->rnk);
	verify(t->r);
    }
}

const int SIZE = 200000;
int A[SIZE];
int cnt[32];

int main() {
    int cc = 0;
    printf("Insert %d\n", SIZE);
    puts("");

    for (int i=0; i<SIZE; i++) A[i] = engine();

    // Empty Tree;
    Tree treap = 0;

    for (int i=0; i<SIZE; i++) {
	int k = treap_rnk();
	treap = treap_insert(treap, A[i], k);
	cc++;
    }
//    for (int i=0; i<SIZE/2; i++) {
//	treap = treap_erase(treap, A[i]);
//	cc--;
//    }

    last = INT_MIN;
    verify(treap);

    printf("Treap\n");
    int h = height(treap);
    printf("Height %d; Ave. depth %f;\n", h, depth_sum / cc);
    puts("");

    // Empty Tree;
    Tree zipTree = 0;

    cc = 0;
    for (int i=0; i<SIZE; i++) {

	int k = zip_tree_rnk();
	cnt[k]++;
	zipTree = unzip_insert(zipTree, A[i], k);
	cc++;
    }
//    for (int i=0; i<SIZE/2; i++) {
//	zipTree = zip_erase(zipTree, A[i]);
//	cc--;
//    }

    // show(zipTree);

    last = INT_MIN;
    verify(zipTree);

    printf("Zip Tree\n");
    depth_sum = 0;
    h = height(zipTree);
    printf("Height %d; Ave. depth %f;\n", h, depth_sum / cc);
    puts("");
    puts("");

    printf("rnk distribution\n");
    for (int i=0; i<32; i++) printf("%d%c", cnt[i], i==31?'\n':' ');


    return 0;
}