ディズニーデラックスを使い始めた

英語のリスニングを向上させるという邪な理由でDisney deluxeを見始めた。31日間無料、以降756円/月。まずはMarvelを古い順に観る。
音声・字幕それぞれ英語・日本語が好きな組み合わせで利用できる。ビデオプレーヤーはシークバー、30秒進める・戻る、フルスクリーン機能がある。再生速度や画質の選択はない。

CODE FESTIVAL 2016 Grand Final(Parallel) J - 123 Pairs 解法

CODE FESTIVAL 2016 Grand Final(Parallel) J - 123 Pairs 解法
https://atcoder.jp/contests/cf16-exhibition-final-open/tasks/cf16_exhibition_final_j

問題

1から2NをN個のペアに分割したい。このとき、差が1のペアA個、差が2のペアB個、差が3のペアC個、N=A+B+Cを満たすものはいくつあるか?

解法

2種類のペアを並べたときに、隙間がなく、それ以上分割できない極小なものを列挙すると、

  • Pattern 1 : 差1のペア1個
  • Pattern 22-k : 差2のペア2個と、差3のペアk個 (k >= 0)
  • Pattern 31 : 差1を1個、差3を1個
  • Pattern 333 : 差3を3個

以上4通り。
Bが奇数の場合は差2が余るので解0。
Pattern 31をs個、Pattern 333をt個使うとすると、Pattern 1 は A-s個。Pattern 22-kは常にB/2個。
これらの並べ方は重複順列である。
\displaystyle \frac{(A+t+B/2)!}{s!t!(A-s)!(B/2)!}

差3の残りC-s-3t個をPattern 22-kの隙間に振り分ける。これは重複組合せ。Bが0の場合だけうまくいかないので注意する。
\displaystyle
\left(
\begin{array}{c}
C-s-3t+B/2-1  \\
B/2-1
\end{array}
\right)

s, tの2重ループでO(n^2)解が求まる。

sとtを分解したい。

\begin{array}{rcl}
f(x) &=& \sum_s \displaystyle \frac{1}{s!(A-s)!} x^s \\
g(x) &=& \sum_t \displaystyle \frac{(A+t+B/2)!}{t!(B/2)!} x^{3t} \\
(f*t)(x) &=& \sum_{u=s+3t} \displaystyle  \frac{(A+t+B/2)!}{s!t!(A-s)!(B/2)!} x^u
\end{array}

f,gの積を高速フーリエ変換と畳み込みで求めてO(n log n)。

CODE FESTIVAL 2016 Grand Final H - AB=C Problem 解法

CODE FESTIVAL 2016 Grand Final H - AB=C Problem
https://atcoder.jp/contests/cf16-exhibition-final-open/tasks/cf16_exhibition_final_h

問題

正方行列Cが与えられる。C=ABを満たすような正方行列のペア(A, B)はいくつあるか求めよ。
行列の要素の演算はすべて2で割った余りを取る。

解法

CとAに対して、行基本変形してもC=ABと(A, B)の個数が維持される。
CとBに対して、列基本変形しても同様。よって解はCのサイズNとランクRにのみ依存する。
Cのランクはガウスの消去法で求めることができる。

Cのj列目のベクトルは、A = (a1, ..., an) のN個の列ベクトルに係数を掛けて足したもので、その係数はBのj列目のベクトルである。
AのN個のベクトルがCの空間を張り、AのランクがSならば、Bのj列目の選び方はpow(2, N-S)通りである。N列あるので、Aが決まったあとのBの選び方はpow(2, N*(N-S))通り。

Cの空間を張るようなAを数える。

dp[i][j][k] := i個のベクトル(a1, ..., ai)の数、ただし

  1. (a1, ..., ai)のランクがj
  2. (a1, ..., ai)とCの共通部分のランクがk

初期値

dp[0][0][0] = 1;

i+1個目のベクトルの pow(2, N)通りは:
(1) j個の基底から作る場合 pow(2, j) 通り

dp[i+1][j][k] += dp[i][j][k] * pow(2, j);

(2) Cとの共通部分がないような線形独立なベクトルを作る場合、「Cの空間とj個の基底」から作ることができないベクトル。pow(2, N) - pow(2, R-k+j) 通り

dp[i+1][j+1][k] += dp[i][j][k] * (pow(2, N) - pow(2, R-k+j));

(3) Cとの共通部分がある線形独立なベクトル。pow(2, R-k+j) - pow(2, j) 通り

dp[i+1][j+1][k+1] += dp[i][j][k] * (pow(2, R-k+j) - pow(2, j));

Bの選び方を考えると、最終的な解は

for (j = R; j <= N; j++)
    ans += dp[N][j][R] * pow(2, N*(N-j));

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