libicpc

チーム kkntkr / Unknown による、ACM-ICPC 向けのアルゴリズムの実装をまとめたページです。

基礎

テンプレート/マクロ

ループなどを簡略して記述するためのマクロなど。バッドノウハウっぽいけど便利だし…

1
2
3

#define REP(i,n) for(int i = 0; i < (int)(n); i++)
#define FOR(it,c) for(__typeof((c).begin()) it = (c).begin(); it != (c).end(); ++it)
#define ALLOF(c) (c).begin(), (c).end()

計算/ビット演算

各種ビット演算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

// 1であるビットの数を数える。
unsigned int __builtin_popcount(unsigned int);
unsigned long long __builtin_popcountll(unsigned long long);

// 末尾の0であるビットの数を数える。
unsigned int __builtin_ctz(unsigned int);
unsigned long long __builtin_ctzll(unsigned long long);

// 最右ビットを抜き出す。
y = x&-x;

// 最右ビットを落とす。
y = x&(x-1);

// 非零部分集合を列挙する。
#define FOR_SUBSET(b, a) for(int b = (a)&-(a); b != 0; b = (((b|~(a))+1)&(a)))
#define FOR_SUBSET(b, a) for(int b = a; b != 0; b = (b-1)&a)

// ビットコンビネーションを列挙する。
int next_combination(int p) {
    int lsb = p&-p;
    int rem = p+lsb;
    int rit = rem&~p;
    return rem|(((rit/lsb)>>1)-1);
}

計算/実数比較

微小な誤差を無視した実数の比較を行う。

1
2
3
4
5

// 実数の比較誤差
#define EPS 1.0e-10

// 実数の符号関数
inline int signum(double x) { return (abs(x) < EPS ? 0 : x > 0 ? 1 : -1); }

幾何

基礎/データ構造

二次元幾何図形のデータ構造。点を complex<double> で表す。

1
2
3
4
5
6
7
8
9
10
11

// 点
typedef complex<double> P;

// 線分・半直線・直線
struct L { P pos, dir; };

// 多角形
typedef vector<P> G;

// 円
struct C { P p; double r; };

検証なし

基礎/内積・外積

ベクトルの内積と外積。他の関数から多用される。

1
2
3
4
5
6
7
8
9

// 二つのベクトルの内積を計算する
inline double inp(const P& a, const P& b) {
    return (conj(a)*b).real();
}

// 二つのベクトルの外積を計算する
inline double outp(const P& a, const P& b) {
    return (conj(a)*b).imag();
}

検証済み

基礎/回転方向関数

線分などの交差判定に用いる回転方向関数。基本的には signum(outp(a, b)) だが、outp(a, b) == 0 のときにa, bの配置によっておせっかいな挙動をする。

1
2
3
4
5
6
7
8
9
10
11

inline int ccw(const P& p, const P& r, const P& s) {
    P a(r-p), b(s-p);
    int sgn = signum(outp(a, b));
    if (sgn != 0)
        return sgn;
    if (a.real()*b.real() < -EPS || a.imag()*b.imag() < -EPS)
        return -1;
    if (norm(a) < norm(b) - EPS)
        return 1;
    return 0;
}

検証済み

  • UVA10709 Intersection is not that easy

基礎/射影

ベクトルや線分などを直線上に射影する。射影を用いると簡単に記述できるコードがけっこうある。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

// ベクトルpをベクトルbに射影したベクトルを計算する
inline P proj(const P& p, const P& b) {
    return b*inp(p,b)/norm(b);
}

// 点pから直線lに引いた垂線の足となる点を計算する
inline P perf(const L& l, const P& p) {
    L m = {l.pos - p, l.dir};
    return (p + (m.pos - proj(m.pos, m.dir)));
}

// 線分sを直線bに射影した線分を計算する
inline L proj(const L& s, const L& b) {
     return (L){perf(b, s.pos), proj(s.dir, b.dir)};
}

未検証

面積・体積/円と円の共通部分

円と円の共通部分の面積を求める。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

double cc_area(const C& c1, const C& c2) {
    double d = abs(c1.p - c2.p);
    if (c1.r + c2.r <= d + EPS) {
        return 0.0;
    } else if (d <= abs(c1.r - c2.r) + EPS) {
        double r = c1.r <? c2.r;
        return r * r * PI;
    } else {
        double rc = (d*d + c1.r*c1.r - c2.r*c2.r) / (2*d);
        double theta = acos(rc / c1.r);
        double phi = acos((d - rc) / c2.r);
        return c1.r*c1.r*theta + c2.r*c2.r*phi - d*c1.r*sin(theta);
    }
}

検証済み

  • PKU2546 Circular Area

面積・体積/多角形の面積

(凸とは限らない)多角形の面積を求める。

1
2
3
4
5
6
7
8
9

double area(G& g) {
    int n = g.size();
    double s = 0.0;
    for(int i = 0; i < n; i++) {
        int j = (i+1)%n;
        s += outp(g[i], g[j])/2;
    }
    return abs(s);
}

未検証

交差/円と円の交点

円と円の交点を計算する。

前提条件

  • 交点が0個または無限個だとうまく動かない
1
2
3
4
5
6
7

pair<P, P> cc_cross(const C& c1, const C& c2) {
    double d = abs(c1.p - c2.p);
    double rc = (d*d + c1.r*c1.r - c2.r*c2.r) / (2*d);
    double rs = sqrt(c1.r*c1.r - rc*rc);
    P diff = (c2.p - c1.p) / d;
    return make_pair(c1.p + diff * P(rc, rs), c1.p + diff * P(rc, -rs));
}

未検証

交差/円と直線の交差判定

円と直線が交差するか判定する。

1
2
3

bool cl_intersects(const C& c, const L& l) {
    return lp_distance(l, c.p) <= c.r + EPS;
}

未検証

交差/円と直線の交点

円と直線の交点を計算する。

1
2
3
4
5
6
7
8
9
10
11
12
13

vector<P> cl_cross(const C& c, const L& l) {
  P h = perf(l, c.p);
  double d = abs(h - c.p);
  vector<P> res;
  if(d < c.r - EPS){
    P x = l.dir / abs(l.dir) * sqrt(c.r*c.r - d*d);
    res.push_back(h + x);
    res.push_back(h - x);
  }else if(d < c.r + EPS){
    res.push_back(h);
  }
  return res;
}

未検証

交差/凸多角形と線分の包含判定

凸多角形と線分が内部で交差しているか判定する。 接しているだけの場合には交差しているとみなさない。

前提条件

  • 多角形は凸でなければいけない。
  • 多角形は反時計回りに与えられなければならない。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

// *** BE CAREFUL OF THE PRECONDITIONS! ***
bool cgs_contains_exclusive(G& r, P from, P to) {
    int m = r.size();

    P dir0 = to - from;

    double lower = 0.0, upper = 1.0;
    REP(i, m) {
        int j = (i+1)%m;
        P a = r[i], b = r[j];
        P ofs = a - from;
        P dir1 = b - a;
        double num = (conj(ofs)*dir1).imag();
        double denom = (conj(dir0)*dir1).imag();
        if (abs(denom) < EPS) {
            if (num < EPS)
                lower = upper = 0.0;
        }
        else {
            double t = num / denom;
            if (denom > 0)
                upper <?= t;
            else
                lower >?= t;
        }
    }

    return (upper-lower > EPS);
}

検証済み

  • PKU1774 Fold Paper Strips

交差/多角形と点の包含判定

(凸とは限らない)多角形の内部に点が包含されているか判定する。

1
2
3
4
5
6
7
8
9
10
11

bool gp_contains(const G& g, const P& p) {
    double sum = 0.0;
    int n = g.size();
    for(int i = 0; i < n; i++) {
        int j = (i+1)%n;
        if (sp_intersects(L(g[i], g[j]-g[i]), p))
            return true;
        sum += arg((g[j]-p)/(g[i]-p));
    }
    return (abs(sum) > 1);
}

未検証

交差/直線と直線の交差判定

1
2
3

bool ll_intersects(const L& l, const L& m) {
    return (abs(outp(l.dir, m.dir)) > EPS || abs(outp(l.dir, m.pos-l.pos)) < EPS);
}

未検証

  • UVA10709 Intersection is not that easy

交差/直線と直線の交点

直線同士の交点を計算する。

1
2
3
4
5

P line_cross(const L& l, const L& m) {
    double num = outp(m.dir, m.pos-l.pos);
    double denom = outp(m.dir, l.dir);
    return P(l.pos + l.dir*num/denom);
}

未検証

メモ

  • denom = 0 となるケース、すなわち 直線が平行となるとき に注意すること。
  • 片方または両方が線分のときに使用する場合は、{ls,ss}_intersects()等と併用して、交点がvalidなものか確認すること。

交差/直線と線分の交差判定

直線と線分が交差するか判定する。

1
2
3
4

bool ls_intersects(const L& l, const L& s) {
    return (signum(outp(l.dir, s.pos-l.pos)) *
            signum(outp(l.dir, s.pos+s.dir-l.pos)) <= 0);
}

検証済み

  • UVA10709 Intersection is not that easy

交差/線分と点の交差判定

線分と点が交差するか判定する。

1
2
3

bool sp_intersects(const L& s, const P& p) {
    return ( abs(s.pos - p) + abs(s.pos + s.dir - p) - abs(s.dir) < EPS );
}

未検証

交差/線分と線分の交差判定

線分と線分が交差するか判定する。

1
2
3
4
5
6

bool ss_intersects(const L& s, const L& t) {
    return (ccw(s.pos, s.pos+s.dir, t.pos) *
            ccw(s.pos, s.pos+s.dir, t.pos+t.dir) <= 0 &&
            ccw(t.pos, t.pos+t.dir, s.pos) *
            ccw(t.pos, t.pos+t.dir, s.pos+s.dir) <= 0);
}

検証済み

  • UVA10709 Intersection is not that easy

距離/最遠点対

凸包の直径を求める。キャリパー法。任意の点集合に対し、凸包を求めてから直径を求めることで最遠点対を求めることができる。

前提条件

  • 点集合は凸包であり、反時計回りに与えられること。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

// BE CAREFUL OF PRECONDITIONS!!!
double farthest(const vector<P>& v) {
    int n = v.size();

    int si = 0, sj = 0;
    REP(t, n) {
        if (v[t].imag() < v[si].imag())
            si = t;
        if (v[t].imag() > v[sj].imag())
            sj = t;
    }

    double res = 0;
    int i = si, j = sj;
    do {
        res >?= abs(v[i]-v[j]);
        P di = v[(i+1)%n] - v[i];
        P dj = v[(j+1)%n] - v[j];
        if (outp(di, dj) >= 0)
            j = (j+1)%n;
        else
            i = (i+1)%n;
    } while(!(i == si && j == sj));

    return res;
}

検証済み

  • PKU2187 Beauty Contest

メモ

  • O(n).

距離/直線と点の距離

1
2
3

double lp_distance(const L& l, const P& p) {
    return abs(outp(l.dir, p-l.pos) / abs(l.dir));
}

検証済み

  • UVA10709 Intersection is not that easy

距離/直線と直線の距離

1
2
3

double ll_distance(const L& l, const L& m) {
    return (ll_intersects(l, m) ? 0 : lp_distance(l, m.pos));
}

検証済み

  • UVA10709 Intersection is not that easy

距離/直線と線分の距離

1
2
3
4
5

double ls_distance(const L& l, const L& s) {
    if (ls_intersects(l, s))
        return 0;
    return min(lp_distance(l, s.pos), lp_distance(l, s.pos+s.dir));
}

検証済み

  • UVA10709 Intersection is not that easy

距離/線分と点の距離

線分と点の距離を計算する。

1
2
3
4
5
6
7
8

double sp_distance(const L& s, const P& p) {
    const P r = perf(s, p);
    const double pos = ((r-s.pos)/s.dir).real();
    if (-EPS <= pos && pos <= 1 + EPS)
        return abs(r - p);
    return min(abs(s.pos - p),
               abs(s.pos+s.dir - p));
}

検証済み

  • UVA10709 Intersection is not that easy

距離/線分と線分の距離

線分同士の距離を計算する。

1
2
3
4
5
6
7
8

double ss_distance(const L& s, const L& t) {
    if (ss_intersects(s, t))
        return 0;
    return (sp_distance(s, t.pos) <?
            sp_distance(s, t.pos+t.dir) <?
            sp_distance(t, s.pos) <?
            sp_distance(t, s.pos+s.dir));
}

検証済み

  • UVA10709 Intersection is not that easy

多角形/凸包

点をy,xの昇順でソートする。この順序通りに反時計回りの凸な折れ線を求めると、下辺＋右側＋上辺の最右点になる。逆順序で同じものを求め、それぞれの最後の点を取り除いて接続すると、凸包になっている。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

bool xy_less(const P& a, const P& b) {
    if (abs(a.imag()-b.imag()) < EPS) return (a.real() < b.real());
    return (a.imag() < b.imag());
}
double ccw(P a, P b, P c) {
    return (conj(b-a)*(c-a)).imag();
}
template<class IN>
void walk_rightside(IN begin, IN end, vector<P>& v) {
    IN cur = begin;
    v.push_back(*cur++);
    vector<P>::size_type s = v.size();
    v.push_back(*cur++);
    while(cur != end) {
        if (v.size() == s || ccw(v[v.size()-2], v.back(), *cur) > EPS)
            v.push_back(*cur++);
        else
            v.pop_back();
    }
    v.pop_back();
}
vector<P> convex_hull(vector<P> v) {
    if (v.size() <= 1)
        return v; // EXCEPTIONAL
    sort(ALLOF(v), xy_less);
    vector<P> cv;
    walk_rightside(v.begin(), v.end(), cv);
    walk_rightside(v.rbegin(), v.rend(), cv);
    return cv;
}

未検証

  • UVA11168 Airport
  • CII3169 Boundary Points

メモ

  • O(n log n)。
  • コードは最小の点数からなる凸包を求める場合。凸包の辺上にある点を全て取りたい場合は ccwの比較を-EPSに変更する。
  • 最大点数を取るように変更した場合、一直線に並んだ点群を処理すると凸包が元の点数より大きくなることに注意。
  • 点がひとつのときに注意。

多角形/凸多角形のクリッピング

凸多角形を直線で切り取る。

前提条件

  • 多角形はcounterclockwiseに与えられていること。
  • 直線の右側が切り取られ、左側が残される。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

G cut_polygon(const G& g, L cut) {
    int n = g.size();

    G res;

    REP(i, n) {
        P from(g[i]), to(g[(i+1)%n]);
        double p1 = (conj(cut.dir)*(from-cut.pos)).imag();
        double p2 = (conj(cut.dir)*(to-cut.pos)).imag();
        if (p1 > -EPS) {
            res.push_back(from);
            if (p2 < -EPS && p1 > EPS)
                res.push_back(from - (to-from)*p1/(conj(cut.dir)*(to-from)).imag());
        }
        else if (p2 > EPS) {
            res.push_back(from - (to-from)*p1/(conj(cut.dir)*(to-from)).imag());
        }
    }

    return res;
}

検証済み

  • UVA10117 Nice Milk
  • PKU1755 Triathlon

その他/アレンジメント

与えられた直線の集合に対し、直線同士の交点からなるグラフを作成する。

線分から作成したいときはlp_intersectsをsp_intersectsに置き換える。ただし、端点もノードにしたい場合は適宜追加すること。

前提条件

  • 重なる直線が存在してはいけない
  • lp_intersectsはNaNである点に対しても正しく動くこと
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

// BE CAREFUL OF NaN POINTS & DUPLICATE LINES
typedef map<P, vector<P> > PGraph;
PGraph arrange(const vector<L>& lines) {
    int m = lines.size();

    set<P> points;
    REP(j, m) REP(i, j) {
        L a = lines[i], b = lines[j];
        double r = outp(b.dir, b.pos-a.pos) / outp(b.dir, a.dir);
        P p = a.pos + a.dir * r;
        if (lp_intersects(a, p) && lp_intersects(b, p))
            points.insert(p);
    }

    PGraph g;
    REP(k, m) {
        vector< pair<double, P> > onlines;
        FOR(it, points) if (lp_intersects(lines[k], *it))
            onlines.push_back(make_pair(inp(*it-lines[k].pos, lines[k].dir), *it));
        sort(ALLOF(onlines));
        REP(i, (int)onlines.size()-1) {
            g[onlines[i+0].second].push_back(onlines[i+1].second);
            g[onlines[i+1].second].push_back(onlines[i+0].second);
        }
    }

    return g;
}

未検証

メモ

  • O(n3).

その他/ダイス

ダイスを回転させる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

struct die_t {
    int top, front, right, left, back, bottom;
};

#define rotate_swap(x,a,b,c,d) swap(x.a, x.b); swap(x.b, x.c); swap(x.c, x.d);
void rotate_right(die_t& x) { rotate_swap(x,top,left,bottom,right); }
void rotate_left(die_t& x) { rotate_swap(x,top,right,bottom,left); }
void rotate_front(die_t& x) { rotate_swap(x,top,back,bottom,front); }
void rotate_back(die_t& x) { rotate_swap(x,top,front,bottom,back); }
void rotate_cw(die_t& x) { rotate_swap(x,back,left,front,right); }
void rotate_ccw(die_t& x) { rotate_swap(x,back,right,front,left); }

void generate_all() {
    REP(i, 6) {
        REP(j, 4) {
            emit();
            rotate_cw(x);
        }
        (i%2 == 0 ? rotate_right : rotate_front)(x);
    }
}

三次元幾何/直線と直線の距離

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

double ll_distance(L& a, L& b) {
    P pa = a.from;
    P pb = b.from;
    P da = a.to - a.from;
    P db = b.to - b.from;
    double num = inp( (db * inp(da, db) - da * inp(db, db)) , pb - pa );
    double denom = inp(da, db) * inp(da, db) - inp(da, da) * inp(db, db);
    double ta;
    if (abs(denom) < EPS)
        ta = 0;
    else
        ta = num / denom;
    P p = pa + da * ta;
    return lp_distance(b, p);
}

検証済み

  • 2007夏合宿Day1C Save the Energy

グラフ

基礎/データ構造

隣接リストによるグラフの表現。重みにはWeight型を用いて、無限大の重みをWEIGHT_INFTYで参照する。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

// 重み
typedef double Weight;
const Weight WEIGHT_INFTY = numeric_limits<Weight>::max() / 4;

// 枝
struct Edge {
    int src, dest;
    Weight weight;
};
bool operator<(const Edge& a, const Edge& b) {
    return (a.weight < b.weight);
}

// グラフ(隣接リスト)
typedef vector<Edge> Edges;
typedef vector<Edges> Graph;

検証なし

最短路/Bellman-Ford

負の重みを持つ枝が存在するときの一対全最短経路を求めるアルゴリズム。negative cycleの存在を確認するためにも用いられる。

枝の重みが非負であるなら、計算量的により優れたDijkstra法を用いるほうがいい。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

Graph g;

bool shortest(int start) {
    int n = g.size();
    vector<Weight> costs(n, WEIGHT_INFTY);
    vector<int> trace(n, -2);
    costs[start] = 0;
    trace[start] = -1;    // TERMINAL

    REP(k, n) REP(a, n) if (costs[a] != WEIGHT_INFTY) FOR(it, g[a]) {
        int b = it->dest;
        Weight w = costs[a] + it->weight;
        if (w < costs[b]) {
            costs[b] = w;
            trace[b] = a;
        }
    }
    REP(a, n) FOR(it, g[a]) // negative cycleのチェック
        if (costs[a] + it->weight < costs[it->dest])
            return false;
    return true;
}

未検証

最短路/Dijkstra

枝の重みが非負であるグラフにおいて、一対全最短経路を求めるアルゴリズム。

負の重みを持つ枝が存在するならBellman-Ford法などのアルゴリズムを用いる。

前提条件

  • パスは非負の重みであること。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

Graph g;
vector<int> trace;

Weight shortest(int start, int goal) {
    int n = g.size();
    trace.assign(n, -2); // UNREACHABLE

    priority_queue<Edge, vector<Edge>, greater<Edge> > q;
    q.push((Edge){-1, start, 0}); // TERMINAL

    while(!q.empty()) {
        Edge e = q.top();
        q.pop();
        if (trace[e.dest] >= 0)
            continue;
        trace[e.dest] = e.src;
        if (e.dest == goal)
            return e.weight;
        FOR(it, g[e.dest])
            if (trace[it->dest] == -2)
                q.push((Edge){it->src, it->dest, e.weight + it->weight});
    }
    return -1;
}

vector<int> buildRoute(int goal) {
    if (trace[goal] == -2)
        return vector<int>(); // UNREACHABLE
    vector<int> route;
    for(int i = goal; i >= 0; i = trace[i])
        route.push_back(i);
    reverse(route.begin(), route.end());
    return route;
}

検証済み

  • UVA10525 NEW TO BANGLADESH?

最短路/Warshall-Floyd

全対全最短路を求めるアルゴリズム。

経路長のみを知りたい場合は数行のループで書けてしまう。ノード数が少ないときに最短路を求めたい場合はこのアルゴリズムを選択すると労力が少なくて良い。

前提条件

  • 自身への距離を0とするかどうかに応じて隣接行列の対角成分を初期化すること。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

Weight g[N][N];
int trace[N][N];

void shortest() {
    REP(i, n) REP(j, n)
        trace[i][j] = -1;
    REP(j, n) REP(i, n) REP(k, n) {
        if (g[i][k] > g[i][j] + g[j][k]) {
            g[i][k] = g[i][j] + g[j][k];
            trace[i][k] = j;
        }
    }
}
void buildRoute(vector<int>& route, int src, int dest, bool rec = false) {
    if (!rec)
        route.clear();
    int inter = trace[src][dest];
    if (inter < 0) {
        route.push_back(src);
    }
    else {
        buildRoute(route, src, inter, true);
        buildRoute(route, inter, dest, true);
    }
    if (!rec)
        route.push_back(dest);
}

検証済み

  • CII3569 World Finals 2005 Degrees of Separation

最小広域木/Kruskal

最小広域木を求めるアルゴリズム。同一連結成分に属しているかの判定を行うためUnion-Findアルゴリズムを必要とする。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

pair<Weight, Edges> kruskal(Graph& g) {
    int n = g.size();
    vector<Edge> edges;
    REP(a, n) FOR(it, g[a])
        if (a < it->dest)
            edges.push_back(*it);
    sort(ALLOF(edges));
    UnionFind uf(n);
    Edges tree;
    Weight weight = 0;
    REP(i, edges.size()) {
        if ((int)tree.size() >= n-1)
            break;
        Edge& e = edges[i];
        if (uf.link(e.src, e.dest)) {
            tree.push_back(e);
            weight += e.weight;
        }
    }
    return make_pair(weight, tree);
}

未検証

  • UVA10600 ACM CONTEST AND BLACKOUT

最小広域木/Prim

最小広域木を求めるアルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

// 戻り値は(木の重み, 木)
pair<Weight, Edges> prim(Graph& g) {
    int n = g.size();
    priority_queue<Edge, vector<Edge>, greater<Edge> > q;
    vector<bool> visited(n, false);
    Edges tree;
    Weight weight = 0;

    q.push((Edge){-1, 0, 0});
    while(!q.empty()) {
        Edge e = q.top();
        q.pop();
        if (visited[e.dest])
            continue;
        visited[e.dest] = true;
        if (e.src >= 0) {
            tree.push_back(e);
            weight += e.weight;
        }
        FOR(it, g[e.dest])
            if (!visited[it->dest])
                q.push(*it);
    }
    return make_pair(weight, tree);
}

検証済み

  • UVA10034 Freckles
  • PKU1251 Jungle Roads

メモ

  • 隣接行列表現を用いたバージョンはO(V3)。

マッチング/Hopcroft-Karp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

vector< vector<int> > g;
int n, m;

vector<bool> visited, matched;
vector<int> levels, matching;

bool augment(int left) {
    if (left == n)
        return true;
    if (visited[left])
        return false;
    visited[left] = true;
    REP(i, g[left].size()) {
        int right = g[left][i];
        int next = matching[right];
        if (levels[next] > levels[left] && augment(next)) {
            matching[right] = left;
            return true;
        }
    }
    return false;
}

int match() {
    matching.assign(m, n);
    matched.assign(n, false);
    bool cont;
    do {
        levels.assign(n+1, -1);
        levels[n] = n;
        queue<int> q;
        REP(left, n) {
            if (!matched[left]) {
                q.push(left);
                levels[left] = 0;
            }
        }
        while(!q.empty()) {
            int left = q.front();
            q.pop();
            REP(i, g[left].size()) {
                int right = g[left][i];
                int next = matching[right];
                if (levels[next] < 0) {
                    levels[next] = levels[left] + 1;
                    q.push(next);
                }
            }
        }
        visited.assign(n, false);
        cont = false;
        REP(left, n)
            if (!matched[left] && augment(left))
                matched[left] = cont = true;
    } while(cont);
    return count(ALLOF(matched), true);
}

検証済み

  • TopCoder TCCC07 Semi1 Hard AncientLanguage

メモ

  • O(sqrt(V)(V+E)).

マッチング/ハンガリアン法

Hungarian methodのO(n3)実装。

重み付きマッチングにおいて、重み行列の各行各列に定数を加えることはマッチングの結果に影響を及ぼさない。

この性質を用いて、まず全ての行について、その行の最小値を同じ行の全ての要素にから引くことを行う。前述の通り、これを行ってもマッチングの重みにオフセットを加えただけであることに注意。すると、全ての重みが非負になるので、もし重み0で完全マッチングを作る事ができれば、それは最小重みマッチングになる。

さて、重み0のマッチングを作るために、重み0の枝からなる辺誘導二部グラフ上でマッチングを構成することを考える。重み無しマッチングの要領で交代路を求めるのだが、誘導グラフ上で最大マッチングが存在するとは限らない。そこで、交代路を伸ばす過程でそれ以上伸ばす事ができなくなったとき、適切に行・列ごとの重みの調整を行うことで重み0の枝を増やす。それが、コードの16～21行目の処理である。

(実際にはコードを簡潔にするために処理に手を加えてある)

前提条件

  • n <= m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

#define residue(i,j) (adj[i][j] + ofsleft[i] + ofsright[j])

vector<int> min_cost_match(vector< vector<int> > adj) {
    int n = adj.size();
    int m = adj[0].size();
    vector<int> toright(n, -1), toleft(m, -1);
    vector<int> ofsleft(n, 0), ofsright(m, 0);

    REP(r, n) {
        vector<bool> left(n, false), right(m, false);
        vector<int> trace(m, -1), ptr(m, r);
        left[r] = true;

        for(;;) {
            int d = numeric_limits<int>::max();
            REP(j, m) if (!right[j])
                d <?= residue(ptr[j], j);
            REP(i, n) if (left[i])
                ofsleft[i] -= d;
            REP(j, m) if (right[j])
                ofsright[j] += d;
            int b = -1;
            REP(j, m) if (!right[j] && residue(ptr[j], j) == 0)
                b = j;
            trace[b] = ptr[b];
            int c = toleft[b];
            if (c < 0) {
                while(b >= 0) {
                    int a = trace[b];
                    int z = toright[a];
                    toleft[b] = a;
                    toright[a] = b;
                    b = z;
                }
                break;
            }
            right[b] = left[c] = true;
            REP(j, m) if (residue(c, j) < residue(ptr[j], j))
                ptr[j] = c;
        }
    }
    return toright;
}

検証済み

  • UVA10888 Warehouse
  • PKU2195 Going Home

メモ

  • O(n2 m)。

マッチング/二部グラフのマッチング

二部グラフの最大マッチングを求めるアルゴリズム。基本的には最大流の特殊化だが、0/1条件からEdmonds-KarpのようにBFSを用いずとも多項式計算量が保証できる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

vector< vector<int> > g;
int n, m;

vector<bool> visited;
vector<int> matching;

bool augment(int left) {
    if (left < 0)
        return true;
    if (visited[left])
        return false;
    visited[left] = true;
    REP(i, g[left].size()) {
        int right = g[left][i];
        if (augment(matching[right])) {
            matching[right] = left;
            return true;
        }
    }
    return false;
}
int match() {
    matching.assign(m, -1);
    int matches = 0;
    REP(left, n) {
        visited.assign(n, false);
        if (augment(left))
            matches++;
    }
    return matches;
}

検証済み

  • UVA10080 GopherII
  • TopCoder TCCC07 Semi1 Hard AncientLanguage

メモ

  • O(V(V+E)).
  • 二部グラフにおいて、最大マッチングの大きさは最小点被覆の大きさに一致する。

マッチング/安定結婚

N人の男性とN人の女性について、各人の希望リストが与えられたときに、安定なマッチングを1つ求める。ある人の希望リストとは、その人の好みで異性に全順序をつけたもの。安定なマッチングとは、現在組んでいるペア(a, b)と(a’, b’)でaがbよりもb’を好みb’がa’よりもaを好むような対(a, b’)(ブロッキングペア)が存在しないマッチング。任意の希望リストについて、安定なマッチングが少なくとも1つ存在することが知られている。

実装はGale-Shapleyのアルゴリズム。 男性優先の解を求める。 返り値は res[womanID] = manID なる割り当て。

前提条件

  • 各人について、優先度の高い相手から順に並んでいること。
  • _man[manID][rank] = womanID, _woman[womanID][rank] = manID
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

inline vector<int> stable_marriage(const vector<vector<int> >& _man,
                 const vector<vector<int> >& _woman) {
  int n = _man.size();
  vector<vector<int> > man(n, n), woman(n, n);
  vector<int> res(n, -1);

  REP(i, n) REP(j, n){
    man[i][j] = _man[i][n-j-1];
    woman[i][_woman[i][j]] = n - j - 1;
  }

  REP(i, n){
    for(int manID = i; manID >= 0; ){
      int womanID = man[manID].back();
      man[manID].pop_back();
      int prevman = res[womanID];
      if(prevman < 0 || woman[womanID][prevman] < woman[womanID][manID]){
  res[womanID] = manID;
  manID = prevman;
      }
    }
  }

  return res;
}

検証済み

  • Seoul 2000 Match Maker(夏合宿2004の2日目問題F)

メモ

  • O(n2)。

最大流/Dinic

二部グラフマッチングのHopcroft-Karpのアルゴリズムとほぼ同様のアイデアによるもの。

ソースからの最短距離(レベル)に基づき、補助ネットワーク内に最短路のみからなるレベルグラフを構成し、その上で極大フローを流す。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

int n, s, t;
int cap[N][N], flow[N][N];
int levels[N];
bool finished[N];
#define residue(i,j) (cap[i][j] - flow[i][j])

void levelize() {
    memset(levels, -1, sizeof(levels));
    queue<int> q;
    levels[s] = 0; q.push(s);
    while(!q.empty()) {
        int here = q.front(); q.pop();
        REP(there, n) if (levels[there] < 0 && residue(here, there) > 0) {
            levels[there] = levels[here] + 1;
            q.push(there);
        }
    }
}
int augment(int here, int cur) {
    if (here == t || cur == 0)
        return cur;
    if (finished[here])
        return 0;
    finished[here] = true;
    REP(there, n) if (levels[there] > levels[here]) {
        int f = augment(there, min(cur, residue(here, there)));
        if (f > 0) {
            flow[here][there] += f;
            flow[there][here] -= f;
            finished[here] = false;
            return f;
        }
    }
    return 0;
}
int maxflow() {
    memset(flow, 0, sizeof(flow));
    int total = 0;
    for(bool cont = true; cont; ) {
        cont = false;
        levelize();
        memset(finished, 0, sizeof(finished));
        for(int f; (f = augment(s, INF)) > 0; cont = true)
            total += f;
    }
    return total;
}

検証済み

  • PKU2391 Ombrophobic Bovines

メモ

  • O(n4).
  • 隣接リスト表現を使えばO(V2E).

最大流/Edmonds-Karp

最大流を求めるEdmonds-Karpのアルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

Weight g[N][N], f[N][N];

#define RESIDUE(from,to) (g[from][to] - f[from][to])
Weight max_flow(int n, int s, int t) {
    REP(i, n) REP(j, n)
        f[i][j] = 0;
    Weight flow = 0;
    for(;;) {
        queue<int> q;
        q.push(s);
        vector<int> trace(n, -1);
        trace[s] = s;
        while(!q.empty() && trace[t] < 0) {
            int i = q.front();
            q.pop();
            REP(j, n) if (trace[j] < 0 && RESIDUE(i, j) > 0) {
                trace[j] = i;
                q.push(j);
            }
        }
        if (trace[t] < 0)
            break;
        Weight w = WEIGHT_INFTY;
        for(int j = t; trace[j] != j; j = trace[j])
            w <?= RESIDUE(trace[j], j);
        for(int j = t; trace[j] != j; j = trace[j]) {
            f[trace[j]][j] += w;
            f[j][trace[j]] -= w;
        }
        flow += w;
    }
    return flow;
}

未検証

  • 旧ライブラリに載っていたものにはバグがある可能性がある(逆方向のフローを修正していない)。
  • PKU1273 Drainage Ditches
  • PKU1459 Power Network

最大流/最小カット

枝に重みが付いた無向グラフを非連結にするために必要な最小のコストを求める。

指定された二つのノードを切り離す場合は最大流を求めてMaxflow-Mincut定理を適用する。

Stoer-Wagnerのアルゴリズム。http://www.cs.dartmouth.edu/~rahul/Teaching/stoerwagner-mincut.pdf

前提条件

  • 非負の重みであること。
  • 隣接行列adjは破壊される。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

Weight adj[N][N];
int n;

int min_cut() {
    Weight res = WEIGHT_INFTY;

    vector<int> v;
    REP(i, n)
        v.push_back(i);

    for(int m = n; m > 1; m--) {
        vector<Weight> ws(m, 0);
        int s, t;
        Weight w;
        REP(k, m) {
            s = t;
            t = max_element(ws.begin(), ws.end()) - ws.begin();
            w = ws[t];
            ws[t] = -1;
            REP(i, m)
                if (ws[i] >= 0)
                    ws[i] += adj[v[t]][v[i]];
        }
        REP(i, m) {
            adj[v[i]][v[s]] += adj[v[i]][v[t]];
            adj[v[s]][v[i]] += adj[v[t]][v[i]];
        }
        v.erase(v.begin()+t);
        res <?= w;
    }
    return res;
}

検証済み

  • PKU2914 Minimum Cut

メモ

  • O(n3)。
  • グラフの隣接リスト表現とpriority queueを使えばO(VElogE)。

最小費用流/Primal-Dual

最小費用最大流を求めるPrimal-Dual法。

前提条件

  • コスト行列は反対称であること。
  • キャパシティ行列は非負であること。
  • キャパシティが正の枝はコストが非負であること。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

// BE CAREFUL OF PRECONDITIONS!!!!
#define residue(i,j) (cap[i][j] - flow[i][j])
#define rcost(i,j) (cost[i][j] + pot[i] - pot[j])
Weight min_cost_flow(vector< vector<Flow> > cap, vector< vector<Weight> > cost, int s, int t) {
    int n = cap.size();

    vector< vector<Flow> > flow(n, n);
    vector<Weight> pot(n);

    Weight res = 0;
    for(;;) {

        vector<Weight> dists(n+1, WEIGHT_INFTY);
        vector<bool> visited(n, false);
        vector<int> trace(n, -1);
        dists[s] = 0;

        for(;;) {
            int cur = n;
            REP(i, n) if (!visited[i] && dists[i] < dists[cur])
                cur = i;
            if (cur == n)
                break;
            visited[cur] = true;
            REP(next, n) if (residue(cur, next) > EPS) {
                if (dists[next] > dists[cur] + rcost(cur, next) && !visited[next]) {
                    dists[next] = dists[cur] + rcost(cur, next);
                    trace[next] = cur;
                }
            }
        }
        if (!visited[t])
            break;

        Flow f = FLOW_INFTY;
        for(int c = t; c != s; c = trace[c])
            f <?= residue(trace[c], c);
        for(int c = t; c != s; c = trace[c]) {
            flow[trace[c]][c] += f;
            flow[c][trace[c]] -= f;
            res += cost[trace[c]][c] * f;
        }

        REP(i, n)
            pot[i] += dists[i];
    }

    return res;
}

検証済み

  • UVA10746 Crime Wave - The Sequel

最小費用流/最小平均長閉路

平均長が最小である単純閉路を探すアルゴリズム。Karp[1]によるBellman-Fordアルゴリズムの変形で、次の事実を用いる:

\mathrm{minavg} = \min_{v \in V} \max_{0 \leq i < n} \frac{d_n(v) - d_i(v)}{n - i}.

このコードで求まるのは閉路の平均長。閉路自体を求めるには、Bellman-Fordの過程で最小値に至るルートを保存して後で辿る。

最小費用循環を効率よく求める際に用いる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

Edges edges;

Weight dp[N+1][N] = { { 0 } };
REP(k, n) {
    REP(i, n)
        dp[k+1][i] = INF;
    REP(i, m) {
        Edge& e = edges[i];
        if (dp[k][e.src] < INF)
            dp[k+1][e.dest] <?= dp[k][e.src] + e.weight;
    }
}
Weight res = INF;
REP(i, n) {
    Weight lo = -INF;
    REP(k, n) if (dp[n][i] < INF && dp[k][i] < INF)
        lo >?= (dp[n][i] - dp[k][i]) / (n - k);
    if (lo > -INF)
        res <?= lo;
}

検証済み

  • PKU2949 Word Rings

メモ

  • O(nm).

参考文献

  • Karp, R. M. A characterization of the minimum cycle mean in a digraph. Discrete Math. 23 (1978), 309–311.

連結成分/Union-Find

オンラインで二つのノードの併合・同一連結成分に属するかどうかの判定を行うアルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

int n;
vector<int> data(n, -1);
bool link(int a, int b) { // 新たな併合を行うとtrue
    int ra = find_root(a);
    int rb = find_root(b);
    if (ra != rb) {
        if (data[rb] < data[ra])
            swap(ra, rb);
        data[ra] += data[rb];
        data[rb] = ra;
    }
    return (ra != rb);
}
bool check(int a, int b) { // 同じ集合ならtrue
    return (find_root(a) == find_root(b));
}
int find_root(int a) { // 代表元を返す
    return ((data[a] < 0) ? a : (data[a] = find_root(data[a])));
}

検証済み

  • UVA10600 ACM CONTEST AND BLACKOUT

メモ

  • -data[find_root(i)]はノードiが属する集合の大きさになっている。

連結成分/二重頂点連結成分

点aを除去したときにm個の連結成分に分割される時、artsにm-1個の数字aが代入される。

bccsには、二重頂点連結成分を誘導する頂点集合の集合が代入される。

前提条件

  • arts, bccsは空にして渡すこと。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

int n;
bool g[N][N];

int bcc_decompose(int a, int depth, vector<int>& labels, vector<int>& comp,
                  vector<int>& arts, vector< vector<int> >& bccs) {
    int children = 0, upward = labels[a] = depth;
    REP(b, n) if (g[a][b]) {
        int u = labels[b];
        if (u < 0) {
            int k = comp.size();
            u = bcc_decompose(b, depth+1, labels, comp, arts, bccs);
            if (u >= depth) {
                comp.push_back(a);
                bccs.push_back(vector<int>(comp.begin()+k, comp.end()));
                comp.erase(comp.begin()+k, comp.end());
                if (depth > 0 || children > 0)
                    arts.push_back(a);
            }
        }
        upward <?= u;
        children++;
    }
    comp.push_back(a);
    if (depth == 0 && children == 0)
        bccs.push_back(comp);
    return upward;
}

void bcc_decompose(vector<int>& arts, vector< vector<int> >& bccs) {
    vector<int> comp, labels(n, -1);
    REP(r, n) if (labels[r] < 0) {
        comp.clear();
        bcc_decompose(r, 0, labels, comp, arts, bccs);
    }
}

検証済み

  • PKU2942 Knights of the Round Table

連結成分/強連結成分分解

Tarjanによる線形時間の強連結成分分解アルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

Graph g;
vector<bool> visited;

vector< vector<int> > scc() {
    int n = g.size();
    Graph r(n);   // make reversed graph
    REP(a, n) FOR(it, g[a])
        r[it->dest].push_back((Edge){it->dest, it->src});
    vector<int> order;
    visited.assign(n, false);
    REP(i, n) dfs(i, order);
    reverse(order.begin(), order.end());
    g.swap(r);
    vector< vector<int> > components;
    visited.assign(n, false);
    REP(i, n)
        if (!visited[order[i]])
            components.push_back(vector<int>()) ,
                dfs(order[i], components.back());
    g.swap(r);
    return components;
}
void dfs(int a, vector<int>& order) {
    if (visited[a])
        return;
    visited[a] = true;
    FOR(it, g[a]) dfs(it->dest, order);
    order.push_back(a);
}

未検証

メモ

  • 出力される強連結成分はDAGに縮約したときのトポロジカル順序に従っている。

連結成分/橋

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

int bridge(Graph& g, int here, int parent, vector<int>& labels,
           vector<int>& current, vector< vector<int> >& comps, int& id) {
    int& myid = labels[here];
    assert(myid < 0);
    myid = id++;

    int res = myid;
    Edges& v = g[here];
    REP(i, v.size()) {
        int there = v[i];
        if (there == parent)
            continue;
        int child;
        if (labels[there] < 0) {
            child = bridge(g, there, here, labels, critical, id);
            if (child > myid)
                critical[here][there] = critical[there][here] = true;
        }
        else {
            child = labels[there];
        }
        res <?= child;
    }

    return res;
}

未検証

連結成分/関節点

取り除いたときにグラフが非連結になるようなノード(関節点)を求めるアルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

typedef vector<int> Edges;
typedef vector<Edges> Graph;

Graph g;
vector<int> depths;
vector<int> counts;

void count() {
    int n = g.size();
    depths.assign(n, -1);
    counts.assign(n, 0);
    REP(a, n) {
        if (depths[a] >= 0)
            continue;
        dfs(a, 0);
        if (counts[a] > 0)
            counts[a]--;
    }
}
int dfs(int a, int cur) {
    depths[a] = cur;
    int upper = cur;
    FOR(it, g[a]) {
        int b = *it;
        int u = depths[b];
        if (u < 0) {
            u = dfs(b, cur+1);
            if (u >= cur)
                counts[a]++;
        }
        upper = min(upper, u);
    }
    return upper;
}

未検証

メモ

  • 戻り値は各ノードが取り除かれたときに新たにできる分裂の数に対応する。
  • TODO: 二重連結成分分解

巡回路/TSP

巡回セールスマン問題。NP困難なので指数時間のアルゴリズムしか知られていないが、DPを用いることによりN<=20程度であれば比較的短時間で求められる。

前提条件

  • 始点sが指定され、sから出発してsに戻ってくる巡回路を求める。
1
2
3
4
5
6
7
8
9
10
11
12

Weight adj[N][N];
Weight tsp[1<<N][N];
REP(i, 1<<n) REP(j, n)
    tsp[i][j] = WEIGHT_INFTY;
tsp[1<<s][s] = 0; // sから始まるpath
REP(i, 1<<n) REP(j, n) if (i&(1<<j))
    REP(k, n)
        tsp[i|(1<<k)][k] <?= tsp[i][j] + adj[j][k];
Weight* last = tsp[(1<<n)-1];
REP(j, n)
    last[j] += adj[j][s]; // sに戻る場合
return *min_element(last, last+n);

検証済み

  • UVA10496 Collecting Beepers
  • PKU2688 Cleaning Robot

メモ

  • O(n2 2n)

巡回路/有向オイラー路

グラフの一筆書き(オイラー路)を求めるアルゴリズム。

有向グラフが辺連結なとき、全ての点において出次数と入次数が等しい場合有向オイラー閉路が存在し、出次数1の点と入次数1の点が各1個でそれ以外の全ての点で出次数と入次数が等しい場合オイラー路が存在する。

前提条件

  • グラフは隣接リストで与えられること
  • 多重辺・セルフループを許す
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

void dfs(Graph& g, int a, vector<int>& route) {
  while(!g[a].empty()){
    int b = g[a].back();
    g[a].pop_back();
    dfs(g, b, route);
  }
  route.push_back(a);
}

vector<int> euler(Graph& g, int start) {
  int nNodes = g.size();
  int nEdges = 0;
  vector<int> deg(nNodes, 0);
  REP(i, nNodes){
    nEdges += g[i].size();
    REP(j, g[i].size())
      deg[g[i][j]]--; // 入次数
    deg[i] += g[i].size(); // 出次数
  }

  vector<int> route;
  int nonzero = nNodes - count(deg.begin(), deg.end(), 0);
  if(nonzero == 0 || (nonzero == 2 && deg[start] == 1)) {
    Graph g_(g);
    dfs(g_, start, route);
  }

  if((int)route.size() != nEdges + 1) // 非連結だった場合
    route.clear();
  reverse(route.begin(), route.end());
  return route;
}

検証済み

  • UVA10129 Play on Words

巡回路/無向オイラー路

グラフの一筆書き(オイラー路)を求めるアルゴリズム。

無向グラフが辺連結なとき、奇点が存在しない場合無向オイラー閉路が存在し、奇点が2個だけ存在する場合無向オイラー路が存在する。

前提条件

  • グラフgは隣接行列で与えられること
  • 多重辺・セルフループを許す
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

int g[N][N];
int adj[N][N];

void dfs(int a, vector<int>& route) {
  REP(b, N){
    if(adj[a][b]){
      adj[a][b]--;
      adj[b][a]--;
      dfs(b, route);
    }
  }
  route.push_back(a);
}

vector<int> euler(int start) {
  int odd = 0;
  int nEdges = 0;
  REP(i, N){
    int deg = accumulate(g[i], g[i] + N, 0);
    if(deg % 2 == 1)
      odd++;
    nEdges += deg;
  }
  nEdges /= 2;

  vector<int> route;
  int startdeg = accumulate(g[start], g[start] + N, 0);
  if(odd == 0 || (odd == 2 && startdeg % 2 == 1)){
    memcpy(adj, g, sizeof(g));
    dfs(start, route);
  }

  if((int)route.size() != nEdges + 1) // 非連結だった場合
    route.clear();
  reverse(route.begin(), route.end());
  return route;
}

検証済み

  • UVA10054 The Necklace

その他/二部グラフの辺彩色

二部グラフの辺に彩色をし、同一ノードに繋がっている辺が同じ色を持たないようにする。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

int n, m;                 // 左右のノード数
vector< vector<int> > g;  // グラフ (index -> [index])
vector< vector<int> > cl; // 色づけ (index -> color -> index)

BGEC(int n, int m) : n(n), m(m), g(n+m) {}
void add_edge(int a, int b) { // 左a、右bの間に辺を張る
    g[a].push_back(n+b);
    g[n+b].push_back(a);
}
int color() {
    int d = 0; // 彩色数
    REP(i, n+m)
        d = max<int>(d, g[i].size());
    cl.assign(n+m, vector<int>(d, -1));
    REP(a, n) REP(i, g[a].size()) {
        int b = g[a][i];
        int ca = min_element(ALLOF(cl[a])) - cl[a].begin();
        int cb = min_element(ALLOF(cl[b])) - cl[b].begin();
        if (ca != cb) {
            augment(b, ca, cb);
            cb = ca;
        }
        cl[a][ca] = b;
        cl[b][cb] = a;
    }
    return d;
}
void augment(int a, int c1, int c2) {
    int b = cl[a][c1];
    if (b >= 0) {
        augment(b, c2, c1);
        cl[b][c2] = a;
        cl[a][c1] = -1;
    }
    cl[a][c2] = b;
}

検証済み

  • UVA10615 Rooks

メモ

  • 最小彩色数はグラフの最大次数に等しい。

その他/範囲の更新

範囲の更新を行う。

前提条件

  • Intervalsには初期値として「負無限大=無色」を代入しておくこと。
1
2
3
4
5
6
7
8
9
10
11
12

typedef map<int, int> Intervals;

// paint [left, right) with color c
void paint(Intervals& v, int left, int right, int c) {
    assert(left < right);
    Intervals::iterator first = v.lower_bound(left), last = v.upper_bound(right);
    int p = (--last)->second;
    last++;
    v.erase(first, last);
    v[left] = c;
    v[right] = p;
}

検証済み

  • PKU2528 Mayor’s posters

メモ

  • Makegumiのライブラリによる
  • ならし時間で1オペレーションあたり平均O(log n)

数論

整数/エラトステネスの篩

1
2
3
4
5
6
7
8
9
10
11
12
13
14

#define MAX_PRIME 1000000
bool isprime[MAX_PRIME+1];
vector<int> primes;
void init_prime() {
    fill(isprime, isprime+MAX_PRIME+1, true);
    isprime[0] = isprime[1] = false;
    for(int i = 2; i <= MAX_PRIME; i++) {
        if (isprime[i]) {
            primes.push_back(i);
            for(int j = i*2; j <= MAX_PRIME; j += i)
                isprime[j] = false;
        }
    }
}

検証済み

  • PKU2689 Prime Distance

メモ

  • 二段目のループの初期条件は j = i*i で良いが、オーバーフローの危険がある。

整数/オイラーのφ関数

φ(n) = {自然数nと互いに素なn以下の数の個数}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

int phi(int n) {
    int res = n;
    for(int i = 2; i*i <= n; i++) {
        if (n % i == 0) {
            res -= res / i;
            while(n % i == 0)
                n /= i;
        }
    }
    if (n > 1) // n is prime
        res -= res / n;
    return res;
}

// phi(1)..phi(N) を求める
void phi_all(int N) {
    int a[N+1], b[N+1];
    REP(i, N+1) {
        a[i] = 1;
        b[i] = i;
    }

    REP(k, N+1) {
        if (b[k] < 2)
            continue;
        for(int n = k; n <= N; n += k) {
            for(int m = k-1; b[n]%k == 0; m = k) {
                a[n] *= m;
                b[n] /= k;
            }
        }
    }
    // a[n] == phi(n), b[n] == 1
}

未検証

メモ

  • 出典: echizenのWiki
  • 公式としては n = p1^r1 * p2^r2 * ... * pk^rk と素因数分解されるとき φ(n) = n * (1 - 1/p1) * (1 - 1/p2) * ... * (1 - 1/pk)

整数/中国剰余定理(1)

中国剰余定理とは、n個の整数 m[0], …, m[n-1] がどの2つの要素も互いに素ならば、与えられる r[0], …, r[n-1] に対して x == r[i] (mod m[i]) を満たすようなxが、Πm[i]を法として唯一つ存在する、というもの。

以下はそのようなxを 0 <= x <= Πm[i] の範囲で求めるアルゴリズム。但し実装時にはオーバーフローに注意すること。

前提条件

  • m[i]はどの2つの要素も互いに素であること
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

// BE CAREFUL OF OVERFLOW!!
integer chinese_remainder(const vector<int>& m, const vector<int>& r) {
  int n = m.size();
  integer prod = 1;
  REP(i, n)
    prod *= m[i];

  integer res = 0;
  REP(i, n){
    integer M = prod / m[i];
    integer a = divide(M, 1, m[i]);
    integer R = r[i] - r[i] / prod * prod;
    if(R < 0)
      R += prod;
    res = (res + M * a * r[i] % prod) % prod;
  }

  return res;
}

未検証

整数/中国剰余定理(2)

2つの制約「K=ax+b」「K=cy+d」に対して、それらを同時に満たすKは「K=ez+f」という一つの式で表せる(か、もしくはそのようなKは存在しない)。ただし、e=lcm(a,c)。aとcは互いに素でなくてもよい。

以下はそのようなe,fを求めるアルゴリズム。fの値の正規化方法は臨機応変に。オーバーフローに注意すること。

N個の制約について求めたい場合はN-1回呼んでください。もしくは、1x+0に対してN回fold。

前提条件

  • 互いに素である必要はない。
  • が、答えがない場合に注意
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

// ax+bのaとb
struct Constraint {
  int mult;
  int base;
};

// BE CAREFUL OF OVERFLOW!!
Constraint chinese_remainder(const Constraint& a, const Constraint& b) {
  Constraint r;
  int g = gcd(a.mult, b.mult);
  int d = b.base - a.base;

  // 解がない場合はこのassertで落ちるので、マズい場合は適宜修正すること
  assert(d % g == 0);
  d /= g;

  // このへんでオーバーフローに注意
  r.mult = a.mult / g * b.mult;
  r.base = a.mult * d * inv(a.mult/g, b.mult/g) + a.base;

  // ここから解の正規化
  // 以下のコードではK=ax+b=cy+d=ez+fに対して、
  // x>=0,y>=0の条件がついていた場合にz>=0を必要十分とするようなfを構成している
  // whileループが長引きそうなら割り算にしたほうがいいかもね
  r.base %= r.mult;
  while(r.base < a.base || r.base < b.base)
    r.base += r.mult;

  return r;
}

検証済み

  • PKU3545 K’ak’-u-pakal and Mayan Script

整数/拡張ユークリッド互除法

整数a, bに対して、ax + by = gcd(a, b) となる整数x, yを求める。

1
2
3
4
5
6
7
8
9
10

void xgcd(integer a, integer b, integer& x, integer& y) {
    if (b == 0) {
        x = 1;
        y = 0;
    }
    else {
        xgcd(b, a%b, y, x);
        y -= a/b*x;
    }
}

検証済み

  • UVA700 Date Bugs

整数/最大公約数・最小公倍数

1
2
3
4
5
6

integer gcd(integer a, integer b) {
    return (b == 0 ? a : gcd(b, a%b));
}
integer lcm(integer a, integer b) {
    return a/gcd(a, b)*b;
}

検証済み

  • UVA700 Date Bugs

整数/素数を法とした逆数

素数pを法とした剰余体上でのaの逆数を求める。

前提条件

  • pは素数
  • 0 < a < p
1
2
3

int inv(int a, int p) {
    return ( a == 1 ? 1 : (1 - p*inv(p%a, a)) / a + p );
}

メモ

  • GCDと同じ計算量
  • pが素数でない場合、aの逆元が存在すればそれを返す。aに0を指定したり、逆元が存在しなければdivision by zeroで落ちる。

整数/線形ディオファントス方程式

an = b mod m となる、非負でかつ最小のnを求める。 もしくは、mを法とした剰余環上でb/aを計算する。 解がない場合はno_solutionをthrowする。 a, b, mは非負である必要がある。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

struct no_solution {};
integer divide(integer a, integer b, integer m) {
    integer g = gcd(a, m);
    if (b%g != 0)
        throw no_solution();
    integer x, y;
    xgcd(a, m, x, y);
    assert(a*x+m*y == gcd(a,m));
    integer n = x*b/g;
    integer dn = m/g;
    n -= n/dn*dn;
    if (n < 0)
        n += dn;
    return n;
}

検証済み

  • UVA700 Date Bugs

行列演算/べき乗

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

matrix_t pow(matrix_t& e, int n, int m) {

    matrix_t r = new_matrix(n);

    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            r[i][j] = ((m&1) == 0 ? (i == j ? 1 : 0) : e[i][j]);

    if (m >= 2) {
        matrix_t u = pow(e, n, m/2);
        matrix_t uu = multiply(u, u, n);
        matrix_t z = multiply(r, uu, n);
        delete_matrix(u, n);
        delete_matrix(uu, n);
        delete_matrix(r, n);
        r = z;
    }

    return r;
}

未検証

行列演算/基本演算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

typedef double* vector_t;
typedef vector_t* matrix_t;
matrix_t new_matrix(int n) {
    matrix_t x = new vector_t[n];
    for(int i = 0; i < n; i++)
        x[i] = new double[n];
    return x;
}
matrix_t dup_matrix(matrix_t x_, int n) {
    matrix_t x = new_matrix(n);
    for(int i = 0; i < n; i++)
        copy(x_[i], x_[i]+n, x[i]);
    return x;
}
void delete_matrix(matrix_t x, int n) {
    for(int i = 0; i < n; i++)
        delete[] x[i];
    delete[] x;
}
matrix_t multiply(matrix_t a, matrix_t b, int n) {
    matrix_t r = new_matrix(n);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            r[i][j] = 0;
            for(int k = 0; k < n; k++)
                r[i][j] += a[i][k] * b[k][j];
        }
    }
    return r;
}

未検証

線形方程式/LU分解

aを破壊しLU形式に変換する。pは行交換の情報を保持する。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

bool lu_decompose(matrix_t& a, int* p, int n) {

    for(int i = 0; i < n; i++)
        p[i] = i;

    for(int k = 0; k < n; k++) {

        int pivot = k;
        for(int i = k+1; i < n; i++)
            if (abs(a[i][k]) > abs(a[pivot][k]))
                pivot = i;

        swap(a[k], a[pivot]);
        swap(p[k], p[pivot]);

        if (abs(a[k][k]) < EP)
            return false;

        for(int i = k+1; i < n; i++) {
            double m = (a[i][k] /= a[k][k]);
            for(int j = k+1; j < n; j++)
                a[i][j] -= a[k][j] * m;
        }

    }

    return true;
}
void lu_solve(matrix_t& a, int* p, vector_t& b, vector_t& x, int n) {

    for(int i = 0; i < n; i++)
        x[i] = b[p[i]];

    for(int k = 0; k < n; k++)
        for(int i = 0; i < k; i++)
            x[k] -= a[k][i] * x[i];

    for(int k = n-1; k >= 0; k--) {
        for(int i = k+1; i < n; i++)
            x[k] -= a[k][i] * x[i];
        x[k] /= a[k][k];
    }

}

検証済み

  • PKU2118 Firepersons
  • ZJU2525 Restore the Polygon

線形方程式/ガウスの消去法

連立方程式 ax=b を解く。

a, bは破壊され、答えがbに代入される。

前提条件

  • invert, moduloを別途定義すること。
  • A, bはmoduloによる正規化が既に行われているものとする。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

void gauss(matrix_t& A, vector_t& b, int n, int m) {
    int pi = 0, pj = 0;
    while(pi < n && pj < m) {
        for(int i = pi+1; i < n; i++) {
            if (abs(A[i][pj]) > abs(A[pi][pj])) {
                swap(A[i], A[pi]);
                swap(b[i], b[pi]);
            }
        }
        if (abs(A[pi][pj]) > 0) {
            int d = invert(A[pi][pj]);
            REP(j, m)
                A[pi][j] = modulo(A[pi][j] * d);
            b[pi] = modulo(b[pi] * d);
            for(int i = pi+1; i < n; i++) {
                int k = A[i][pj];
                REP(j, m)
                    A[i][j] = modulo(A[i][j] - k * A[pi][j]);
                b[i] = modulo(b[i] - k * b[pi]);
            }
            pi++;
        }
        pj++;
    }
    for(int i = pi; i < n; i++)
        if (abs(b[i]) > 0)
            throw Inconsistent();
    if (pi < m || pj < m)
        throw Ambiguous();
    for(int j = m-1; j >= 0; j--)
        REP(i, j)
            b[i] = modulo(b[i] - b[j] * A[i][j]);
}

検証済み

  • PKU2497 Widget Factory

メモ

  • O(n3).

線形方程式/単体法

二段階単体法により線形計画問題を解く。

ここで言う線形計画問題とは、Ax=bのもとでc.xを最小化する問題。

前提条件

  • 解が存在しない/最適値が発散する場合は vector_t() を返す。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

const double EPS = 1.0e-10;
const double INF = numeric_limits<double>::infinity();

typedef vector<double> vector_t;
typedef vector<vector_t> matrix_t;

vector_t simplex(matrix_t A, vector_t b, vector_t c) {

    const int n = c.size(), m = b.size();

    // modify b to non-negative
    REP(i, m) if (b[i] < 0) {
        REP(j, n)
            A[i][j] *= -1;
        b[i] *= -1;
    }

    // list of base/independent variable ids
    vector<int> bx(m), nx(n);
    REP(i, m)
        bx[i] = n+i;
    REP(i, n)
        nx[i] = i;

    // extend A, b
    A.resize(m+2);
    REP(i, m+2)
        A[i].resize(n+m, 0);
    REP(i, m)
        A[i][n+i] = 1;
    REP(i, m) REP(j, n)
        A[m][j] += A[i][j];
    b.push_back(accumulate(ALLOF(b), (double)0.0));
    REP(j, n)
        A[m+1][j] = -c[j];
    REP(i, m)
        A[m+1][n+i] = -INF;
    b.push_back(0);

    // main optimization
    REP(phase, 2) {

        for(;;) {

            // select an independent variable
            int ni = -1;
            REP(i, n)
                if (A[m][nx[i]] > EPS && (ni < 0 || nx[i] < nx[ni]))
                    ni = i;
            if (ni < 0)
                break;

            int nv = nx[ni];

            // select a base variable
            vector_t bound(m);
            REP(i, m)
                bound[i] = (A[i][nv] < EPS ? INF : b[i] / A[i][nv]);
            if (!(*min_element(ALLOF(bound)) < INF))
                return vector_t(); // -infinity

            int bi = 0;
            REP(i, m)
                if (bound[i] < bound[bi]-EPS || (bound[i] < bound[bi]+EPS && bx[i] < bx[bi]))
                    bi = i;

            // pivot
            double pd = A[bi][nv];
            REP(j, n+m)
                A[bi][j] /= pd;
            b[bi] /= pd;

            REP(i, m+2) if (i != bi) {
                double pn = A[i][nv];
                REP(j, n+m)
                    A[i][j] -= A[bi][j] * pn;
                b[i] -= b[bi] * pn;
            }

            swap(nx[ni], bx[bi]);
        }

        if (phase == 0 && abs(b[m]) > EPS)
            return vector_t(); // no solution

        A[m].swap(A[m+1]);
        swap(b[m], b[m+1]);
    }

    vector_t x(n+m, 0);
    REP(i, m)
        x[bx[i]] = b[i];
    x.resize(n);

    return x;
}

検証済み

  • UVA10498 Happiness
  • NETLIB Test Problems for LPの一部

メモ

  • 単体法は最悪の場合で指数時間がかかるので、あくまで最終手段。これを使う前に、ほかの方法が使えないか十分考えること。
  • 特に、2変数の不等式制約は二次元の凸多角形クリッピングとして捉えることができる。

文字列

検索/Aho-Corasick

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

struct AC {
    AC* fail;
    AC* next[4];
    vector<int> accepts;

    AC() : fail(0) {
        memset(next, 0, sizeof(next));
    }
    void insert(const int* p, int id) {
        if (*p < 0)
            accepts.push_back(id);
        else
            (next[*p] ?: (next[*p] = new AC()))->insert(p+1, id);
    }
    void make() {
        this->fail = 0; // this is root
        queue<AC*> q;
        q.push(this);
        while(!q.empty()) {
            AC* cur = q.front();
            q.pop();
            REP(k, 4) {
                AC*& next = cur->next[k];
                if (!next) {
                    next = (cur->fail ? cur->fail->next[k] : this);
                }
                else {
                    AC* fail = cur->fail;
                    while(fail && !fail->next[k])
                        fail = fail->fail;
                    fail = (fail ? fail->next[k] : this);
                    next->fail = fail;
                    next->accepts.insert(next->accepts.end(), ALLOF(fail->accepts));
                    q.push(next);
                }
            }
        }
    }
};

検証済み

  • PKU2778 DNA Sequence

検索/KMP

Knuth-Morris-Prattによる線形時間の文字列検索アルゴリズム。

検索するパターンについてあらかじめ線形時間でテーブルを計算しておくことで、テキストを線形時間で検索することができる。

ナイーブに文字列検索を行う場合、テキストとパターンの照合中に文字の不一致が見つかった時はテキストのポインタをひとつすすめ、パターンのポインタを先頭に戻すことになる。一方KMP法では、テキストのポインタを動かさずに、パターンのポインタを少し戻すことによってさらなるマッチングを試みる。ポインタを戻す距離は、パターンのポインタ以前とテキストがマッチしているような配置の中で動く距離が一番短いものを選ぶ。パターンの先頭文字とも一致しないような場合にのみテキストのポインタを進める。

KMP法の前処理では以上の処理に使う「文字が一致しなかったときにパターンのポインタをどこに戻すか」というテーブル(パターン中の不一致位置 -> パターンのポインタを戻す位置)を作成する。このテーブルは、実は「パターン中の先頭n文字を考えたときに、そのprefixとsuffixが一致するような部分長さ L (0 <= L < n) のうち最長のもの」のテーブルと同値になる。(L < n であることに注意。L = nのような自明なprefix/suffixは考えない。)

このテーブルTの計算方法を考えると、以下のようなDPが成立することがわかる:

1
2
3
4
5

T[p] = max (q + 1)
 s.t. s[p-1] = s[q]
      q = {q0, T[q0], T[T[q0]], ...}
      q0 = T[p-1]
T[0] = -1

ただし、T[0]は計算上の便宜的なもので、検索時には参照されない。

1
2
3
4
5
6
7
8
9
10
11

vector<int> kmp(string s) {
    int n = s.size();
    vector<int> table(n+1);
    int q = table[0] = -1;
    for(int p = 1; p <= n; p++) {
        while(q >= 0 && s[q] != s[p-1])
            q = table[q];
        table[p] = ++q;
    }
    return table;
}

検証済み

  • PKU1961 Period
  • PKU2406 Power Strings
  • PKU2752 Seek the Name, Seek the Fame

メモ

  • O(n + m)

Suffix Array/Karkkainen Sanders

  1. Kärkkäinen と P. Sanders による線形時間のSuffix Array構築アルゴリズム。

前提条件

  • テキストサイズ >= 2.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

// Derived from:
// J. Karkkainen and P. Sanders: Simple Linear Work Suffix Array Construction.
// In 30th International Colloquium on Automata, Languages and Programming,
// number 2719 in LNCS, pp. 943--955, Springer, 2003.

void radix_sort(int* src, int* dest, int* rank, int n, int m) {
    int* cnt = new int[m+1];
    memset(cnt, 0, sizeof(int)*(m+1));
    for(int i = 0; i < n; i++)
        cnt[rank[src[i]]]++;
    for(int i = 1; i <= m; i++)
        cnt[i] += cnt[i-1];
    for(int i = n-1; i >= 0; i--)
        dest[--cnt[rank[src[i]]]] = src[i];
}

bool triple_less(int a0, int a1, int a2, int b0, int b1, int b2) {
    if (a0 != b0)
        return (a0 < b0);
    if (a1 != b1)
        return (a1 < b1);
    return (a2 < b2);
}

void suffix_array(int* text, int* sa, int n, int m) {
    assert(n >= 2);

    int n0 = (n+2)/3, n1 = (n+1)/3, n2 = n/3, n02 = n0 + n2;

    int* rank12 = new int[n02+3]; rank12[n02+0] = rank12[n02+1] = rank12[n02+2] = 0;
    int* sa12 = new int[n02+3]; sa12[n02+0] = sa12[n02+1] = sa12[n02+2] = 0;
    int* tmp0 = new int[n0];
    int* sa0 = new int[n0];

    for(int i = 0; i < n02; i++)
        rank12[i] = i*3/2+1;
    radix_sort(rank12, sa12, text+2, n02, m);
    radix_sort(sa12, rank12, text+1, n02, m);
    radix_sort(rank12, sa12, text+0, n02, m);

    int name = 0;
    for(int i = 0; i < n02; i++) {
        if (name == 0 ||
            text[sa12[i]+0] != text[sa12[i-1]+0] ||
            text[sa12[i]+1] != text[sa12[i-1]+1] ||
            text[sa12[i]+2] != text[sa12[i-1]+2])
            name++;
        rank12[sa12[i]/3 + (sa12[i]%3 == 1 ? 0 : n0)] = name;
    }

    if (name < n02) {
        suffix_array(rank12, sa12, n02, name);
        for(int i = 0; i < n02; i++)
            rank12[sa12[i]] = i+1;
    }
    else {
        for(int i = 0; i < n02; i++)
            sa12[rank12[i]-1] = i;
    }

    for(int i = 0, j = 0; i < n02; i++)
        if (sa12[i] < n0)
            tmp0[j++] = 3*sa12[i];
    radix_sort(tmp0, sa0, text, n0, m);

    for(int i = 0, j = n0-n1, k = 0; k < n; ) {
        int a = sa0[i];
        int b = (sa12[j] < n0 ? sa12[j]*3+1 : (sa12[j]-n0)*3+2);
        if (b%3 == 1 ? triple_less(text[a], rank12[a/3], 0,
                                   text[b], rank12[b/3+n0], 0) :
                       triple_less(text[a], text[a+1], rank12[a/3+n0],
                                   text[b], text[b+1], rank12[b/3+1])) {
            sa[k++] = a;
            if (++i == n0)
                while(k < n)
                    sa[k++] = (sa12[j] < n0 ? sa12[j++]*3+1 : (sa12[j++]-n0)*3+2);
        }
        else {
            sa[k++] = b;
            if (++j == n02)
                while(k < n)
                    sa[k++] = sa0[i++];
        }
    }

    delete[] rank12; delete[] sa12; delete[] tmp0; delete[] sa0;
}

検証済み PKU3581 Sequence

メモ

Suffix Array/Larsson Sadakane

Larsson と Sadakane によるSuffix Array構築アルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

int* suffix_array(unsigned char* str) {
    int n = strlen((char*)str);
    int *v, *u, *g, *b, r[256];
    v = new int[n+1]; u = new int[n+1]; g = new int[n+1]; b = new int[n+1];
    for(int i = 0; i <= n; i++) {
        v[i] = i;
        g[i] = str[i];
    }
    for(int h = 1; ; h *= 2) {
        int m = *max_element(g, g+n+1);
        for(int k = h; k >= 0; k -= h) {
            for(int ord = 0; m >> ord; ord += 8) {
                memset(r, 0, sizeof(r));
                for(int i = 0; i <= n; i++)
                    r[(g[min(v[i]+k, n)] >> ord) & 0xff]++;
                for(int i = 1; i < 256; i++)
                    r[i] += r[i-1];
                for(int i = n; i >= 0; i--)
                    u[ --r[(g[min(v[i]+k, n)] >> ord) & 0xff] ] = v[i];
                swap(u, v);
            }
        }
        b[0] = 0;
        for(int i = 1; i <= n; i++)
            b[i] = b[i-1] + ((g[v[i-1]] != g[v[i]] ? g[v[i-1]] < g[v[i]] : g[v[i-1]+h] < g[v[i]+h]) ? 1 : 0);
        if (b[n] == n)
            break;
        for(int i = 0; i <= n; i++)
            g[v[i]] = b[i];
    }
    delete[] g; delete[] b; delete[] u;
    return v;
}

未検証

メモ

  • たぶん O(n log n).

Suffix Array/高さ配列

Kasaiらによる、構築されたSuffix Arrayについて、順位が隣り合った二つのsuffixのlongest common prefixを求めるアルゴリズム。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

int* height_array(char* text, int* sa, int n) {
    int* height = new int[n];
    int* rank = new int[n];
    REP(i, n)
        rank[sa[i]] = i;
    int h = 0;
    height[0] = 0;
    REP(i, n) {
        if (rank[i] > 0) {
            int k = sa[rank[i]-1];
            while(text[i+h] == text[k+h])
                h++;
            height[rank[i]] = h;
            if (h > 0)
                h--;
        }
    }
    delete[] rank;
    return height;
}

未検証

メモ

その他

平衡木/Range Tree

区間に対する加算、区間中の最小値のクエリ、各要素の値のクエリを対数程度の時間で行う木。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

struct range_tree {
    int m;
    vector< pair<int, int> > tree;

    range_tree(int n) {
        m = n;
        while(m&(m-1))
            m += m&-m;
        tree.assign(2*m, make_pair(0, 0));
        for(int k = n; k < m; k++)
            tree[k] = make_pair(INF, INF);
        for(int k = m; k < 2*m-1; k++)
            tree[k] = make_pair(0,
                                min(tree[((k&~m)<<1)^0].second,
                                    tree[((k&~m)<<1)^1].second));
    }
    void range_add(int a, int b, int d) {
        while(a < b) {
            int r = 1, k = a;
            while((a & r) == 0 && a + (r<<1) <= b) {
                r <<= 1;
                k = (k >> 1) | m;
            }
            tree[k].first += d;
            do {
                tree[k].second = tree[k].first +
                    (k < m ? zero :
                     min(tree[((k&~m)<<1)^0].second,
                         tree[((k&~m)<<1)^1].second));
                k = (k >> 1) | m;
            } while(k < 2*m-1);
            a += r;
        }
    }
    int range_min(int a, int b) {
        int res = INF;
        while(a < b) {
            int r = 1, k = a;
            while((a & r) == 0 && a + (r<<1) <= b) {
                r <<= 1;
                k = (k >> 1) | m;
            }
            res = min(res, tree[k].second);
            a += r;
        }
        return res;
    }
    int query(int k) {
        int res = 0;
        while(k < 2*m-1) {
            res += tree[k].first;
            k = (k >> 1) | m;
        }
        return res;
    }
};

未検証

メモ

  • range_add に O((log n)2).
  • range_min に O(log n).
  • query に O(log n).

平衡木/Treap

ノードにkeyと一緒にpriorityなる値を割り当てる。Treapはkeyに関して一般のBSTと同様の順序関係を維持するが、それに加えてpriorityに関してheap条件を満たすようにする。任意のノード集合(key/priorityペアの集合)に対して、二つの条件を同時に満たすような木は必ず存在する。また、keyおよびpriorityがそれぞれdistinctである場合は、木は一意に定まる。

その性質から、priorityをランダムに割り当てることによって木の形をランダムにし、確率的に平衡させることができる。

Treapのふたつの条件を維持するような回転操作は1段の右回転・左回転操作のみで行うことができる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

template<class Key, class Value, class Cache>
struct Treap {
    Key key;
    Value value;
    int prio;
    Treap* ch[2];
    bool cached;
    Cache cache;
    Treap(const Key& key, const Value& value) : key(key), value(value), prio(rand()), cached(false) {
        ch[0] = ch[1] = 0; // is this necessary?
    }
    Treap* insert(const Key& newkey, const Value& newvalue) {
        if (!this)
            return new Treap(newkey, newvalue);
        if (newkey == key)
            return this;
        int side = (newkey < key ? 0 : 1);
        ch[side] = ch[side]->insert(newkey, newvalue);
        cached = false;
        return rotate(side);
    }
    Treap* rotate(int side) {
        if (ch[side] && ch[side]->prio > prio) {
            Treap* rot = ch[side];
            this->ch[side] = rot->ch[side^1];
            rot->ch[side^1] = this;
            this->cached = rot->cached = false;
            return rot;
        }
        return this;
    }
    void clear() {
        if (this) {
            ch[0]->clear();
            ch[1]->clear();
            delete this;
        }
    }
    Treap* remove(const Key& oldkey) {
        if (!this)
            return this;
        if (key == oldkey) {
            if (!ch[1]) {
                Treap* res = ch[0];
                delete this;
                return res;
            }
            pair<Treap*, Treap*> res = ch[1]->delete_min();
            res.first->ch[0] = this->ch[0];
            res.first->ch[1] = res.second;
            res.first->cached = false;
            delete this;
            return res.first->balance();
        }
        int side = (oldkey < key ? 0 : 1);
        ch[side] = ch[side]->remove(oldkey);
        cached = false;
        return this;
    }
    pair<Treap*, Treap*> delete_min() {
        cached = false;
        if (!ch[0])
            return make_pair(this, ch[1]);
        pair<Treap*, Treap*> res = ch[0]->delete_min();
        ch[0] = res.second;
        return make_pair(res.first, this);
    }
    inline int priority() {
        return (this ? prio : -1);
    }
    Treap* balance() {
        int side = (ch[0]->priority() > ch[1]->priority() ? 0 : 1);
        Treap* rot = rotate(side);
        if (rot == this)
            return this;
        return rot->balance();
    }
    Cache eval() {
        if (!this)
            return Cache();
        if (!cached)
            cache = Cache(key, value, ch[0]->eval(), ch[1]->eval());
        cached = true;
        return cache;
    }
};

// 小さいほうからn番目の要素をクエリできるようにする場合の例
struct SizeCache {
    int size;
    SizeCache() : size(0) {}
    template<class Key, class Value>
    SizeCache(const Key& key, const Value& value, const SizeCache& left, const SizeCache& right)
     : size(left.size+right.size+1) {}
};

template<class Key, class Value>
pair<Key, Value> nth(Treap<Key, Value, SizeCache>* root, int k) {
    int l = root->ch[0]->eval().size;
    if (k < l)
        return nth<Key, Value>(root->ch[0], k);
    if (k == l)
        return make_pair(root->key, root->value);
    return nth<Key, Value>(root->ch[1], k-(l+1));
}

検証済み

  • PKU2182 Lost Cows

メモ

  • サブツリーに対する操作を行ってはいけない。

列/Fenwick Tree

配列のpartial sumの取得と要素の書き換えをそれぞれ対数時間で行うデータ構造。別名Binary Indexed Tree。

n次元にも容易に拡張できる。

1
2
3
4
5
6
7
8
9
10
11
12
13

T bitquery(const vector<T>& bit, int from, int to) { // [from, to)
    if (from > 0)
        return bitquery(bit, 0, to) - bitquery(bit, 0, from);
    T res = T();
    for(int k = to-1; k >= 0; k = (k & (k+1)) - 1)
        res += bit[k];
    return res;
}

void bitupdate(vector<T>& bit, int pos, const T& delta) {
    for(const int n = bit.size(); pos < n; pos |= pos+1)
        bit[pos] += delta;
}

検証済み

  • PKU1195 Mobile Phone

メモ

  • 必要メモリはO(n). 普通に配列を持つのと同じサイズ。
  • 和の取得・要素の更新ともにO(log n).

列/Range Minimum Query

列のある区間の最小値を対数時間で求めるデータ構造。

前提条件

  • 配列サイズは２のべき乗であること。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

// 配列を拡張してRMQに対応させる
void rmq_ext(vector<T>& v) {
    int n = v.size();
    assert(__builtin_popcount(n) == 1); // 長さは2のべき乗でなければならない
    v.resize(n*2);
    for(int i = n; i < 2*n; i++)
        v[i] = min(v[(i-n)*2+0], v[(i-n)*2+1]);
}
// 列の要素を書き換える
void rmq_update(vector<T>& rmq, int pos, const T& value) {
    int n = rmq.size() / 2;
    rmq[pos] = value;
    while(pos < 2*n-1) {
        rmq[pos/2+n] = min(rmq[pos], rmq[pos^1]);
        pos = pos/2+n;
    }
}
// [from, to)の最小値を取り出す
T rmq_query(const vector<T>& rmq, int from, int to) {
    int n = rmq.size() / 2;
    int p = min((from == 0 ? 32 : __builtin_ctz(from)), 31-__builtin_clz(to-from));
    T x = rmq[(from>>p)|((n*2*((1<<p)-1))>>p)];
    from += 1<<p;
    if (from < to)
        x = min(x, rmq_query(rmq, from, to));
    return x;
}

未検証

メモ

列/最長共通部分列

Longest Common Sequenceは、二つの列に共通な部分列のうち最長のもののこと。

1
2
3
4
5
6

int lcs[N+1][M+1];
REP(i, n+1) REP(j, m+1) lcs[i][j] = 0;
REP(i, n) REP(j, m)
    lcs[i+1][j+1] = (str1[i] == str2[j] ? lcs[i][j]+1 : 0)
                    >? lcs[i+1][j] >? lcs[i][j+1];
return lcs[n][m];

検証済み

  • UVA10192 Vacation

メモ

  • 二つの列の長さをN, Mとして、O(NM)。

列/最長増加部分列

最長増加部分列(Longest Increasing Sequence)とは、数列の部分列のうち単調増加であるような最長のもの。

1
2
3
4

int lis[N];
REP(i, n) lis[i] = INF;
REP(i, n) *upper_bound(lis, lis+n, v[i]) = v[i];
return find(lis, lis+n, INF) - lis;

メモ

  • O(n log n)。
  • 減少部分列を求めたい場合はあらかじめ数列をnegateするか、二分探索のcomparatorにgreater<T>()を渡せば良い。
  • コードは厳密には最長非減少部分列の場合。strictな増加列を求める場合はlower_boundを使えば良い。

日付/Zellerの公式

グレゴリオ暦における指定された日付の曜日を求める公式。

1
2
3
4
5
6
7

int zeller(int y, int m, int d) {
    if (m <= 2) {
        y--;
        m += 12;
    }
    return (y + y/4 - y/100 + y/400 + (13 * m + 8) / 5 + d) % 7;
}

未検証

ゲーム/Grundy数

盤面ひとつひとつに数字を割り当てる。後手必勝が0、先手必勝が1以上。

選択はmex。平行していくつかのゲームを進める場合はxor。

ゲーム/αβ枝狩り

なんだかんだで毎回書くのに苦労するあれ。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

// [alpha..beta]の値以外を返しても意味を持たない、という意味
int search(int alpha = -INF, int beta = INF) {
    if (finished())
        return 0;
    if (enemy_turn())
        return -search(-beta, -alpha);
    for(move a : possible moves) {
        int p = gain(a);
        alpha >?= p + search(alpha-p, beta-p);
        if (alpha >= beta)
            break;
    }
    return alpha;
}

組み合わせ最適化/マトロイド交差

マトロイド交差の最大独立集合を求める。

前提条件

  • M1, M2について独立性判定ができること(M::appendable)。
  • M2について、独立な集合に要素を加えて非独立になったとき、存在する一意の回路を計算できること(M2::circuit)。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

template<class E, class M1, class M2>
set<E> augment(M1 m1, M2 m2, set<E> g, set<E> s) {
    map<E, E> trace;

    vector<E> rights;
    FOR(i, g) {
        if (m1.appendable(s, *i)) {
            rights.push_back(*i);
            trace[*i] = *i;
        }
    }
    while(!rights.empty()) {
        vector<E> lefts;
        REP(i, rights.size()) {
            E e = rights[i];
            if (m2.appendable(s, e)) {
                while(trace[e] != e) {
                    s.insert(e); e = trace[e];
                    s.erase(e); e = trace[e];
                }
                s.insert(e);
                return s;
            }
            vector<E> c = m2.circuit(s, e);
            REP(j, c.size()) {
                E f = c[j];
                if (trace.count(f) == 0) {
                    trace[f] = e;
                    lefts.push_back(f);
                }
            }
        }
        rights.clear();
        REP(i, lefts.size()) {
            E f = lefts[i];
            s.erase(f);
            FOR(i, g) {
                E e = *i;
                if (m1.appendable(s, e) && trace.count(e) == 0) {
                    trace[e] = f;
                    rights.push_back(e);
                }
            }
            s.insert(f);
        }
    }
    return s;
}

未検証

  • ACM-ICPC OB/OGの会 冬合宿2007 Day3 E Networking Company

メモ

  • やってみたかっただけ。コードは役に立たない。

NP完全問題/ナップサック問題

各アイテムの価値と大きさと容器のキャパシティが与えられ、価値の合計を最大にするようなアイテムの選び方を探す問題をナップサック問題と呼ぶ。

NP完全問題であるが、アイテムの大きさが整数である場合に擬多項式時間アルゴリズムが存在する。

前提条件

  • アイテムは各1個とする。
1
2
3
4
5

Weight v[C+1];
REP(i, c+1) v[i] = 0;
REP(k, n) for(int i = c; i >= costs[k]; i--)
    v[i] >?= v[i - costs[k]] + values[k];
return *max_element(v, v+c+1);

検証済み

  • UVA10130 Supersale

メモ

  • アイテム数をn, キャパシティをCとして O(nC)。
  • アイテムの数が無制限な場合はテーブルの走査を順方向に行えば良い。

チェックリスト

  • とりあえず深呼吸。
  • 思いつかないとき
    • これフローじゃね？
    • これ二分探索じゃね？
    • これ凸じゃね？
  • 最小/最大/特異ケースのテスト
  • freopen!!!
  • 問題の条件をよみなおせ！
    • 入力の形式はあっているか？
    • 入力ファイルを修正してそのままにしていないか
    • スペースを含む文字列
    • 0-origin, 1-origin
    • (i,j) <-> (x,y) 変換をはじめとする座標変換
    • 改行コードwww
  • 出力の形式はあっているか？
    • スペルミス
    • -0.000
    • 空行のはさみ方
    • デバッグ出力の削りのこし
    • 解が無い時
  • オーバーフロー
    • 答えがintに収まらない
    • numeric_limits<int>::max() * 2
    • bool型の変数にint/doubleを代入、関数の返り値型がbool (特に、DP/メモの型を後で変更した場合)
  • ポインタの進め忘れ
  • 参照/イテレータの無効化
  • 実数
    • EPSを使わずに比較していないか
    • 問題に合ったEPSの設定
    • NaN (/, sqrt, asin, acos, etc…)
  • switch
    • breakし忘れていないか
    • defaultにassert(false)
  • ライブラリ関数
    • 全順序関係でないようなoperator<を使ったソート
    • accumulate/inner_productの初期値は必ず型をキャストで明示すること
  • 配列
    • サイズは十分か？
    • 初期化を忘れていないか？
    • 番兵のせいでサイズ・インデックス・初期化が間違っていないか？
    • 添字を間違っていないか？(jと書くべきところにi、など)
  • コピペ
    • コピペ修正ミスは無いか？(y += dx みたいな)
  • 考え方
    • 問題を分解したときに独立にとけないこともあるので柔軟に考えるといいよ

謝辞

このライブラリでは数多くのサイト/ライブラリを参考にしています。各著者の方々に感謝します。(以下敬称略)

参考にしているサイト

参考にしているライブラリ