2016-09-11

SA-IS

Suffix Arrayを空間・時間計算量共に線形で構築出来るというアルゴリズムSA-ISを, 以前ふわっと読んだことはあったけれど真面目に読みなおしました. よく出来過ぎだと思う. 細かいところの理解が間違っていたら教えて下さい. ざっくりとした流れしか書きません. ちゃんとした証明とかは元論文を見てください. というかこんな記事よりiwiさんのブログとかがわかりやすいと思うけど自分で忘れないために書きました.

(元論文: Linear Suffix Array Construction by Almost Pure Induced-Sorting)

文字列の最後には最小の番兵を置いておく.

記号を色々定義しておくと,まず
Suf(S,i):=Sのi文字目以降のsuffix
Suf(S,i)がS-type...Suf(S,i)Suf(S,i+1)
Suf(S,i+1)とはSuf(S,i)の先頭を除いたものであることに注意(当たり前だけどはっきり意識した方がわかりやすいと思った)

LMS...Suf(S,i)がS-typeかつSuf(S,i-1)がL-type
LMS-substring...LMSから次のLMSまでの部分文字列(両端含む) また最後の番兵はそれだけでLMS-substringということにする

ここで,LMSたちのsuffix arrayが求まっていたら全体のsuffix arrayが計算出来ることがわかる(Induced-sorting)
まず,考えると出来上がったsuffix arrayでは(先頭の文字,type)でブロックに分けることが出来るのがわかる(L-typeがS-typeの前に来る)

とりあえずLMSたちをその小さい順にブロックに入れていく. LMSの順番だけが重要だとわかるので場所は確定出来なくてよい.
次に,SA[i]を小さい方から見ていって,SA[i]が埋まっていてSA[i]-1がL-typeならばS[SA[i]-1]のブロックの後ろにSA[i]-1を入れる.
というのもSA[i]-1がL-typeならばSA[i]はL-typeまたはLMSで,いずれの場合もSA[i-1]はSA[i]より大きく,またs < tの時'x'+s < 'x'+tであるから小さい方からやるとうまくいく. これでL-typeのsuffix arrayが求まる.
次に,S-typeを埋める.これは今度はSA[i]を大きい方から見てほぼ同じことをする.LMSを入れた場所は上書きしていいことに注意.

以上より,LMSがソート出来ていればsuffix arrayが求まる.

次に,上の手順を行うためにLMSたちをソートする方法を考える.
これは上で定義したLMS-substringたちのランクを文字として扱って(LMS-substringを比較する時には文字にtypeの属性も含める), suffix arrayを求めれば(ここで再帰が呼ばれる)出来ることがわかる.
よってLMS-substringの順位付けを考える. ここにもInduced-sortingが使えるというのがSA-ISのすごいところらしい.

実はこれは,LMSたちの先頭でとりあえずバケットに分けて(同じバケットの中では開始位置が後ろのやつから入れる) induced-sortingすればよい.
これはその文字列のk文字目までソートしてあるとk+1文字目までの情報が...という感じになるのでLMS-substringたちが正しい順番で並んでくれるため
は?という感じだがよく考えると上手く出来ている
sortされたLMS-substringたちのrenameは,まあLMS-substringの和はN未満なので愚直にやってよい(1つ前と同じかどうか調べる)

おしまい

参考:
SA-IS - (iwi) ｛反省します - TopCoder部
 SA-IS - Shogo Computing Laboratory

2016-09-06

JAG夏合宿2016参加記

ICPC

参加記を書き書き(典型)

去年に引き続きJAGの夏合宿に参加してきたので感想を書きます(合宿代安くて有り難すぎる)

ICPCアジア地区のチーム(catsatmat)のうち自分だけの参加となりました(悲しい)

・day1
合宿初日、コンテストは無し。懇親しなかったのでtozangezan/DEGwerとMTKに行きました。あとyosupoさんがいた

・day2
朝食を食べ損なう。この日はsnukeさんとatahunというチームで出ました。snukeさんはこのセットはOpencupで参加したことがあったらしく記憶がある問題は自分が考えることに。

自分がまずCを書く。朝に弱いのでクソコードを書いてRE。
コードを見てもらっている間に朝なのでH(Good morning!)を通す。
Iが簡単そうだったので書く。WA。(見た目はミスってなさそうだったのでsnukeさんに後でchallengeケースを作ってもらってなるほどなあ...ってなった。実際簡単だけど丁寧にやらないとダメなことが後でわかった)
Iを直す前にCの方針を変えて綺麗に書くと通る(ごめんなさい...)
Dが簡単とのことなのでsnukeさんに通してもらう。
Iをちょっと直してAC。
Fはrng問題のdictionaryみたいな感じだけど3種類しかないので愚直にdpをヒイヒイ言いながら書く。実際TLがきつかったので前計算出来るところをするようにするもTLE。もしやと思ってdpの初期化パートをmemsetから真面目に必要な分だけすることにするとAC。テストケース数を明示してほしい....
Lを考えているとKがわかったとsnukeさんが言うので書いてもらって単独AC。sugoy
Lが区間dpで解けることがわかるも時間切れ。なんか細かいところをすぐミスるなあ...という感じ。Fに時間を吸われすぎた。夜はsnukeさんとMTKに行った。

・day3
朝食を食べ損なう。(2回目)この日はヒカリエのLINE会場。ICPCチームの指圧マットが来てくれたので2人で出ました。

Aを指圧マットに読んでもらって簡単そうということで提出(難易度順なのかなと思った)。WA。数字が連続していないといけないのを見落としていたらしい。直してAC。
その間に自分がBを読んでゴールまでの最短距離が他のより小さいならOKということがわかったので書く。簡単だけど綺麗だと思った。
AC。
Cが地獄っぽいので指圧マットに任せて(最低)D以降を読んで考えていく。
Dを読むが少しややこしそう...なのでEも読むがまともな解法がありそうに見えなかった。F,Gも読むがほぼ思考せず。幾何は余裕がなさそうかなあという感じだった。(結果的には簡単でしたね...)
そうこうしているうちにCが書けたというので提出してWA。実はサンプル1つ答えが合ってなかったらしいがコンテスト終了まで気付かず大反省。Cのコードをチェックしていてもらう。
Dがわかったつもりになったので書いて提出してWAを貰う。2/3くらいは通っていたのでなんでやという気分になる。
どこがまずいかわからないしEはハッシュでやると子同士で重ねあわせ出来ていい感じっぽいなあと思ったので愚直ハッシュを書いて送ったらACが貰えてびっくりした。
Dの方針を説明して指圧マットにも考えてもらう。自分は")))...)"の列に"("たちを挿入してコストがどうなるかを考察した感じになった(指圧マットはグリッド上の経路として考えていたが結果的な方針は同じだと思ったが...
数分で書けるようなコードだったので指圧マットの方針で書いてもらってAC。激ごめんなさいという感じに...
まだCのWAが取れていなかったのでCを見直してもらってFを考える。最後に使うやつを全部試して後を高速化する方針だとRMQが必要な感じになった。(解説の方針の方がバグりにくそうだったなあ.....)添字が1ずれたりして苦戦しながらなんとかAC。
結局Cはiteratorのインクリメントの前置後置に問題があったらしい。かなしいなあ... サンプルのWAに気付かないのは流石にダメ。

夜はsugim48さん作のAGC004があった。開いてBを少し考えるも問題の解釈とサンプルが合わないのでDを考える。木なので葉からdfsするだけで良いとわかる。Cを考えるもわからない。Bを読み直すと意味がわかる。綺麗。Aを書く。Cはどう頑張ってもわからない....これはまあわからない時はわからないやつだなあという感じになったのでEを考える。状態HHWWなdpが出来そうなことがわかるが遷移にこまる。終了(絶望)

夜MTKにみんなで行くも閉まっていて悲しかった。

・day4
朝食を食べ損なう(3回目) Klabでラブラブ。この日はHecさんとYazatenさんとチームMILK_TEAで出た。(メモ:2015-2016 XVI Open Cup, Grand Prix of Bashkortostan, SKB Kontur Cup Stage 2)

時系列で書くの疲れたので...
A: ColorfulBoardじゃんと思ったとみんな言っていて面白かった。
B: 無理ゲーらしい
C: 公比1かそれ以外で場合分けするも公比1の時1+1+...+1以外のパターンがあるのになかなか気付かずハマった。ごめんなさい...
D: 2分探索出来ないのに注意。1/iの和がNlogNなのでRMQを書いてAC。
E: 読んでない。
F: 問題を見た感じ嘘をつきやすそう。
G: 読んでない
H: ICPCという感じの問題。紙コーディングで詰めたの初めてだった。暦の仕様を読み飛ばして1WAもらった(なんで100日あるんや)
I: 手元でテストすると40回も質問すればx1が一意に定まることがわかったので特定した後はやるだけ。
J: しばらく考えたらdp[i][j]:i番目まで見てAの和=jの時のBの和の最大値としてdpしてj < pについてdp[n][j]がqより小さければOK
(p,qとA,Bをswapして2回やる) 復元が面倒...
K: しらない
L: aの長さで漸化式が立つことがわかる. 意外とiが膨れ上がることがあるらしくmapとlong longにしたら通った。ちゃんとテストしようね...
M: しらない

なんか家でやる時とオンサイトなコンテストだと全然速度が違って悲しくなる3日間だった。まあ焦るしね

さいごに: クソなぞなぞ: ~ に早く飽きないとまずい

2016-08-02

Atcoder Grand Contest No.2

AtCoder

コンテストを寝過ごしたので解いた. Dの一般的なテク感すき. きれい. Bが面白い感あった

Eはとりあえずグリッド上を動くゲームなことまではわかったけどあんま考えてない.

Fの方針が解説と少し違ったので書いておこう.

1~Nの1番左のもの(0は無視)がこの順番にあることにして最後にN!をかける. ここで0は出てくる順番に1~Nと結びつけて良いことがわかる

よって0NN..NN (NはK-1個)から始めて, 0(N-1)(N-1)(N-1)...(N-1), ...., 01...1を間に追加していくことが出来ないか考える.

まず上の考察により新しい0は1番左に. そして0iii,,,iを足す時, 1番左のiは, 現在の列の1番左の(i+1)より左, つまり0でない要素のうち1番左より左に入れないといけないことがわかるが, これは先頭の0の個数を状態として持てばOK.

よってdp[i][j]: 大きい方からi種類追加して先頭にj個0があるような条件を満たす並べ方の個数とおくと, 状態はO(N^2)で
遷移も1番左のiを置く場所を全部試すことで遷移O(N), 合計O(N^3)で計算可能. (実際は二項係数を求めるのにlogがつく)

しかし, 係数の式を考えるとそれは新しい1番左のiの場所にしか依存しないことがわかるので, dp[i][j]のjの大きい方から累積和を持っておくことでO(N^2log hoge)で解くことができる.

解説みたいにグラフにするとわかりやすくて感動した.

追記: 1<=i<=Mについてi!の逆元を求めたい時, M!の逆元を求めた後逆向きにかけていけばO(N^2)で計算出来るの今まで気付いてなかった...というわけでO(N^2)で解けます

#include <bits/stdc++.h>
using namespace std;

typedef pair<int, int> pii;
typedef long long ll;
typedef vector<int> vi;

#define pb push_back
#define eb emplace_back
#define mp make_pair
#define fi first
#define se second
#define rep(i,n) rep2(i,0,n)
#define rep2(i,m,n) for(int i=m;i<(n);i++)
#define ALL(c) (c).begin(),(c).end()

const ll MOD = 1000000007;

ll f[4000010];
ll invf[4000010];

ll mod_pow(ll x, ll k)
{
    ll res = 1;
    for (; k; x = x * x % MOD, k /= 2) if (k & 1) res = res * x % MOD;
    return res;
}

ll comb(int n, int k)
{
    if (k < 0 || k > n) return 0;
    return f[n] * invf[k] % MOD * invf[n-k] % MOD;
}

int N, K;

ll dp[2010][2010]; // j = num of front 0
ll acm[2010];

int main() {
    f[0] = invf[0] = 1;

    for (int i = 1; i <= 4000000; ++i) {
	f[i] = f[i-1] * i % MOD;
	invf[i] = mod_pow(f[i], MOD - 2);
    }
    
    cin >> N >> K;

    if (K == 1) {
	puts("1");
	return 0;
    }

    dp[1][1] = 1;

    for (int i = 2; i <= N; ++i) {
	memset(acm, 0, sizeof(acm));

	for (int j = 2000; j >= 0; --j) {
	    acm[j] = (acm[j + 1] + dp[i-1][j]) % MOD;
	}

	for (int j = 1; j <= i; ++j) {
	    int pos = (i - 1) * K - j + 1;

	    ll t = comb(pos + K - 2, K - 2);
	    dp[i][j] = (dp[i][j] + acm[j - 1] * t) % MOD;
	}
    }

    ll ret = 0;
    rep(i, 2010) ret = (ret + dp[N][i]) % MOD;
    ret = ret * f[N] % MOD;

    cout << ret << endl;

    return 0;
}

2016-07-26

Codeforces #296 Div1 D. Fuzzy Search

Codeforces

問題:
http://codeforces.com/contest/528/problem/D

問題概要は省略します.

典型なのかもしれないけど, あんまりこういうの解いたことなかったので面白かった

解法:

流石にまとめてマッチ箇所を求めるのは大変そうなので, 置く場所は全部試すことにする.
そこで判定を速く出来ればそれでOK

文字が4種類しかないので, 各場所で各文字について判定出来ればOK

まず, Sのi番目に文字Pをあわせるとき, [i-k,i+k]番目にPがあれば良い.
Sの各場所にPを合わせてよいかどうかの配列はimos法などで線形で計算出来る.
そしてTについても各文字がPか否かのbit列にする.

左端をi (0 <= i < |S|-|T|)番目に置いた時, bit演算風に書くとS[i,i+|T|-1] & T[0,|T|-1]がT[0,|T|-1]に等しければその文字Pについてマッチ可能.

よく考えるとこの処理は内積みたいな感じで,
sum S[i+j]*T[j] (0<=j<|T|)がT中の1の個数に等しければOKです.

ここでTをひっくり返したものをT'とすると, T'[|T|-1-j] = T[j]より
sum S[i+j]*T[j]=S[i+j]*T'[|T|-1-j] (0<=j<|T|)となり, これは添字の和が一定なので畳込みになっている
これはFFTで高速に計算出来る(非再帰 FFTの中身は相変わらず知りません...)

計算量は結局0(NlogN) (Nは文字列長さ)

コード:

#include <bits/stdc++.h>
using namespace std;

typedef pair<int, int> pii;
typedef long long ll;
typedef vector<int> vi;

#define pb push_back
#define eb emplace_back
#define mp make_pair
#define fi first
#define se second
#define rep(i,n) rep2(i,0,n)
#define rep2(i,m,n) for(int i=m;i<(n);i++)
#define ALL(c) (c).begin(),(c).end()
 
typedef double Num;
const Num PI = acos(-1.0);
typedef complex<Num> Complex;

void fft_main(int n, Num theta, Complex a[])
{
    for (int m = n; m >= 2; m >>= 1) {
        int mh = m >> 1;
        Complex thetaI = Complex(0, theta);
        rep(i, mh) {
            Complex w = exp((Num)i * thetaI);
            for (int j = i; j < n; j += m) {
                int k = j + mh;
                Complex x = a[j] - a[k];
                a[j] += a[k];
                a[k] = w * x;
            }
        }
        theta *= 2;
    }

    int i = 0;
    for (int j = 1; j < n - 1; ++j) {
        for (int k = n >> 1; k > (i ^= k); k >>= 1) ;
        if (j < i) swap(a[i], a[j]);
    }
}
 
void fft(int n, Complex a[]) { fft_main(n, 2 * PI / n, a); }
void inverse_fft(int n, Complex a[]) { fft_main(n, -2 * PI / n, a); }
 
void convolution(vector<Complex> &v, vector<Complex> &w)
{
    int n = 1, vwn = v.size() + w.size() - 1;
    while (n < vwn) n <<= 1;
    v.resize(n), w.resize(n);
    fft(n, &v[0]);
    fft(n, &w[0]);
    rep(i, n) v[i] *= w[i];
    inverse_fft(n, &v[0]);
    rep(i, n) v[i] /= n;
}
 
string S, T;
int w;
int l1, l2;
int imos[200010];
int ok[200010];

const string gen = "AGCT";

int main() {
    cin >> l1 >> l2 >> w;
    cin >> S >> T;
    reverse(ALL(T));

    for (char c : gen) {
	memset(imos, 0, sizeof(imos));

	rep(i, l1) {
	    if (S[i] == c) {
		int l = max(i - w, 0);
		int r = min(i + w + 1, l1);
		++imos[l];
		--imos[r];
	    }
	}

	vector<Complex> v1, v2;

	rep(i, l1) {
	    imos[i + 1] += imos[i];
	    int f = !!imos[i];
	    v1.pb(f);
	}

	int num = 0;

	rep(i, l2) {
	    int f = T[i] == c;
	    v2.pb(f);
	    num += f;
	}
	convolution(v1, v2);

	for (int i = 0; i <= l1 - l2; ++i) {
	    ok[i] += (int)(v1[i + l2 - 1].real() + 0.5) == num;
	}
    }

    int ret = 0;

    for (int i = 0; i <= l1 - l2; ++i) {
	ret += (ok[i] == 4);
    }

    cout << ret << endl;

    return 0;
}

2016-01-10

SRM456 Div1 hard FunctionalEquation

Topcoder

rngさんの問題表を適当に解いていたら聞いたことあるのが出てきたので解いた. 自力で解けて嬉しい

というか実質初めて解いたhardがこんな重いのになってしまった(平均的難易度を知らないのでなんとも言えない, hardは怖い)

関数方程式という問題名, 一体競技プログラミングとは何だったのか

解説はeditorialがあるはずなので思考過程だけ(正しいものを思いついた順に, 数オリとかやってる/やってた人はもっと賢く解けるのかもしれない)
数式書けないので許してください, メモなので必要十分も適当.

以下ネタバレ

f:Z->Zに注意

まず自由度を減らしたいので色々代入してみる.
f(2f(x) - x + 1) = f(x) + Cのxに2f(x) - x + 1を代入する, 左の式の結果を用いてf(x + 2C) = f(x) + 2Cがわかる. よってx mod 2Cで値が決まる.

少し考えやすくなったが条件が足りなそう. (例えばfが1次関数だと仮定するとf(x) = x + (C - 1) / 2しかだめっぽいし(Cが奇数じゃないとこれもだめ))

試しにf(0) ~ f(2C-1)を文字(a0, a1, ... , a2C-1)で置いてみると, f(x) = ak + (x - k) (k = x mod 2C) = x + (ak - k)
とかける(ak - k = bkと置き直す)

これを与式に代入してみると, f(x + 2bk + 1) = x + bk + C (k = x mod 2C)となる

これx + 2bk + 1が違うx mod 2Cで一緒になったら辛くない??という気分になって

x + 2bk + 1 = y + 2bl + 1 (mod 2C)としてみる. この時 x + bk = y + bl mod 2C も成り立ち, bk = bl, x = yがわかる.

よってk -> k + 1 + 2bk (k=0, ..., 2C-1)の右側もmod 2Cで互いに異なるっぽい.

またここで, l = k + 1 + 2bk mod 2Cとすると,
f(x + 1 + 2bk) = x + bk + Cは, x + 1 + 2bk + bl = x + bk + Cとなって, bk + bl = C - 1 (これは普通の意味での等号)
(x + 1 + 2bk) + 1 + 2bl = x + 2bk + 2 + 2 * (C-1-k) = x+2C (= x mod 2C)で戻ってきた. 思えばx ---- f(x) ---- 2f(x) - xなのであった.

よって0, 2, ..., 2Cと1, 3, ..., 2C - 1で2部グラフと見なせるっぽい.
もうこの対応関係があれば, fは勝手に成り立ちそう.

コストもx[i] mod 2Cで別に考えれば良くて, modを決めればコストが最適化出来るならもうbit dpで計算できる.

ここまで来てコスト計算に苦労してしまった(こっちの方が出来るべきでは...)

2m と2n + 1(0 <= m, n < C)がセットとして, x[i] = 2m or 2n+1 mod 2Cの部分だけ考える.
b2m = tとおく, これらが対応することから2m+2t+1=2n+1 mod 2C つまり2(n-m)=2t mod 2C, n-m=t mod C
(つまりtはn-mからCごとの値を取る)

b2n+1 = C-1-tで,
また|f(x)-y| = |bk-(y-x)| (k = x mod 2C)だから2種類のkについて|bk-t|みたいなのをたくさん足した形になってる

これは三分探索でOK

(editorialに載ってるコスト計算頭良すぎでは)(初期値でsortして絶対値の中身の正負の分かれ目を上手いことする)

2015-09-04

AOJ2636 Distance Sum

問題:
Distance Sum | Aizu Online Judge

ブログ書くのめんどくさい.
問題文はとても読みやすいので読んでください.

解説:

普通に順番に1からNまで追加してシミュレーションします.

まず, i番目を追加した時, i-1番目の時和を最小化する頂点(vとする)とのパス上に新しい候補点がありそう.

vからそっち方向に辺を1つ移動した点に変えた時のcostの変化を考えると, それよりこっち側の部分木の頂点数と向こう側の頂点数(今まで追加された分)(それぞれnum1, num2とする)とcost * (num2 - num1 + 1)減る(1は新しい頂点の分) costは全体にかかるので難しい計算はしなくてよい. i-1番目まででの最小性より num2 - num1 <= 0になっている. num1 + num2 = i - 1 また, 今までの点かlcaっぽいところでnum2 - num1が減少し, 1回変化すると最小となる.

この辺の操作は行き掛け順で管理してsegtreeとかBITとかでサポートするとなんとかなる. (こっちは慣れっぽい)