Algorithm library by Tomoki Imai documentation¶

テンプレート¶

C++¶

// compile in C++11. use -std=c++11.
#include <iostream>
#include <iomanip>
#include <vector>
#include <valarray>
#include <map>
#include <set>
#include <list>
#include <queue>
#include <stack>
#include <bitset>
#include <utility>
#include <numeric>
#include <algorithm>
#include <functional>
#include <complex>
#include <string>
#include <sstream>

#include <cstdio>
#include <cstdlib>
#include <cctype>
#include <cstring>

// this require C++11
#include <unordered_set>
#include <unordered_map>
#include <random>

using namespace std;

#define all(c) c.begin(),c.end()
#define repeat(i,n) for(int i=0;i<static_cast<int>(n);i++)
#define debug(x) #x << "=" << (x)
#define dump(x) cerr << debug(x) << " (L:" << __LINE__ << ")"<< endl

typedef long long ll;
typedef vector<int> vi;

template<typename T>
ostream& operator<<(ostream& os,const vector<T>& vec){
    os << "[";
    for(const auto& v : vec){
        os << v << ",";
    }
    os << "]";
    return os;
}

template<typename T>
T input(){
    T t;cin >> t;
    return t;
}

template<typename T>
vector<T> input(const int N){
    vector<T> v(N);
    repeat(i,N) cin >> v[i];
    return v;
}



int main(){
    return 0;
}

Python¶

#!/usr/bin/env python
#coding:utf-8

from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

from math import *
from string import *
from fractions import *
from itertools import *

def main():
    pass

if __name__ == "__main__":
    main()

算術型¶

int¶

基本中の基本。$10^{9}$くらい。こわいときにはlong longを使うことを推奨。

long long¶

大きい整数。$10^{18}$?くらい。

ビット演算¶

ビットDPとかに使う。

using ll = long long;
// 下から2番目のビットを立てる。1llと書くことに注意
x |= (1ll << 1);

// 集合演算
// 和
ll z = x | y;
// 積
ll z = x & y;
// xor
ll z = x ^ y;
// ループ
for(int i=0;i<N;i++){
    // i番目が立っていれば
    if(x & (1ll << i)){
        //something.
    }
}

ll dp[1 << N][N];
dp[0][0] = 1;
// Bitのほうが先!
for(int i=0;i<(1<<N);i++){
    for(int j=0;j<N;j++){
        //something.
    }
}

double¶

floatは使ってはだめ。

char¶

-128 ~ 127くらい。ちいさい。基本的には文字を入れるのに使う。vector<char>をvector<bool>の代わりに使ってもいい。

char c = 'a';
//cctypeが必要。大文字に変換する。すでに大文字のときは何も起こらない。
c = toupper(c);
//小文字に変換する。
c = tolower(c);
//vector<bool>の代わり。
vector<char> used(10,false);

bool¶

true(==1)とかfalse(==0)を入れるためだけに使う。ただし vector<bool> は特殊化されているため、注意が必要。 https://learn.microsoft.com/ja-jp/cpp/standard-library/vector-bool-class?view=msvc-170

補助関数¶

上記の型に関する便利な関数。

床関数、天井関数、および四捨五入。返り値はdouble。cmathをincludeする。

double x = 0.3;
int f = floor(x); // -> 0
int c = ceil(x); // -> 1
int r = round(x) // -> 0

任意の場所で四捨五入したいときには、$10^{n}$をかけて、roundした後に、 $10^{n}$で割ればいい。

double x = 0.123456789;
//0.123
double r = round(x*1000) / 1000.0;

入出力¶

基本はcin,coutを使おう。

cin,cout¶

iostream,iomanipをincludeしておくこと。

cin¶

基本的な使い方について。

int n;
cin >> n;
vector<int> V(n);
rep(i,n) cin >> V[i];

こうすると、短く書ける。

入力の最後まで読む。

int n;
while(cin >> n){
    //処理
}

n=0のとき終わりとかの場合は、条件に&&n!=0とかをつける。数値をカンマ区切りで読み込む。

int x,y;char c;
//cにカンマが入る
cin >> x >> c >> y;

冗長かもだけど、一番楽。空白とか含めて一行読む。

string s;
getline(cin,s);

改行文字は、sに入らず、かつ読み捨てされる。 cinでは、改行文字は読み捨てないことに注意しよう。つまり、数値<改行>文字列<改行>を読みたいときには、

int n;string s;
// 数値
cin >> n;
// 改行よみとばし
cin.ignore();
// 文字列
getline(cin,s);

とする。cinは改行文字を残すので、ignoreでそれを読み捨てないといけない。また、ignoreの引数は読み捨てる文字数。引数なしの場合は1を渡したのと同等の効果がある。

cout¶

有効数字等が設定されている問題は、必ず多めに出力すること。多めに出す分には大丈夫。

基本の使い方

int n;vector<int> V(n);
cout << n << endl;
rep(i,n) cout << V[i] << " ";
cout << endl;

以下主なiomanipの使い方

int n = 123;double d = 1.23;

//10進数 −> 123
cout << dec << n << endl ;

//8進数 −> 173
cout << oct << n << endl ;

//16進数 −> 7b
cout << hex << n << endl ;

//16進数かつ、大文字 −> 7B
cout << hex << uppercase << n << endl;

//10進数に戻す
cout << dec;

//幅が10になるようにする。デフォルトは右寄せ
// -> xxxxxxx123 (default)
cout << setfill('x') << setw(10) << right << n << endl;

// -> 123xxxxxxx
cout << setfill('x') << setw(10) << left << n << endl;
// -> 123yyyyyyy
cout << setfill('y') << setw(10) << n << endl;

//小数点以下10桁表示に。
cout << fixed << setprecision(10);

// -> 1.2300000000
cout << d << endl;
// -> 12.3000000000
cout << 10*d << endl;

//小数点の表示を元に戻す
cout.unsetf(ios::fixed);
// -> 1.23
cout << d << endl;

基本的には、引数のあるマニピュレータの効果は保存される。

高速化¶

以下のコードをmain関数の最初に書くことで、cin,coutの速度が2倍程度になる。

ios::sync_with_stdio(false);
cin.tie(0);

ただし、このコードはstdioとの同期を切るという意味なので、これを使うときにはprintfやscanfを使用してはだめ。

scanf,printf¶

C++では、cstdioをinclude。複雑な書式とかが必要なときにはこっちを使うといいかもしれない。

scanf¶

基本的な使い方

int n;char tmp[256];
scanf("%d\n",&n);
gets(tmp);

stringに直接いれるのはだめ。scanfはcinと同様に改行を残す。getlineするなら cin.ignore。getsするなら、直前のscanfで改行を読んでおく必要がある。また、scanfで改行を読むのでなく、直後にgetc(stdin)してもいい。

printf¶

基本的な使い方

int n = 100;
printf("n is %d\n",n);

scanfとほとんど同様の使い方ができる。

指定子	出力書式
%c	文字
%s	文字列
%d	符号付き10進整数
%u	符号なし10進
%f	10進浮動小数点数
%o	符号なし8進
%x	符号なし16進(Xなら大文字)
%%	%記号

std::vector¶

ここでは、配列のSTL版である、vectorの使いかたについて書く。ここに書かれている関数は、string等にも用いることができるものが多い。 vector<bool>は特殊化されているため注意が必要。また、all(vector)は、vector.begin(),vector.end()とdefineしている。

基本¶

//要素数10で、初期値は-1にする。
vector<int> vec(10,-1);
//vecの最初から3つの要素をコピーする。
vector<int> newvec(vec.begin(),vec.begin()+3);
//vecの最初から3つの要素を削除する。
vec.erase(vec.begin(),vec.begin()+3);

並び換え¶

sort¶

C++のsortは、 $O(n \log n)$ で、introsort。何も指定しない場合には昇順にソートされる。注意すべきなのは、C++11では、最悪ケースで$O(n log n)$となること。C++03では特に何も制限はないが、g++ならば$O(n log n)$である。

//昇順 (sort(vec.begin(),vec.end())) (2,1,3) -> (1,2,3)
sort(all(vec));
//降順 (ただ単純にreverseしてもいい) (2,1,3) -> (3,2,1)
sort(all(vec),greater<int>());

第三引数には関数を渡すことができる。

bool comp(const int& a ,const int& b){
    return abs(a) < abs(b);
}
int main(){
    vector<int> vec(10);
    //絶対値が小さい順にソート
    sort(all(vec),comp);
    sort(all(vec),[](int l,int r) {return abs(l)<abs(r)});
}

stable_sort¶

sortとちがって、同じ優先順位の要素の順番は保存される。最悪計算量は$O(n log ^ 2 n)$である。

stable_sort(all(vec),comp);

unique¶

隣あう同じ要素を一つにまとめる。eraseすることに注意。

int ints[] = {1,1,2,1,1};
vector<int> vec(ints,ints+5);
vec.erase(unique(all(vec)),vec.end());
// 1 2 1
repeat(i,vec.size()) cout << vec[i] << endl;

rotate¶

rotateは第二引数の場所を先頭にするように回転する。

int[] ints = {1,2,3,4,5};
vector<int> vec(ints,ints+5);
rotate(vec.begin(),vec.begin()+1),vec.end()); //2,3,4,5,1
rotate(vec.begin(),vec.end()-1,vec.end()); //5,1,2,3,4

next_permutation¶

順列をすべて列挙する。$N!$個なので、それなりの勢いで大きくなる。章末の付録参照。

vector<int> V = {3,2,1,4};
//ソートすること。
sort(all(V));

do{
   for(int i=0;i<V.size();i++){
       cout << V[i] << " ";
   }
   cout << endl;
}while(next_permutation(all(V)));

座標圧縮¶

// 一次元の座標圧縮を行う、 vector<T> をとって圧縮後の vector と元の値への　map (vector<T>) を返す
// O(N log N)
template<typename T>
std::pair<std::vector<size_t>, std::vector<T>> compress_vector(const std::vector<T>& v)
{
    // ソートし、重複を削除した配列の中で二分探索を行うことで何番目の大きさの要素かを求める
    std::vector<T> sorted_v = v;
    std::sort(sorted_v.begin(), sorted_v.end());
    sorted_v.erase(unique(sorted_v.begin(), sorted_v.end()), sorted_v.end());

    std::vector<size_t> compressed(v.size());
    std::vector<T> compress_to_original(sorted_v.size());

    for (size_t i = 0; i < v.size(); i++) {
        compressed[i] = std::lower_bound(sorted_v.begin(), sorted_v.end(), v[i]) - sorted_v.begin();
        compress_to_original[compressed[i]] = v[i];
    }
    return { std::move(compressed), std::move(compress_to_original) };
}

int my_main([[maybe_unused]] int argc, [[maybe_unused]] char** argv)
{
    vector<int> a = {8, 100, 33, 33, 33, 12, 6, 1211};
    auto [compressed, to_original] = compress_vector(a);

    // 圧縮後の結果と元の値を出力
    for(const auto& c : compressed) {
        std::cout << c << " => " << to_original[c] << std::endl;
    }
    return 0;
}

#include <iostream>
#include <vector>
#include <iomanip>
#include <algorithm>
#define rep(x,n) for(int x=0;x<(n);x++)
#define all(vec) vec.begin(),vec.end()


using namespace std;
template<typename T>
ostream& operator<<(ostream& os,const vector<T>& vec){
    os << "[";
    for(const auto& v : vec){
        os << v << ",";
    }
    os << "]";
    return os;
}

// 重なる長方形がカバーしている面積
int main(){
    int k = 0;
    while(true){
        k++;
        int n; cin >> n;
        if(n == 0) break;
        // (X1,Y1) is upper left.
        // (X1,Y1) is lower right.
        vector<double> X1(n),X2(n),Y1(n),Y2(n);
        rep(i,n){
            double x,y,r;
            cin >> x >> y >> r;
            X1[i] = x-r;
            X2[i] = x+r;
            Y1[i] = y-r;
            Y2[i] = y+r;
        }
        // X,Yは全てのX座標,Y座標
        vector<double> X,Y;
        X.insert(X.end(),all(X1));
        X.insert(X.end(),all(X2));
        Y.insert(Y.end(),all(Y1));
        Y.insert(Y.end(),all(Y2));

        // sortして重複を削除
        sort(all(X));
        sort(all(Y));
        X.erase(unique(all(X)),X.end());
        Y.erase(unique(all(Y)),Y.end());

        // X1,X2,Y1,Y2の値をX,Y中のindexにする。この時点でintにすべきかも。
        rep(i,n){
            X1[i] = lower_bound(all(X),X1[i]) - X.begin();
            X2[i] = lower_bound(all(X),X2[i]) - X.begin();
            Y1[i] = lower_bound(all(Y),Y1[i]) - Y.begin();
            Y2[i] = lower_bound(all(Y),Y2[i]) - Y.begin();
        }
        vector<vector<char>> is_on(X.size(),vector<char>(Y.size()));
        for(int i=0;i<n;i++){
            for(int sx=X1[i];sx<X2[i];sx++){
                for(int sy=Y1[i];sy<Y2[i];sy++){
                    is_on[sx][sy] = true;
                }
            }
        }
        double ans = 0;
        for(int i=0;i<X.size()-1;i++){
            for(int j=0;j<Y.size()-1;j++){
                if(is_on[i][j]){
                    ans += (X[i+1]-X[i]) * (Y[j+1]-Y[j]);
                }
            }
        }
        cout << fixed << setprecision(2);
        cout << k << " " << ans << endl;
    }
}

探索¶

全探索¶

全部しらべる。

int linear_search(vector<int> V,int val){
    repeat(i,V.size()){
        if(V[i] == val) return i;
    }
    return -1;
}

二分探索¶

単調なものに関してある条件を満たす点を探す。

template<typename T>
T binary_method(const std::function<bool(T)>& pred, T ok, T fail)
{
    static_assert(std::is_integral_v<T> || std::is_floating_point_v<T>);
    assert(pred(ok) == true);
    assert(pred(fail) == false);
    T EPS;
    if constexpr (std::is_integral_v<T>) EPS = 1;
    else if (std::is_floating_point_v<T>) EPS = 1e-9;

    while (abs(ok - fail) > EPS) {
        T mid = std::midpoint(ok, fail);
        if (pred(mid)) ok = mid;
        else fail = mid;
    }

    return ok;
}

対象がvectorの場合は以下のように書ける。 lower_bound, upper_bound を使える。

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

int main(){
    vector<int> v = {1,2,3,1,2,3,1,2,3};
    // ソートする必要あり。
    //  i  0 1 2 3 4 5 6 7 8
    // -> [1,1,1,2,2,2,3,3,3]
    sort(v.begin(),v.end());
    // 2以上の数値が初めて現れる場所 -> 3
    int lower = lower_bound(v.begin(),v.end(),2) - v.begin();
    // 2より大きい数値が初めて表われる場所 -> 6
    int upper = upper_bound(v.begin(),v.end(),2) - v.begin();
    // 2の個数
    int number_of_2 = upper - lower;
}

三分探索¶

凸関数の極値を求める

// 凸関数の極大な点をもとめる
template<typename T, typename F>
T ternary_search(function<F(T)> f, T left, T right)
{
    static_assert(std::is_integral_v<T> || std::is_floating_point_v<T>);
    assert(left <= right);
    // 整数の場合は候補を 3 つまでに絞る
    constexpr T EPS = (std::is_integral_v<T>) ? 2 : 1e-9;
    while (abs(left - right) > EPS) {
        T l = (2*left + right) / 3;
        T r = (left + 2*right) / 3;
        if (f(l) < f(r))
            left = l;
        else
            right = r;
    }
    if constexpr (std::is_integral_v<T>) {
        // 3 つの候補のうち一番大きいものを返す
        auto max_index = left;
        auto max_value = f(left);
        for (auto i = left + 1; i <= right; i++) {
            auto iv = f(i);
            if (iv > max_value) {
                max_value = iv;
                max_index = i;
            }
        }
        return max_index;
    } else
        return (left + right) / 2;
}

// 凹関数の極小な・を求める
template<typename F,typename T>
T ternary_search_concave(F f,T left,T right){
    return ternary_search([f](T x){return -f(x);},left,right);
}

黄金分割探索¶

三分探索と同様に凸関数の極大値を求めるが、分割点を使い回すことによって f の呼び出し回数が少ない。

// 黄金分割探索、凸関数の極大な点をもとめる
// f の呼び出し回数が少ないので、 f の呼び出しが重い時はこちらを使う
// この関数は定義域が実数の場合 (T が double など) のみ使える
// Ref: https://qiita.com/tanaka-a/items/f380257328da421c6584
template<typename T, typename F>
T golden_selection_search(function<F(T)> f, T left, T right)
{
    static_assert(std::is_floating_point_v<T>);
    assert(left <= right);

    double golden_ratio = (3 - std::sqrt(5)) / 2.0;
    constexpr T EPS = 1e-9;

    // 初期分割点
    T l = left + (right - left) * golden_ratio;
    T r = right - (right - left) * golden_ratio;
    T fl = f(l);
    T fr = f(r);
    while (abs(left - right) > EPS) {
        // 分割を更新する際に l, r を使い回す
        if (fl < fr) {
            left = l;
            l = r;
            fl = fr;
            r = right - (right - left) * golden_ratio;
            fr = f(r);
        } else {
            right = r;
            r = l;
            fr = fl;
            l = left + (right - left) * golden_ratio;
            fl = f(l);
        }
    }
    return (left + right) / 2;
}

フィボナッチ探索¶

黄金分割探索の定義域が整数版。

// フィボナッチ探索、凸関数の極大な点をもとめる
// f の呼び出し回数が少ないので、 f の呼び出しが重い時はこちらを使う (内部でメモする)
// この関数は定義域が整数の場合 (T が int など) のみ使える
// Ref: https://qiita.com/tanaka-a/items/f380257328da421c6584
// https://atcoder.jp/contests/typical90/tasks/typical90_ba
template<typename T, typename F>
T fibonacci_search(function<F(T)> f, T left, T right, const F MINF = std::numeric_limits<F>::min())
{
    static_assert(std::is_integral_v<T>);

    const T original_left = left;
    const T original_right = right;
    const T extended_left = left - 1;
    T extended_right = right;
    vector<T> fib;
    fib.push_back(1);
    fib.push_back(1);
    fib.push_back(2);
    fib.push_back(3);
    // 間の数をフィボナッチ数にしたいので、弊区間の長さはフィボナッチ数　+ 1 に調整
    while (fib.back() + 1 < (extended_right - extended_left + 1)) {
        fib.push_back(fib[fib.size() - 1] + fib[fib.size() - 2]);
    }
    extended_right = (fib.back() + 1) + extended_left - 1;
    assert(extended_right - extended_left + 1 == fib.back() + 1);

    map<T, F> memo; // f の呼び出しが遅い可能性があるのでメモしておく。
    auto extended_f = [&](T x) {
        if (memo.contains(x))
            return memo[x];
        if (x < original_left) {
            // 単調増加にするため、 extended_left で MINF にするように調整
            return memo[x] = MINF + (x - extended_left);
        } else if (original_right < x) {
            // 単調減少にするために extended_right で MINF にするように調整
            return memo[x] = MINF + (extended_right - x);
        } else {
            return memo[x] = f(x);
        }
    };

    left = extended_left;
    right = extended_right;

    // 初期分割点
    int cur_fib_index = fib.size() - 1;
    T l = left + fib[cur_fib_index - 2];
    T r = right - fib[cur_fib_index - 2];
    while (cur_fib_index >= 3) {
        // 分割を更新する際に l, r を使い回す
        --cur_fib_index;
        if (extended_f(l) < extended_f(r)) {
            left = l;
            l = r;
            r = right - fib[cur_fib_index - 2];
        } else {
            right = r;
            r = l;
            l = left + fib[cur_fib_index - 2];
        }
    }
    {
        // 最後には 3 つの要素が残っているはず。
        auto max_index = left;
        auto max_value = extended_f(max_index);
        for (auto i = left + 1; i <= right; ++i) {
            auto iv = extended_f(i);
            if (iv > max_value) {
                max_value = iv;
                max_index = i;
            }
        }
        return max_index;
    }
}

合計の重さXを達成するような組み合わせの数え上げ¶

#include <vector>
using namespace std;

// if we have objects that each has w_{i} weight.
//  we calc how many combination that has total weight X.
//  O(2^(N/2) logN)?
template<typename T>
int howmany_combination_such_that_total_weight_is_X(vector<T> W,int X){
    // sort(all(W));
    int N = W.size();
    vector<vector<T>> parts(2);
    copy(W.begin(),W.begin()+N/2,back_inserter(parts[0]));
    copy(W.begin()+N/2,W.end(),back_inserter(parts[1]));

    int ret = 0;
    vector<vector<T>> sums(2,{0});
    for(size_t p=0;p<parts.size();p++){
        for(size_t i=0;i<parts[p].size();i++){
            vector<T> nex;
            for(size_t j=0;j<sums[p].size();j++){
                nex.push_back(sums[p][j]);
                nex.push_back(sums[p][j]+parts[p][i]);
            }
            sums[p].swap(nex);
        }
    }
    for(vector<T>& sum : sums){
        sort(sum.begin(),sum.end());
    }

    for(T s : sums[0]){
        ret += upper_bound(sums[1].begin(),sums[1].end(),X-s)
             - lower_bound(sums[1].begin(),sums[1].end(),X-s);
    }
    return ret;
}

ぐるぐる¶

いつかつかうかも。 http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1189

// WH/2,WH/2 を中心としたぐるぐる。
// まんなかは1
vector<vector<int>> guruguru(int N){
    const int WH = sqrt(N) + 2;
    vector<vector<int>> field(WH,vector<int>(WH,N + 100));
    {
        int cx=WH/2,cy=WH/2;
        field[cy][cx] = 1;
        int c = 2;
        for(int i=0;;i++){
            int sx = cx + i+1;
            int sy = cy + i;
            for(int j=0;j<2*i+2;j++){
                field[sy-j][sx] = c;
                c++;
                if(c >= MAXN) goto end;
            }
            vector<int> dx = {-1,0,1};
            vector<int> dy = {0,1,0};
            // left
            int ny = sy-2*i-1;
            int nx = sx;
            for(int k=0;k<3;k++){
                for(int j=0;j<2*(i+1);j++){
                    ny = ny + dy[k];
                    nx = nx + dx[k];
                    field[ny][nx] = c;
                    c++;
                    if(c >= MAXN) goto end;
                }
            }
        }
end:;
    }
    return field;
}

文字列操作¶

stringをincludeする。cctypeもいるかも。

std::string¶

charをラップしたテンプレートクラス。基本的な使い方について

部分列¶

          //0123456789
string str("abcdefghij");
// 5番目以降
str.substr(5);    // "fghij"
// 5番目から3つ
str.substr(5,3); // "fgh"
//全部小文字にする
transform(s.begin(),s.end(),s.begin(),::tolower);

substrは一つの引数でそこから先全部、二つの場合は第一引数の位置から、第二引数の数だけ持ってくる。

検索¶

stringには、いくつかのfindが定義されている。線形検索なので、早い検索が必要なときには後述するKMP法やBM法を用いること。

find 引数が最初に現れる位置
rfind 引数が最後に表われる位置
find_first_of 引数の文字のどれかが最初に表われる位置
find_last_of 引数の文字のどれかが最後に表われる位置
find_first_not_of 引数の文字のどれかではない文字が最初に表われる位置
find_first_not_of 引数の文字のどれかではない文字が最後に表われる位置

第二引数として探すための最初の位置を指定できる。

Boyer Moore法¶

TBD

KMP法¶

TBD

stringstream¶

cinやcoutのようなstreamをstringを元に作成したりする。基本的には、string をフォーマットしたり、intやlongに、intやlongから変換するために使用する。

stringstream ss;
ss << 102;
string s;
ss >> s;

再帰下降構文解析¶

BNFを書いて、それにしたがっていく。左再帰除去を忘れずに。

#include <iostream>
#include <string>
#include <cctype>
#include <sstream>

using namespace std;

typedef string::const_iterator Cursor;
class ParseError{};

// <四則演算の式> ::= <乗算除算の式> (+ or -) <乗算除算の式> (+ or -) ...
// <乗算除算の式> ::= <括弧か数> (* or /) <括弧か数> (* or /) ...
// <括弧か数>     ::= '(' <四則演算の式> ')' or <数>
// <数>           ::= (0|1|2|3|4|5|6|7|8|9)+

int expression(Cursor&);
int term(Cursor&);
int factor(Cursor&);
int number(Cursor&);

// <四則演算の式> ::= <乗算除算の式> (+ or -) <乗算除算の式> (+ or -) ...
int expression(Cursor &c){
    int ret = term(c);
    while(*c == '+' or *c == '-'){
        if(*c == '+'){
            c++;
            ret += term(c);
        }else{
            c++;
            ret -= term(c);
        }
    }
    return ret;
}

// <乗算除算の式> ::= <括弧か数> (* or /) <括弧か数> (* or /) ...
int term(Cursor &c){
    int ret = factor(c);
    while(*c == '*' or *c == '/'){
        if(*c == '*'){
            c++;
            ret *= factor(c);
        }else{
            c++;
            ret /= factor(c);
        }
    }
    return ret;
}

// <括弧か数>     ::= '(' <四則演算の式> ')' or <数>
int factor(Cursor &c){
    if(*c == '('){
        c++;
        int ret = expression(c);
        c++; // ')'
        return ret;
    }else{
        return number(c);
    }
}

// <数>           ::= (0|1|2|3|4|5|6|7|8|9)+
int number(Cursor &c){
    stringstream ss;
    while(isdigit(*c)){
        ss << *c;
        c++;
    }
    int ret;
    ss >> ret;
    return ret;
}

int main(){
    int test_case;
    cin >> test_case;
    cin.ignore();
    for(int i=0;i<test_case;i++){
        string s;
        getline(cin,s);
        Cursor c = s.begin();
        cout << expression(c) << endl;
    }
    return 0;
}

LL1文法(== 再帰下降で書けるかどうか)¶

予測型構文解析表を埋めて重複をチュックする。テストしたい。一般に、"X -> Xb,X -> a"なる生成規則は、"X->aX',X'->bX',X'->"とできる。これは前者がab*なる形の文字列を導出するからである。

#include <iostream>
#include <vector>
#include <set>
#include <map>

using namespace std;

struct Rule{
    char from;
    vector<char> to;
    Rule(char from_,vector<char> to_)
        : from(from_),to(to_){
    }
};

bool operator==(const Rule& lhs, const Rule& rhs){
    if(lhs.from != rhs.from or lhs.to.size() != rhs.to.size()){
        return false;
    }
    for(size_t i=0;i<lhs.to.size();i++){
        if(lhs.to[i] != rhs.to[i]){
            return false;
        }
    }
    return true;
}

bool operator<(const Rule& lhs, const Rule& rhs){
    if(lhs.from != rhs.from){
        return lhs.from < rhs.from;
    }
    if(lhs.to.size() != rhs.to.size()){
        return lhs.to.size() < rhs.to.size();
    }
    for(size_t i=0;i<lhs.to.size();i++){
        if(lhs.to[i] != rhs.to[i]){
            return lhs.to[i] < rhs.to[i];
        }
    }
    return false;
}

bool operator>(const Rule& lhs, const Rule& rhs){
    return not(lhs == rhs) and not(lhs < rhs);
}

ostream& operator<<(ostream& os,Rule r){
    os << r.from << "->";
    for(char c : r.to){
        os << " " << c;
    }
    return os;
}

template<typename T>
ostream& operator<<(ostream& os,set<T> s){
    os << "{";
    for(T v : s){
        os << v << ",";
    }
    os << "}";
    return os;
}

// can be written in recursive descent parsing?
//  see tiger book(p.46).
bool can_be_written_in_recursive_descent_parsing(const set<char>& terminal,
                                                 const set<char>& non_terminal,
                                                 const vector<Rule>& rules){
    for(Rule r : rules){
        cerr << r << endl;
    }
    map<char,bool> nullable;
    map<char,set<char>> first,follow;
    for(char c : non_terminal){
        nullable[c] = false;
        first[c] = {};
        follow[c] = {};
    }
    for(char c : terminal){
        nullable[c] = false;
        first[c] = {c};
        follow = {};
    }

    {
        bool changed;
        do{
            changed = false;
            for(Rule r : rules){
                bool all_nullable = true;
                for(char t : r.to){
                    all_nullable = all_nullable and nullable[t];
                }
                if(all_nullable and not nullable[r.from]){
                    nullable[r.from] = true;
                    changed = true;
                }
                for(size_t i=0;i<r.to.size();i++){
                    bool nul = true;
                    for(size_t k=0;k<i;k++){
                        nul = nul and nullable[r.to[k]];
                    }
                    if(nul){
                        size_t bef = first[r.from].size();
                        first[r.from].insert(first[r.to[i]].begin(),
                                             first[r.to[i]].end());
                        changed = changed or (bef != first[r.from].size());
                    }
                }
                for(size_t i=0;i<r.to.size();i++){
                    // nullable[y[k]] for k <- [i+1,]
                    bool nul = true;
                    for(int k=i+1;k<(int)r.to.size();k++){
                        nul = nul and nullable[r.to[k]];
                    }
                    if(nul){
                        size_t bef = follow[r.to[i]].size();
                        follow[r.to[i]].insert(follow[r.from].begin(),
                                               follow[r.from].end());
                        changed = changed or (bef != follow[r.to[i]].size());
                    }
                }
                for(size_t i=0;i<r.to.size();i++){
                    for(size_t j=i+1;j<r.to.size();j++){
                        // nullable[y[k]] for k <- [i+1,j-1]
                        bool nul = true;
                        for(size_t k=i+1;k<j;k++){
                            nul = nul and nullable[r.to[k]];
                        }
                        if(nul){
                            size_t bef = follow[r.to[i]].size();
                            follow[r.to[i]].insert(first[r.to[j]].begin(),
                                                   first[r.to[j]].end());
                            changed = changed or (bef != follow[r.to[i]].size());
                        }
                    }
                }
            }
        }while(changed);
    }
    for(char c : non_terminal){
        cerr << c << endl
             << "nullable? " << nullable[c] << endl
             << "first " << first[c] << endl
             << "follow " << follow[c] << endl;
    }
    map<char,map<char,set<Rule>>> maps;
    for(Rule r : rules){
        // check r.to is nullable or not.
        //  and collect first.
        set<char> to_first;
        bool to_nullable = true;
        for(char t : r.to){
            to_first.insert(first[t].begin(),first[t].end());
            to_nullable = to_nullable and nullable[t];
            if(not to_nullable) break;
        }

        for(char tf : to_first){
            maps[r.from][tf].insert(r);
        }
        if(to_nullable){
            for(char t : follow[r.from]){
                maps[r.from][t].insert(r);
            }
        }
    }

    cerr << "Map" << endl;
    bool dup = false;
    for(const pair<char,map<char,set<Rule>>> m : maps){
        int from = m.first;
        for(const pair<char,set<Rule>> p : m.second){
            int to = p.first;
            cerr << char(from) << " " << char(to) << " " << p.second << endl;
            dup = dup or p.second.size() > 1;
        }
    }
    return not dup;
}


int main(){
    set<char> terminal;
    set<char> non_terminal;

    {
        // terminal;
        int n;cin >> n;
        for(int i=0;i<n;i++){
            char c;cin >> c;
            terminal.insert(c);
        }
    }

    {
        // non_terminal;
        int n;cin >> n;
        for(int i=0;i<n;i++){
            char c;cin >> c;
            non_terminal.insert(c);
        }
    }

    // number of rule.
    int n;
    cin >> n;
    vector<Rule> rules;
    for(int i=0;i<n;i++){
        int k;cin >> k;
        char from;cin >> from;
        vector<char> to(k);
        for(int j=0;j<k;j++){
            cin >> to[j];
        }
        rules.push_back(Rule(from,to));
    }
    cout << can_be_written_in_recursive_descent_parsing(terminal,non_terminal,rules) << endl;
    return 0;
}

              // 9                  // 9
// 3          // $ + - * / i n ( )  // $ + - * / i n ( )
// a c d      // 3                  // 6
// 3          // S E T              // S E T D U F
// X Y Z      // 10                 // 12
// 6          // 2 S E $            // 2 S E $
// 1 X Y      // 3 E E p T          // 2 E T D
// 1 X a      // 3 E E m T          // 3 D + T D
// 0 Y        // 1 E T              // 3 D - T D
// 1 Y c      // 3 T T * F          // 0 D
// 1 Z d      // 3 T T / F          // 2 T F U
// 3 Z X Y Z  // 1 T F              // 3 U * F U
//  -> 0      // 1 F i              // 3 U / F U
              // 1 F n              // 0 U
              // 3 F ( E )          // 1 F i
              //  -> 0              // 1 F n
                                    // 3 F ( E )
                                    //  -> 1

Z algorithm¶

文字列のすべての部分文字列に対し、それと全体の共通接頭文字列の長さのテーブルを作る。 O(N)

// z-algorithm (O(|S|)
// Return V where V[i] = length of the largest common prefix of S and S[i]
// Example:
//  S: aaabaaaab
//  V: 921034210
// Ref: https://ivanyu.me/blog/2013/10/15/z-algorithm/
vector<int> zAlgorithm(const std::string& s)
{
    std::vector<int> z(s.size());
    z[0] = s.size();

    // left-bound, right-bound of current z-box.
    int lt = 0, rt = 0;

    for (int k = 1; k < s.size(); k++) {
        if (k > rt) {
            // If k is outside the current z-box, then we do naive computation.
            int n = 0;
            while (n + k < s.size() && s[n] == s[n + k])
                n++;
            z[k] = n;
            if (n > 0) {
                lt = k;
                rt = k + n - 1;
            }
        } else {
            // If k is inside the current z-box, then consider two cases
            int p = k - lt; // pair index
            // rightPartLen is a length of k..rt 
            int rightPartLen = rt - k + 1;

            if (z[p] < rightPartLen) {
                z[k] = z[p];
            } else {
                // s[i - k] is a A[i] where A starts from k of S.
                int i = rt + 1;
                while (i < s.size() && s[i] == s[i - k])
                    i++;
                z[k] = i - k;
                lt = k;
                rt = i - 1;
            }
        }
    }
    return z;
}

Rolling Hash¶

// TODO: string 以外にも対応する
struct RollingHash {
    // ハッシュを外から与えるコンストラクタ、複数の文字列でハッシュを共有したい時はこちら
    // RollingHash(t, other_rollinghash.r()) のように使う
    RollingHash(std::string s, const uint64_t exp_r)
    : array(std::move(s)),
    length(array.size()),
    exp_r(exp_r)
    {
        calc_pows();
        calc_hash();
    }
    // ハッシュを自動的に決定するコンストラクタ、一つの文字列内で部分文字列のハッシュを比較する場合はこちら
    explicit RollingHash(std::string s)
        : RollingHash(std::move(s), decide_exp_r()) {}

    // 半開区間 [l, r) のハッシュ値を求める
    uint64_t hash(int l, int r)
    {
        assert(l <= r);
        assert(0 <= l && l <= static_cast<int>(length));
        assert(0 <= r && r <= static_cast<int>(length));
        // return pre_hash[r] - pre_hash[l];
        uint64_t sub = (pows[r-l] * pre_hash[l]) % MOD;
        if (pre_hash[r] >= sub) {
            return (pre_hash[r] - sub) % MOD;
        } else {
            return (pre_hash[r] + MOD - sub) % MOD;
        }
    }
    // [l1, r1) == [l2, r) かどうかチェックする
    // ハッシュしかチェックしないので外す可能性もある。
    uint64_t probably_equal(int l1, int r1, int l2, int r2)
    {
        assert(l1 <= r1 && l2 <= r2);
        if (r1 - l1 != r2 - l2) return false; // 長さが違ったら false
        if (l1 == l2 && r1 == r2) return true;
        if (hash(l1, r1) != hash(l2, r2)) return false;
        return true;
    }
    // [l1, r1) == [l2, r) かどうかチェックする
    // ハッシュが一致した時に完全な文字列比較を行う
    uint64_t exactly_qual(int l1, int r1, int l2, int r2)
    {
        if (!probably_equal(l1, r1, l2, r2)) return false;
        for (int i = 0; i < (r1-l1); i++) {
            if (array[l1+i] != array[l2+i]) return false;
        }
        return true;
    }
    uint64_t r() { return exp_r; }
private:
    void calc_pows()
    {
        pows.resize(length+1);
        pows[0] = 1;
        for (size_t i = 1; i <= length; i++) {
            pows[i] = (pows[i-1] * exp_r) % MOD;
        }
    }
    void calc_hash()
    {
        pre_hash.resize(length+1);
        pre_hash[0] = 0;
        for (size_t i = 1; i <= length; i++) {
            pre_hash[i] = ((exp_r * pre_hash[i-1]) % MOD + to_num(array[i-1])) % MOD;
        }
    }
    uint64_t to_num(char c)
    {
        return c;
    }
    // ハックを防ぐために exp_r は乱数にしておく
    static uint64_t decide_exp_r()
    {
        std::random_device seed_gen;
        std::mt19937_64 mt64(seed_gen());
        std::uniform_int_distribution<uint64_t> get_rand_uni_int(2, MOD-2);
        return get_rand_uni_int(mt64);
    }
    static constexpr uint64_t MOD = 1000000007;
    static_assert(MOD < (1ull << 32));
    // hash[i] = [0..i) の hash = array[0] * r^(i-1) + array[1] * r^(i-2) + ... array[i-1] * r^0.
    std::vector<uint64_t> pre_hash;
    // pows[i] = r^i
    std::vector<uint64_t> pows;
    // ハッシュ衝突時チェック用にコピーしておく
    std::string array;
    size_t length;
    uint64_t exp_r;
};

// https://atcoder.jp/contests/tessoku-book/tasks/tessoku_book_bd
int my_main([[maybe_unused]] int argc, [[maybe_unused]] char** argv)
{
    int n, q; cin >> n >> q;
    string s; cin >> s;

    std::array<RollingHash, 2> hash = { RollingHash(s), RollingHash(s)};

    rep(_, q)
    {
        int a, b, c, d; cin >> a >> b >> c >> d;
        a--; b--; c--; d--;
        bool all_match = true;
        repv(h, hash) all_match &= h.probably_equal(a, b+1, c, d+1);
        cout << (all_match ? "Yes" : "No") << endl;
    }
    return 0;
}

Suffix Array¶

// 接尾辞をソートして保持する
// 空文字列を含む
// 実装は Manber & Myers の O(n logn^2) アルゴリズム
// 例: missisippi
// sa[0]  = 10 φ
// sa[1]  = 9 i
// sa[2]  = 6 ippi
// sa[3]  = 4 isippi
// sa[4]  = 1 issisippi
// sa[5]  = 0 missisippi
// sa[6]  = 8 pi
// sa[7]  = 7 ppi
// sa[8]  = 5 sippi
// sa[9]  = 3 sisippi
// sa[10] = 2 ssisippi
struct SuffixArray {
    explicit SuffixArray(std::string s)
        : m_string(std::move(s)),
          m_size(m_string.size()),
          m_rank(m_size+1),
          m_prefix_index(m_size+1)
    {
        const int n = static_cast<int>(m_size);
        m_rank[n] = -1; // n のケースは空文字列扱いで一番早いことにする
        m_prefix_index[n] = n;
        for (int i = 0; i < n; i++) {
            m_prefix_index[i] = i;
            m_rank[i] = static_cast<int>(m_string[i]); // 文字コードをランクにする
        }
        // k 文字のランクが計算されていたとして、　k*2 文字のランクを計算する
        for (int k = 1; k <= n; k *= 2) {
            // (rank[i], rank[i+k]), (rank[j], rank[j+k]) から k*2 文字の比較を行う関数
            auto compare = [this, k, n](int i, int j) {
                // 最初の k 文字ですでに順序がつく場合はそれを採用
                if (m_rank[i] != m_rank[j]) return m_rank[i] < m_rank[j];
                // 次の k 文字で比較する。もしオーバーしている (空文字列) ならそれを優先する。
                const int ri = (i+k) <= n ? m_rank[i+k] : -1;
                const int rj = (j+k) <= n ? m_rank[j+k] : -1;
                return ri < rj;
            };
            std::sort(m_prefix_index.begin(), m_prefix_index.end(), compare);
            std::vector<int> tmp_rank(n+1);
            // 空文字列のランクは 0 にして
            tmp_rank[m_prefix_index[0]] = 0;
            // それを基準に一致していたら +0、文字列が一致していなかったら +1 していく。
            // suffix_array はソート済みのため、 compare が false を返すケースは一致しているケース
            for (int i = 1; i <= n; i++) {
                tmp_rank[m_prefix_index[i]] = tmp_rank[m_prefix_index[i-1]] + (compare(m_prefix_index[i-1], m_prefix_index[i]) ? 1 : 0);
            }
            swap(m_rank, tmp_rank); // m_rank に中身を移す
        }
    }
    // s が文字列の中に何度現れるかを返す。 O(logn)
    // https://atcoder.jp/contests/abc362/tasks/abc362_g
    int count(const std::string& s)
    {
        if (s.size() == 0) return 0;
        return upper_bound(s) - lower_bound(s);
    }
    // 接尾時の先頭 s.size() 文字を見て、辞書式順で s 以上のものが初めて現れる場所を返す
    int lower_bound(const std::string& s)
    {
        if (s.size() == 0) return 0;
        // fail = 0 は空文字列であり、必ず s より辞書順で小さい。
        // succ = m_size+1
        int fail = 0, succ = m_size+1;
        while (abs(fail - succ) > 1) {
            int mid = (fail + succ) / 2;
            if (s <= suffix(mid).substr(0, min(s.size(), suffix(mid).size()))) succ = mid;
            else fail = mid;
        }

        return succ;
    }
    // 接尾時の先頭 s.size() 文字を見て、辞書式順で s より大きいものが初めて現れる場所を返す
    int upper_bound(const std::string& s)
    {
        if (s.size() == 0) return 1;
        // fail = 0 は空文字列であり、必ず s より辞書順で小さい。
        // succ = m_size+1
        int fail = 0, succ = m_size+1;
        while (abs(fail - succ) > 1) {
            int mid = (fail + succ) / 2;
            if (s < suffix(mid).substr(0, min(s.size(), suffix(mid).size()))) succ = mid;
            else fail = mid;
        }

        return succ;
    }
    [[nodiscard]] size_t size() const { return m_size+1; }
    size_t operator[](int i) const { return m_prefix_index[i]; }
    [[nodiscard]] std::string_view suffix(int i) const { return std::string_view(m_string).substr(m_prefix_index[i]); }
    std::vector<int>& ref() { return m_prefix_index; }

private:
    // 念の為コピーしておく
    std::string m_string;
    size_t m_size;
    std::vector<int> m_rank;
    // その文字列が何文字目から始まっているか
    std::vector<int> m_prefix_index;
};

ST = TS の条件¶

文字列 S と t が与えられる。長さ t の文字列 T であって、 ST = TS を満たすものを出力せよ。ただし T がない可能性もある。参考: https://atcoder.jp/contests/arc181/tasks/arc181_b

|S| >= 1 かつ |T| >= 1 の時、 ST = TS ならば g = gcd(|S|, |T|) として S および T は g の長さのパターンで構成される。

例

|S| = 2, |T| = 3 なら gcd(2, 3) = 1 なので |S| および |T| は 1 文字の繰り返しで構成される。 S = aa T = a など。
|S| = 4, |T| = 2 なら gcd(2, 4) = 2 なので |S| および |T| は 2 文字の繰り返しで構成される。 S = abab, T = ab など。

// 長さ t の文字列 T であり、 ST = TS を満たす文字列を返す。
// なかったら nullopt
// 例:
// "A", 1 -> "A",
// "A", 4 -> "AAAA"
// "ab", 1 -> nullopt
// "ab", 2 -> "ab"
// "ab", 3 -> nullopt
// "ab", 4 -> "abab"
// "abcabc", 3 -> "abc"
// "abcabc", 2 -> nullopt
std::optional<std::string> solve_ST_eq_TS(const std::string& s, size_t t)
{
    // ST = TS なる S および T は必ず gcd(|S|, |T|) の長さのパターンの繰り返しである
    int g = std::gcd(s.size(), t);
    // もしパターンの条件を S が満たさないなら nullopt
    for (int i = 0; i + g < s.size(); i++) {
        if (s[i] != s[i+g]) return std::nullopt;
    }
    const std::string pattern = s.substr(0, g);
    std::string ret;
    while (ret.size() < t) ret += pattern;
    return ret;
}

説明: |S| <= |T| とする。この時 T の先頭 |S| 文字は S と一致するので、 T = S + A 書くことができる。

ST = SSA
TS = SAS

以上で SA = AS なる似た問題が出てくる。もし |A| >= |S| なら同様に A = S + B とかくことができる。

ST = SSA = SSSB
TS = SAS = SSBS

この手続きを |S| > |C| となるまで繰り返す。

ST = SSSSSSC
TS = SSSSSCS

さらに SC = CS という似た問題で S = C + D とかく、ことを繰り返す。

以上の変形を数式で表すことにすると T = q * S + C であった。 q は |T| / |S|、 |C| は　|T| % |S| である。この時 |S| と |T| の最大公約数を g とし |X| と |Y| が互いに素な文字列 X, Y, Z があると

T = g * X
S = g * Y
C = g * Z (※)

と表すことができる。つまり T および S は長さ g = gcd(|S|, |T|) のパターンの繰り返しである。

(※) |T| - |S| = g|X| - g|Y| = g(|X| - |Y|) となるが、この時 (|X| - |Y|) と |X| は互いに素であるため |T| および |S| と、 |T| - |S| は最大公約数 g を持つ。これを繰り返すと　|T| - |S| - |S| = g (x - y - y) は |T| - |S| と最大公約数 g を持つので、結局 |T| - q |S| = |C| は |S| および |T| と最大公約数　g を持つ。

以上の計算はユークリッドの互除法と同じことをしている。

行列¶

// T 型を要素とする行列型
// Matrix<int> m = {{1, 2}, {3, 4}}
template<typename T>
struct Matrix {
    Matrix(size_t height, size_t width) : m_height(height), m_width(width), m_value(height, vector<T>(width)) { }
    Matrix(const vector<vector<T>>& v) :  m_height(v.size()), m_width(v[0].size()), m_value(v) { } // NOLINT
    Matrix(initializer_list<initializer_list<T>> init) : m_height(init.size()), m_width(init.begin()->size()), m_value(m_height, vector<T>(m_width))
    {
        size_t r = 0;
        for (const auto& row : init) {
            size_t c = 0;
            for (const auto& col : row) {
                m_value[r][c] = col;
                c++;
            }
            assert(c == m_width);
            r++;
        }
    }
    // i 行目を取り出す
    vector<T>& operator[](size_t i) { return m_value[i]; }
    const vector<T>& operator[](size_t i) const { return m_value[i]; }

    Matrix& operator+=(const Matrix& rhs)
    {
        assert(m_width == rhs.m_width);
        assert(m_height == rhs.m_height);
        for (size_t i = 0; i < m_height; i++)
            for (size_t j = 0; j < m_width; j++)
                m_value[i][j] += rhs.m_value[i][j];
        return *this;
    }
    Matrix& operator-=(const Matrix& rhs)
    {
        assert(m_width == rhs.m_width);
        assert(m_height == rhs.m_height);
        for (size_t i = 0; i < m_height; i++)
            for (size_t j = 0; j < m_width; j++)
                m_value[i][j] -= rhs.m_value[i][j];
        return *this;
    }
    friend Matrix operator+(Matrix lhs, const Matrix& rhs) { lhs += rhs; return lhs;}
    friend Matrix operator-(Matrix lhs, const Matrix& rhs) { lhs -= rhs; return lhs;}
    friend Matrix operator*(T lhs, Matrix rhs)
    {
        for (size_t i = 0; i < rhs.m_height; i++)
            for (size_t j = 0; j < rhs.m_width; j++)
                rhs.m_value[i][j] *= lhs;
        return rhs;
    }
    // AxB の行列と BxC の行列を受け取って、 AxC の行列を返す
    friend Matrix operator*(Matrix lhs, const Matrix& rhs)
    {
        assert(lhs.m_width == rhs.m_height);
        // TODO: 最適化する。 r -> k -> c の方が早いはず？？
        Matrix ret(lhs.m_height, rhs.m_width);
        for (size_t r = 0; r < lhs.m_height; r++)
            for (size_t c = 0; c < rhs.m_width; c++)
                for (size_t k = 0; k < lhs.m_width; k++)
                    ret.m_value[r][c] += lhs.m_value[r][k] * rhs.m_value[k][c];
        return ret;
    }
    // 一次元ベクトルとの掛け算
    friend vector<T> operator*(Matrix lhs, const vector<T>& rhs)
    {
        assert(lhs.m_width == rhs.size());
        vector<T> ret(lhs.m_height);
        for (size_t r = 0; r < lhs.m_height; r++) {
            for (size_t c = 0; c < lhs.m_width; c++) {
                ret[r] += lhs[r][c] * rhs[c];
            }
        }
        return ret;
    }
    // 行列を n 乗
    Matrix pow(uint64_t n)
    {
        assert(m_height == m_width);
        if (n == 0) {
            Matrix E(m_height, m_width);
            for (size_t i = 0; i < m_height; i++) E[i][i] = 1;
            return E;
        }
        Matrix ret = (*this * *this).pow(n/2);
        if (n % 2 == 1) ret = ret * *this;
        return ret;
    }
    friend ostream& operator<<(ostream& os, const Matrix& rhs)
    {
        for (size_t i = 0; i < rhs.m_height; i++) {
            for (size_t j = 0; j < rhs.m_width; j++) {
                if (j != 0) os << " ";
                os << rhs.m_value[i][j];
            }
            if (i != rhs.m_height - 1) os << endl;
        }
        return os;
    }
    // 行基本変形によって上三角化 (標準形) への変形を行う
    // https://drken1215.hatenablog.com/entry/2019/03/20/202800
    int triangler_inplace(bool is_extended = false)
    {
        assert(m_width >= 1 && m_height >= 1);
        int rank = 0;
        // 連立方程式を解く際には係数拡大行列を作る、その際にはその部分は掃き出し法の対象にしない
        for (size_t col = 0; col < (is_extended ? (m_width-1) : m_width) ; col++) {
            std::optional<int> pivot = std::nullopt;
            for (size_t row = rank; row < m_height; row++) {
                if constexpr (std::is_floating_point_v<T>) {
                    if ((!pivot.has_value() && abs(m_value[row][col]) > 1e-10) || (pivot.has_value() && abs(m_value[row][col]) > abs(m_value[pivot.value()][col])))
                        pivot = row;;
                } else {
                    if (!pivot.has_value() && m_value[row][col] != 0) {
                        pivot = row;
                        break;
                    }
                }
            }
            // ピボットが見つからなかったら次の列へ
            if (!pivot.has_value()) continue;
            swap(m_value[pivot.value()], m_value[rank]);

            // ピボットの値を 1 にする
            for (int j = m_width - 1; j >= static_cast<int>(col); j--) {
                m_value[rank][j] /= m_value[rank][col];
            }

            // 掃き出しによって他の列から
            for (size_t row = 0; row < m_height; row++) {
                if (static_cast<int>(row) == rank) continue;
                auto m = m_value[row][col];
                if constexpr (std::is_floating_point_v<T>) {
                    if (abs(m) < 1e-10) continue;
                } else {
                    if (m == 0) continue;
                }
                for (size_t j = rank; j < m_width; j++) {
                    m_value[row][j] -= m_value[rank][j] * m;
                }
            }
            rank++;
        }

        return rank;
    }
    // 行基本変形によって上三角化を行う
    Matrix trianglar(bool is_extended = false)
    {
        Matrix ret = *this;
        ret.triangler_inplace(is_extended);
        return ret;
    }
    [[nodiscard]] size_t width() const { return m_width; }
    [[nodiscard]] size_t height() const { return m_height; }
private:
    size_t m_height;
    size_t m_width;
    vector<vector<T>> m_value;
};

Gauss-Jordan¶

掃き出し法によって連立一次方程式を解く (https://drken1215.hatenablog.com/entry/2019/03/20/202800)

// Gauss-Jordan 法によって一次連立方程式を解く
template <typename T>
struct GaussJordanResult {
    enum class SolutionType {
        none, // 解なし
        one, // 一意に定まる
        some  // いくつか存在する
    };
    SolutionType type;
    Matrix<T> extended_matrix {0, 0};
    int rank;
};

// mat x = v を解く
// result.type == none なら解なし (一次方程式の中に矛盾が存在)
// result.type == one なら解が一つだけ
// result.type == some の時、実際の解の数は T に依存する
//   自由変数 (パラメータ) の数が mat.width() - rank 個存在するため、
//   - もし T が mod_int なら解の数は 2 ** (mat.width() - rank)
//   - もし T が実数 (double など) なら無限
template<typename T>
GaussJordanResult<T> solve_linear_equations(const Matrix<T>& mat, const vector<T>& v)
{
    assert(mat.height() == v.size());
    // 拡大係数行列を作成する (右側に v をくっつける)
    GaussJordanResult<T> result;

    result.extended_matrix = Matrix<T>(mat.height(), mat.width() + 1);

    size_t extended_col = mat.width();
    for (size_t r = 0; r < mat.height(); r++) {
        for (size_t c = 0; c < mat.width(); c++) {
            result.extended_matrix[r][c] = mat[r][c];
        }
        result.extended_matrix[r][extended_col] = v[r];
    }

    result.rank = result.extended_matrix.triangler_inplace(true);

    for (size_t row = result.rank; row < mat.height(); row++) {
        // 矛盾があるので解なし
        bool no_answer = false;
        if constexpr (std::is_floating_point_v<T>)
            no_answer = abs(result.extended_matrix[row][extended_col]) > 1e-10;
        else
            no_answer = result.extended_matrix[row][extended_col] != 0;

        if (no_answer) {
            result.type = GaussJordanResult<T>::SolutionType::none;
            return result;
        }
    }
    if (result.rank == static_cast<int>(result.extended_matrix.width() - 1)) {
        result.type = GaussJordanResult<T>::SolutionType::one;
    } else {
        // some の時、実際の解の数は T に依存する
        // 自由変数 (パラメータ) の数が mat.width() - rank 個存在するため、
        // - もし T が mod_int なら解の数は 2 ** (mat.width() - rank)
        // - もし T が実数 (double など) なら無限
        result.type = GaussJordanResult<T>::SolutionType::some;
    }
    return result;
}

// https://atcoder.jp/contests/typical90/tasks/typical90_be
int main([[maybe_unused]] int argc, [[maybe_unused]] char** argv)
{
    int n, m; cin >> n >> m;

    Matrix<mod_int<2>> mat(m, n);
    for (int i = 0; i < n; i++) {
        int t; cin >> t;
        for (int j = 0; j < t; j++) {
            int a; cin >> a;
            a--;
            mat[a][i] = 1;
        }
    }

    vector<mod_int<2>> v(m);
    for (int i = 0; i < m; i++) {
        int s; cin >> s;
        v[i] = s;
    }

    auto r = solve_linear_equations(mat, v);

    if (r.type == GaussJordanResult<mod_int<2>>::SolutionType::one) {
        cout << 1 << endl;
        return 0;
    } else if (r.type == GaussJordanResult<mod_int<2>>::SolutionType::some) {
        mod_int<998244353> ret = 2;
        ret = powi(ret, n - r.rank);
        cout << ret << endl;
    } else {
        cout << 0 << endl;
    }

    return 0;
}

データ構造¶

Union-Find¶

グループを管理するデータ構造、グループの統合や同じグループにいるかの判定が高速に可能

struct UnionFind {
    explicit UnionFind(int n)
    : data(n, -1)
    , number_of_group_var(n) { }
    // 親を探す
    int root(int x) {
        if (data[x] < 0) return x;
        return data[x] = root(data[x]);
    }
    // x,yの含まれる集合を併合する
    // 併合が発生したら true, しなかったら false
    bool unite(int x,int y) {
        x = root(x);
        y = root(y);
        if (x == y) return false; // すでに同じところにいる
        // 大きいほうに追加する。
        if(data[y] < data[x]) swap(x,y);
        data[x] += data[y];
        data[y] = x;
        number_of_group_var--;
        return true;
    }
    // 同じ集合にいるかどうか
    bool is_same_group(int x, int y) {
        return root(x) == root(y);
    }
    [[nodiscard]] int number_of_group() const {
        return number_of_group_var;
    }
    int group_size(int x) {
        return -data[root(x)];
    }
private:
    // data は親ならはそのグループのサイズに -1 をかけたもの、子なら親の値
    vector<int> data;
    int number_of_group_var;
};

ヒープ¶

#include <queue>
#include <iostream>

using namespace std;

typedef pair<int,int> pii;

struct Comp{
    bool operator()(pii left,pii right){
        if(left.second < right.second) return true;
        else if(left.second == right.second and left.first > right.first) return true;
        else return false;
    };
};

struct Robot{
    int y,x,dir,step;
    Robot(int y,int x,int dir,int step) : y(y),x(x),dir(dir),step(step) {};
};

// <,>を定義すればless<Robot>みたいに扱える。
bool operator<(const Robot& lhs,const Robot& rhs){
    return lhs.step < rhs.step;
}
bool operator>(const Robot& lhs,const Robot& rhs){
    return lhs.step > rhs.step;
}

int main(){
    // 何も書かないと降順。(おっきい方からでてくる。)
    // これは昇順(ちいさいほうから出てくる)にしたもの。
    priority_queue<int,vector<int>,greater<int> > Qi;
    //関数オブジェクトを使っていい感じにもできる。
    priority_queue<pii,vector<pii>,Comp> Q;
    // 自作クラスの場合はこんな感じ
    priority_queue<Robot,vector<Robot>,greater<Robot> > que;

    Q.push(make_pair(1,2));
    Q.push(make_pair(2,2));
    Q.push(make_pair(3,2));
    while(not Q.empty()){
        cout << Q.top().first << " " << Q.top().second << endl;
        Q.pop();
    }
}

bitset¶

限られた大きさのvector<bool>を使いたいときに、bitsetを使うことができる。

#include <iostream>
#include <bitset>

using namespace std;
const int N = 10;

struct bit_cmp{
    bool operator() (const bitset<N> &left,const bitset<N> &right) {
        for(int i=N-1;i>=0;i--){
            if(left[i] < right[i]) return true;
            if(left[i] > right[i]) return false;
        }
        return false;
    }
};


int main(){
    //定数じゃないとダメ。最初は全部false
    bitset<N> bits;
    // すべての要素をtrue -> 1111111111
    bits.set();
    if(bits.all()) cout << "all" << endl;
    // 立っているbitの数 -> 10
    cout << bits.count() << endl;
    // すべての要素をfalse -> 0000000000
    bits.reset();
    if(bits.none()) cout << "none" << endl;

    //1番目の要素をtrue -> 0100000000
    bits.set(1);
    if(bits.any()) cout << "any" << endl;

    // 0110000000
    bits.set(2);
    //1番目の要素をfalse -> 0010000000
    bits.reset(1);

    if(bits[2]) cout << 2 << endl;
    cout << bits << endl;

    bitset<N> newbits;
    // 和を取る
    bits |= newbits;
    // 積を取る
    bits &= newbits;

    // 関数オブジェクトを作る必要アリ
    map<bitset<N>,int,bit_cmp> M;
}

分数¶

テストちゅう.

typedef long long ll;
ll gcd(ll a,ll b){
    return b==0 ? a : gcd(b,a%b);
}

ll lcm(ll a,ll b){
    if(a < 0) a *= -1;
    if(b < 0) b *= -1;
    return a*b / gcd(a,b);
}

struct Fraction{
    ll n,d;
    Fraction(ll _n,ll _d){
        ll c = lcm(_n,_d);
        n = c / _d;
        d = c / _n;
        if(d < 0){
            n *= -1;
            d *= -1;
        }
    }

    Fraction(ll _n){
        Fraction(_n,1);
    }

    bool operator<(const Fraction& r) const{
        ll c_d = lcm(d,r.d);
        return n*(c_d/d)< r.n*(c_d/r.d);
    }
    bool operator>(const Fraction& r) const{
        return not ((*this) < r or (*this) == r);
    }
    bool operator<=(const Fraction& r) const{
        return (*this) < r or (*this) == r;
    }
    bool operator>=(const Fraction& r) const{
        return (*this) > r or (*this) == r;
    }
    bool operator==(const Fraction& r) const{
        return n == r.n and d == r.d;
    }
    Fraction operator+(const Fraction& r) const{
        ll c_d = lcm(d,r.d);
        return Fraction(n*(c_d/d)+r.n*(c_d/r.d),c_d);
    }
    Fraction operator-(const Fraction& r) const{
        return (*this) + (-r);
    }
    Fraction operator*(const Fraction& r) const{
        return Fraction(n*r.n,d*r.d);
    }
    Fraction operator/(const Fraction& r) const{
        return (*this) * Fraction(r.d,r.n);
    }
    Fraction operator+() const{
        return Fraction(n,d);
    }
    Fraction operator-() const{
        return (*this) * -1;
    }
    Fraction operator*(const ll& a) const{
        return Fraction(a*n,d) ;
    }
};

ostream& operator<<(ostream &os,const Fraction& f){
    os << f.n << "/" << f.d;
    return os;
}

セグメント木¶

普通の¶

RMQ (POJ 3264)

template<typename T, auto op, auto e>
class SegmentTree {
    static_assert(std::is_convertible_v<decltype(op), std::function<T(T, T)>>,
              "op must work as T(T, T)");
    static_assert(std::is_convertible_v<decltype(e), std::function<T()>>,
              "e should be T()");
public:
    explicit SegmentTree(int _n) : original_n(_n)
    {
        n = 1;
        while (n < _n) n *= 2;
        data = vector<T>(2 * n - 1, e());
    }
    template<class InputIterator>
    explicit SegmentTree(InputIterator first, InputIterator last)
    {
        n = 1;
        original_n = distance(first, last);
        while (n < original_n) n *= 2;
        data = vector<T>(2 * n - 1, e());
        size_t i = 0;
        for (auto it = first; it != last; ++it, i++) {
            update(i, *it);
        }
    }
    void update(int i, T v)
    {
        assert(0 <= i && i < original_n);
        int k = i + n - 1;
        data[k] = v;
        while (k > 0) {
            k = get_parent_index(k);
            auto left = data[get_left_index(k)];
            auto right = data[get_right_index(k)];
            data[k] = op(left, right);
        }
    }

    // [a, b) 半開区間の値を取得する、 a == b なら e() を返す
    [[nodiscard]] T get(int a, int b) const
    {
        assert(0 <= a && a <= b && b <= original_n);
        // ノード 0 の管理区間は [0, n)
        return get(a, b, 0, 0, n);
    }
    [[nodiscard]] size_t size() const { return n; }
private:
    static int get_left_index(int k) { return (k << 1) + 1; }
    static int get_right_index(int k) { return (k << 1) + 2; };
    static int get_parent_index(int k) { return (k - 1) >> 1; }
    // l, r はノード k の管理する区間
    [[nodiscard]] T get(int a, int b, int k, int l, int r) const
    {
        if (r <= a || b <= l) return e();
        if (a <= l && r <= b) return data[k];
        const int left_index = get_left_index(k);
        const int right_index = get_right_index(k);
        const int mid = (l+r) >> 1;
        const auto left = get(a, b, left_index, l, mid);
        const auto right = get(a, b, right_index, mid, r);
        return op(left, right);
    }

    int n = 0;
    int original_n = 0;
    vector<T> data;
};


int main(int argc, char** argv) {
    // RMQ.
    const auto op = [](int a, int b) {return max(a, b); };
    const auto e = []() { return 0; };
    vector<int> a;
    SegmentTree<int, op, e> st(all(a));
}

Lazy Segment Tree¶

// T = 各要素および区間取得の結果の型
// op = T(T, T), 区間取得を行う関数 (例えば RMQ なら min)
// e = T(), op における単位元
// F = "Lazy" 遅延している操作の型
// mapping = T(F, T) 遅延している操作 F を区間取得の結果 T に対し行う関数
// composition = F(F f, F g)、遅延している操作 g に対し、追加で f を行った時に生成する遅延操作を生成する
//   f(g(x)) となるので g がさきに行われる処理
// id = F(), composition の恒等写像 (遅延している操作の単位元)
// F(op(x, y)) = op(F(x), F(y)) を満たす必要がある
template<typename T, auto op, auto e, typename F, auto mapping, auto composition, auto id>
class LazySegmentTree {
    static_assert(std::is_convertible_v<decltype(op), std::function<T(T, T)>>, "op should be T(T, T)");
    static_assert(std::is_convertible_v<decltype(e), std::function<T()>>, "e should be T()");
    static_assert(std::is_convertible_v<decltype(mapping), std::function<T(F, T)>>, "mapping should be T(F, T)");
    static_assert(std::is_convertible_v<decltype(composition), std::function<F(F, F)>>, "composition should be F(F, F)");
    static_assert(std::is_convertible_v<decltype(id), std::function<F()>>, "id should be F()");

public:
    explicit LazySegmentTree(int _n)
    {
        n = 1;
        while (n < _n) n *= 2;
        data = vector<T>(2 * n - 1, e());
        lazy = vector<F>(2 * n - 1, id());

    }
    template<class InputIterator>
    explicit LazySegmentTree(InputIterator first, InputIterator last) : LazySegmentTree(distance(first, last))
    {
        size_t i = 0;
        for (auto it = first; it != last; ++it, i++) {
            update(i, *it);
        }
    }
    void update(int i, T v)
    {
        int k = i + n - 1;
        data[k] = v;
        while (k > 0) {
            k = get_parent_index(k);
            auto left = data[get_left_index(k)];
            auto righ = data[get_right_index(k)];
            data[k] = op(left, righ);
        }
    }
    // [a, b) に操作 f を行う
    void apply(int a, int b, F f)
    {
        assert(0 <= a && a <= b && b <= n);
        return apply(a, b, f, 0, 0, n);
    }

    // [a, b) 半開区間の値を取得する
    [[nodiscard]] T get(int a, int b)
    {
        assert(0 <= a && a < b && b <= n); // FIXME: should check the original n?
        // ノード 0 の管理区間は [0, n)
        return get(a, b, 0, 0, n);
    }
    [[nodiscard]] size_t size() const { return n; }
private:
    static int get_left_index(int k) { return 2 * k + 1; }
    static int get_right_index(int k) { return 2 * k + 2; };
    static int get_parent_index(int k) { return (k - 1) / 2; }
    // k 番目のノードの遅延操作を反映する
    void apply_lazy(int k, int l, int r)
    {
        if (lazy[k] == id()) return;
        data[k] = mapping(lazy[k], data[k]);
        // 子供に遅延評価を伝搬する
        if (r - l > 1) {
            const int left_index = get_left_index(k);
            const int right_index = get_right_index(k);
            lazy[left_index] = composition(lazy[k], lazy[left_index]);
            lazy[right_index] = composition(lazy[k], lazy[right_index]);
        }
        lazy[k] = id();
    }

    // l, r はノード k の管理する区間
    [[nodiscard]] T get(int a, int b, int k, int l, int r)
    {
        if (r <= a || b <= l) return e();
        apply_lazy(k, l, r);
        if (a <= l && r <= b) return data[k];
        const int left_index = get_left_index(k);
        const int right_index = get_right_index(k);
        const auto left = get(a, b, left_index, l, (l+r)/2);
        const auto right = get(a, b, right_index, (l+r)/2, r);
        return op(left, right);
    }
    void apply(int a, int b, F f, int k, int l, int r)
    {
        apply_lazy(k, l, r);
        if (r <= a || b <= l) return;
        if (a <= l && r <= b) {
            lazy[k] = composition(f, lazy[k]);
            apply_lazy(k, l, r);
        } else {
            apply(a, b, f, get_left_index(k), l, (l+r)/2);
            apply(a, b, f, get_right_index(k), (l+r)/2, r);
            data[k] = op(data[get_left_index(k)], data[get_right_index(k)]);
        }
    }

    int n = 0;
    vector<T> data;
    vector<F> lazy;
};

範囲セット、範囲の和の計算をするセグメント木¶

// https://atcoder.jp/contests/abc371/tasks/abc371_f
struct Group {
    long long value; // ノードの合計値
    long long size; // ノードに含まれる要素の数
};

ostream& operator<<(ostream& os, const Group& g)
{
    return os << "{" << "," << g.value << "," << g.size << "}";
}

using T = Group;
auto op = [](const T& a, const T& b) {
    return T {a.value + b.value, a.size + b.size};
};
auto e = []() { return T {0, 0}; };
using F = optional<long long>;

auto mapping = [](const F& f, const T& g) {
    if (!f.has_value()) return g;
    // ノードに対するアップデートであるため、その範囲の合計値計算にサイズが必要
    return T {f.value() * g.size, g.size };
};
auto composition = [](const F& f, const F& g) {
    return f;
};

auto id = []() -> F { return std::nullopt; };
LazySegmentTree<T, op, e, F, mapping, composition, id> lst(n);

// 初期化はサイズ 1 の要素アップデートとして行う
for (int i = 0; i < n; i++) lst.update(i, Group {x[i], 1});

範囲への足し算、範囲の和の計算をするセグメント木¶

// https://atcoder.jp/contests/abc371/tasks/abc371_f
struct Group {
    long long value; // ノードの合計値
    long long size; // ノードに含まれる要素の数
};

using T = Group;
auto op = [](const T& a, const T& b) {
    return T {a.value + b.value, a.size + b.size};
};
auto e = []() { return T {0, 0}; };
using F = optional<long long>;

auto mapping = [](const F& f, const T& g) {
    if (!f.has_value()) return g;
    // ノードに対するアップデートであるため、その範囲の合計値計算にサイズが必要
    return T {f.value() * g.size + g.value, g.size };
};
auto composition = [](const F& f, const F& g) {
    if (g.has_value())
        return f.value() + g.value();
    return f.value();
};
auto id = []() -> F { return std::nullopt; };
LazySegmentTree<T, op, e, F, mapping, composition, id> lst(n);

// 初期化はサイズ 1 の要素アップデートとして行う
for (int i = 0; i < n; i++) lst.update(i, Group {x[i], 1});

Binary-Indexed Tree (FenwickTree)¶

オペレーションとしては SegmentTree の下位互換だが、高速な構造。

// SegmentTree はモノイドであれば良かったが、それに加えて可換であることと逆元 inv が必要になる。
// 例: 値への加算と部分集合の和を計算する
// const auto op = [](int a, int b) {return a+b; };
// const auto e = []() { return 0; };
// const auto inv = [](int a) { return -a;};
// FenwickTree<int, op, e, inv> ft(5);
template<typename T, auto op, auto e, auto inv>
class FenwickTree {
    static_assert(std::is_convertible_v<decltype(op), std::function<T(T, T)>>, "op must work as T(T, T)");
    static_assert(std::is_convertible_v<decltype(e), std::function<T()>>, "e should be T()");
    static_assert(std::is_convertible_v<decltype(inv), std::function<T(T)>>, "inv must work as T(T)");
public:
    explicit FenwickTree(int _n) : original_n(_n) { data = vector<T>(_n + 1); }

    // a[i] に　v を加算 (op) する。 (a[i] = op(a[i], v)
    void add(int i, T v)
    {
        assert(0 <= i && i < original_n);
        i++; // 1-indexed に変換
        for (size_t j = i; j < data.size(); j += (j & -j)) {
            data[j] = op(data[j], v);
        }
    }
    // [a, b) 半開区間の値を取得する、 a == b なら e() を返す
    [[nodiscard]] T get(int a, int b) const
    {
        assert(0 <= a && a <= b && b <= original_n);
        return op(get(b), inv(get(a)));
    }
    // [0, b) 半開区間の値を取得する、 b == original_n なら　e() を返す
    [[nodiscard]] T get(int b) const
    {
        assert(b <= original_n);
        int a = 1; // a を 1-indexed に。
        // b++, b--; b を 1-indexed に、 さらに閉区間に。
        T ret = e();
        for (int j = b; j > 0; j -= (j & -j)) {
            ret = op(ret, data[j]);
        }
        return ret;
    }
    vector<T> data;
    size_t original_n;
};

定数個のみを保持する priority_queue¶

ビームサーチとかに使える？

#include <vector>
#include <queue>
#include <iostream>


template<class T,
         class Container=std::vector<T>,
         class Compare=std::less<typename Container::value_type>>
struct lens_queue : public std::priority_queue<T,Container,Compare>{
    typedef std::priority_queue<T,Container,Compare> super;
    lens_queue(int _cap)
        : cap(_cap){
    }
    inline std::size_t capacity() const {return cap;}
    void push(const T& x){
        if(should_push(x)){
            super::push(x);erase_overflow();
        }
    }
    void push(T&& x){
        if(should_push(x)){
            super::push(std::move(x));erase_overflow();
        }
    }
    template <class... Args>
    void emplace(Args&&... args){
        this->c.emplace_back(std::forward<Args>(args)...);
        if(should_push(this->c.back())){
            std::push_heap(this->c.begin(),this->c.end(),this->comp);
            erase_overflow();
        }else{
            this->c.pop_back();
        }
    }
    // this crash lens_queue.
    Container&& take(){
        return std::move(this->c);
    }
    template<class T2,class Container2,class Compare2>
    friend std::ostream& operator<<(std::ostream&,
                                    const lens_queue<T2,Container2,Compare2>&);
private:
    std::size_t cap;
    void erase_overflow(){
        while(this->size() > capacity()){
            this->pop();
        }
    }
    inline bool should_push(const T& x){
        return this->size()<capacity() || !this->comp(this->top(),x);
    }
};

template<class T,class Container,class Compare>
std::ostream& operator<<(std::ostream& os,const lens_queue<T,Container,Compare>& que){
    os << "{";
    for(const auto& v : que.c){
        os << v << ",";
    }
    os << "}";
    return os;
}



signed main() {
    int n,k;char c;
    std::cin >> n >> k >> c;
    if(c == 'g'){
        // k = 3
        // 4 2 3 1 5 -> 3 4 5
        lens_queue<int,std::vector<int>,std::greater<int>> que(k);
        for(int i=0;i<n;i++){
            int v;std::cin >> v;
            que.push(v);
            // que.emplace(v);
        }
        vector<int> taken = que.take();
        // while(not que.empty()){
        //     std::cout << que.top() << std::endl;
        //     que.pop();
        // }
    }else if(c == 'l'){
        // k = 3
        // 4 2 3 1 5 -> 3 2 1
        lens_queue<int> que(k);
        for(int i=0;i<n;i++){
            int v;std::cin >> v;
            que.push(v);
            // que.emplace(v);
        }
        while(not que.empty()){
            std::cout << que.top() << std::endl;
            que.pop();
        }
    }
    return 0;
}

Set で区間を管理するテクニック¶

// 「区間を set で管理するテクニック」
// 半開区間 [a, b) を pair(a, b) として set にいれる。
// ただし set の中の区間は重複を持たないよう管理する
// https://codeforces.com/contest/915/problem/E
struct RangeSet {
    using T = int64_t;
    RangeSet()
    {
        // 番兵
        segs.emplace(std::numeric_limits<T>::min(), std::numeric_limits<T>::min());
        segs.emplace(std::numeric_limits<T>::max(), std::numeric_limits<T>::max());
    }
    // [l, r) を完全に含む区間が存在するならその区間を返す
    std::optional<std::pair<T, T>> containedBy(T l, T r)
    {
        // assert(l <= r);
        // [l, r) が含まれているとしたら一つ前のもの
        auto it = std::upper_bound(segs.begin(), segs.end(), make_pair(l, r));
        --it;
        if (it->first <= l && r <= it->second) return *it;
        return std::nullopt;
    }
    // 区間の数を返す
    size_t size()
    {
        return segs.size() - 2; // 番兵のぶんをひく
    }
    // 区間 [l, r) を追加する。区間のうち、新規に追加された要素を返す
    T insert(T l, T r)
    {
        // assert(l <= r);
        auto it = segs.upper_bound(make_pair(l, r));
        --it;
        // カバーされているので 0
        if (it->first <= l && r <= it->second) return 0;
        // 以下では l, r を書き換えていくが、この時既存の要素でカバーされていたものを "erased" とする
        // そして最後に追加された要素を (l-r) - erased で計算。
        T erased = 0;
        // 左の区間と一部被っている場合は合体する
        if (it->first <= l && l <= it->second) {
            erased += it->second - it->first;
            l = it->first;
            it = segs.erase(it);
        } else {
            ++it; // 左の区間とは被っていないので次へ
        }
        // 右の区間と共通部分があれば合体していく (it = upper_bound 相当)
        while (it->first <= r) {
            if (it->second <= r) {
                // 新しくできた要素で完全にカバーされている
                erased += it->second - it->first;
                it = segs.erase(it);
            } else {
                // l <= it->first <= r < it->second
                erased += it->second - it->first;
                r = it->second;
                segs.erase(it);
                // 既存の要素の方が r が長いのであればそれ以上に長い要素は存在しない
                break;
            }
        }
        segs.emplace(l, r);
        return (r - l) - erased;
    }
    // 区間 [l, r) を削除する。区間のうち、削除された要素を返す
    T erase(T l, T r)
    {
        auto it = segs.upper_bound(make_pair(l, r));
        --it;
        // 完全にカバーされている場合はその区間を分割する
        if (it->first <= l && r <= it->second) {
            if (it->first != l)
                segs.emplace(it->first, l);
            if (r != it->second)
                segs.emplace(r, it->second);
            segs.erase(it);
            return r - l;
        }
        T erased = 0;
        // 左の区間と一部被っている場合は左の区間の削除を行う
        if (it->first <= l && l < it->second) {
            erased += it->second - l;
            if (it->first != l) {
                segs.emplace(it->first, l);
            }
            it = segs.erase(it);
        } else {
            ++it;
        }
        // 右の区間と共通部分があれば削除していく (it = upper_bound 相当)
        while (it->first < r) {
            if (it->second <= r) {
                // 完全に削除される
                erased += it->second - it->first;
                it = segs.erase(it);
            } else {
                // l <= it->first <= r < it->second
                erased += r - it->first;
                segs.emplace(r, it->second);
                segs.erase(it);
                break;
            }
        }
        return erased;
    }
    // x 以上のカバーされていない最小の値を返す
    T mex(T x)
    {
        auto it = segs.lower_bound(make_pair(x, x+1));
        if (it->first <= x && x <= it->second) return it->second;
        return x;
    }
private:
    std::set<std::pair<T, T>> segs;
    friend ostream& operator<<(ostream&, const RangeSet&);
};
ostream& operator<<(ostream& os, const RangeSet& v) {
    bool first = true;
    os << "{";
    for (auto it = v.segs.begin(); it != v.segs.end(); ++it) {
        // 最初と最後は番兵なので skip
        if (it == v.segs.begin() || it == (std::prev(v.segs.end()))) continue;
        if (!first) {
            os << ", ";
        }
        os << *it;
        first = false;
    }
    return os << "}";
}

ZobristHash¶

集合をハッシュで管理する構造。要素を hash 化して集合をその hash の xor で扱う。 (xor で扱う場合は集合に重複する要素を許さない。許したい場合は演算を + で行う)

なお、アップデートが必要であったりする場合は vector に関する問題なら SegmentTree を使うこともできる。 (コード参照)

// 集合の一致判定を確率的に行う
// ローリングハッシュのように vector を受け取って、その一部分の集合のハッシュを計算する
// ハッシュに xor を使う場合は多重集合が扱えないため、 関数を + にする必要がある。
// 指定するハッシュ関数は T によって変わる。 T が文字列などの場合は std::hash も使えるが、 long long などの場合は自分で定義する必要がある (乱数を割り当てれば良い)
template<typename T, bool MULTISET=false>
struct ZobristHash {
    explicit ZobristHash(const std::vector<T>& arr, const std::function<uint64_t(T)>& myhash) : hashes(arr.size()), hashes_acc(arr.size() + 1)
    {
        for (size_t i = 0; i < arr.size(); i++) {
            hashes[i] = myhash(arr[i]);
        }
#if DEBUG
        {
            // MULTISET ではないなら hash は (できるだけ) 全て違う必要がある。
            if (!MULTISET) {
                assert(std::set(all(hashes)).size() == arr.size());
            }
        }
#endif
        hashes_acc[0] = zero();
        for (size_t i = 1; i <= hashes.size(); i++) {
            hashes_acc[i] = op(hashes_acc[i-1], hashes[i-1]);
        }
    }
    // [l, r) のハッシュ値を返す
    [[nodiscard]] uint64_t hash(const size_t l, const size_t r) const
    {
        assert(l <= r);
        return op2(hashes_acc[r], hashes_acc[l]);
    }
    static uint64_t op(const uint64_t l, const uint64_t r)
    {
        if (MULTISET) return l + r;
        else          return l ^ r;
    }
    static uint64_t op2(const uint64_t l, const uint64_t r)
    {
        if (MULTISET) return l - r;
        else return l ^ r;
    }
    static uint64_t zero()
    {
        return 0;
    }
    std::vector<uint64_t> hashes;
    // i 番目の要素は 集合 [0, i) のハッシュ
    std::vector<uint64_t> hashes_acc;
};
template<typename T>
using MultiZobristHash = ZobristHash<T, true>;

int main(int argc, char **argv) {
    auto a = input_vector<int>(n);
    auto b = input_vector<int>(n);

    // 対応表になる unordered_map を作り乱数を入れておく
    std::random_device seed_gen;
    std::mt19937_64 engine(seed_gen());
    std::uniform_int_distribution<uint64_t> distribution;
    std::unordered_map<int, uint64_t> hash_map;
    auto myhash = [&hash_map, &distribution, &engine](int x) -> uint64_t {
        if (hash_map.contains(x)) return hash_map[x];
        return hash_map[x] = distribution(engine);
    };
    MultiZobristHash<int> ah(a, myhash);
    MultiZobristHash<int> bh(b, myhash);

    bool probably_same = ah.hash(la, ra+1) == bh.hash(lb, rb+1);

    // セグメントツリーを使う版。重複を許さない場合は + ではなく xor を使う。
    auto op = [](uint64_t left, uint64_t right) { return left + right; };
    auto e = []() -> uint64_t { return 0;};
    SegmentTree<uint64_t, op, e> seg_a(n);
    for (int i = 0; i < n; i++) {
        seg_a.update(i, myhash(a[i]));
    }
    SegmentTree<uint64_t, op, e> seg_b(n);
    for (int i = 0; i < n; i++) {
        seg_b.update(i, myhash(b[i]));
    }
    bool probably_same = seg_a.get(la, ra+1) == seg_b.get(lb, rb+1);
    // セグメントツリー使う版ここまで。

}

動的計画法的なやつ¶

Longest Common Sequence (LCS)¶

O(NM)

string lcs(const string& a,const string& b){
    const int n = a.size(),m = b.size();
    vector<vector<int> > x(n+1,vector<int>(m+1));
    vector<vector<int> > y(n+1,vector<int>(m+1));

    for(int i=0;i<n;i++){
        for(int j=0;j<m;j++){
            if(a[i] == b[j]){
                x[i+1][j+1] = x[i][j]+1;
                y[i+1][j+1] = 0;
            }else if(x[i+1][j] < x[i][j+1]){
                x[i+1][j+1] = x[i][j+1];
                y[i+1][j+1] = 1;
            }else{
                x[i+1][j+1] = x[i+1][j];
                y[i+1][j+1] = -1;
            }
        }
    }
    string ret;
    for(int i=n,j=m;i>0 and j>0;){
        if(y[i][j] > 0){
            i--;
        }else if(y[i][j] < 0){
            j--;
        }else{
            ret.push_back(a[i-1]);
            i--;j--;
        }
    }
    reverse(all(ret));
    return ret;
}

Longest Increasing Subsequence (LIS)¶

与えられた数列の最長増加増加部分列 (部分列であって、その要素が単調に増加するもの) を計算する。

// https://atcoder.jp/contests/abc354/editorial/10027
// O(n log n) で最長増加部分列を求める。 strict = true なら同じ値は許さない、 false なら同じ値も許す
template<typename T>
std::pair<size_t, std::vector<size_t>> longest_increasing_sequence(const std::vector<T>& v, const bool strict=true) {
    static_assert(std::is_signed_v<T> , "T should be signed value"); // 符号付きだと dp 表の番兵と被る可能性があるので許さないことにしておく
#ifdef DEBUG
    for (const auto& t : v) assert(t != std::numeric_limits<T>::min());
#endif

    // 長さの計算と一緒に、その要素がどの長さの場合を更新したかを覚えていく
    std::vector<size_t> updated_index(v.size());
    // dp[x] = x の長さを実現した際の数列の最後尾の値
    std::vector<T> dp(v.size()+1, std::numeric_limits<T>::max());
    dp[0] = std::numeric_limits<T>::min();
    for (int i = 0; i < v.size(); i++) {
        auto it = strict ? std::lower_bound(all(dp), v[i]) : std::upper_bound(all(dp), v[i]);
        *it = std::min(v[i], *it);
        updated_index[i] = it - dp.begin();
    }
    size_t length = std::lower_bound(all(dp), std::numeric_limits<T>::max()) - dp.begin() - 1;
    return {length, updated_index };
}

void lis_example()
{
    std::vector<int> v = {1, 3, 4, 2, 1, 5};
    auto [length, updated_index] = longest_increasing_sequence(v, false);
    // 最長増加部分列を実現する index
    vector<int> lis_index;
    {
        int cur_len = length;
        for (int i = static_cast<int>(updated_index.size()) -1; i >= 0; i--) {
            if (cur_len == updated_index[i]) {
                lis_index.push_back(i);
                cur_len--;
            }
        }
        reverse(all(lis_index));
    }
}

巡回セールスマン問題¶

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DPL_2_A

ある点からスタートして、全ての点を辿って最初の点に戻る時の最短距離を計算する。

// dp[a][b] = 現在 b にいて、 集合 a をすでに辿っている時の最短距離
vector<vector<int>> dp(1 << v, vector<int>(v, INF));
// スタート地点にいて、かつ辿った集合なしの状態を距離 0 にする。
// この時辿った集合に s を入れないことによって最後に s に戻ってくるようにする。
dp[0][s] = 0;
for (uint32_t g = 0; g < (1 << v); g++) {
    // dp[g][p] から配る DP で遷移
    for (int p = 0; p < v; p++) {
        for (const auto& e : graph[p]) {
            if (g & (1 << e.to)) continue; // すでに遷移済みならスキップ
            dp[g | (1 << e.to)][e.to] = min(dp[g | (1 << e.to)][e.to], dp[g][p] + e.value);
        }
    }
}
int ret = dp[(1 << v)-1][s];

部分集合の列挙¶

https://atcoder.jp/contests/typical90/tasks/typical90_as より。

bit で集合を表すことにしたとき、その空集合を除く部分集合の列挙は以下のような for 文で行うことができる。

uint32_t g = 0b0110111;
for (uint32_t a = g; a != 0; a = (a - 1) & g) {
    for (int i = 0; i < 32; i++) {
        cout << ((a & (1 << i)) ? 1 : 0);
    }
    cout << endl;
}

これを利用して、部分集合の部分集合を列挙して行って bit DP を行うことがある。以下は全体を k 個の部分集合に分割するDP。 O(k*3^n)。

vector<vector<ll>> dp(k + 1, vector<ll>(1 << n, 1ll << 62));
dp[0][0] = 0;

for (int cnt = 1; cnt <= k; cnt++) {
    for (int g = 1; g < (1 << n); g++) {
        for (int a = g; a != 0; a = (a - 1) & g) {
            chmin(dp[cnt][g], max(dp[cnt-1][g-a], max_dist[a]));
        }
    }
}

ナップサック問題¶

個数制限なしのとき

#include <iostream>
#include <vector>

using namespace std;

#define rep(i,n) for(int i=0;i<(int)n;i++)

typedef vector<int> vi;

template<class T>
inline int len(const T& t){
    return t.size();
}

struct Treasure{
    int value,weight;
    Treasure():
        value(0),weight(0) {};
    Treasure(int _value,int _weight)
        : value(_value),weight(_weight) {}
};

int main(){
    int W,N;
    cin >> W;
    cin >> N;
    vector<Treasure> v(N);
    repeat(i,N){
        cin >> v[i].value >> v[i].weight;
    }
    // もしも複数入れて良いなら、直接valueを更新する。
    vi value(W+1);
    repeat(n,N){
        vi tmp(W+1);
        repeat(i,W+1){
            tmp[i] = max(tmp[i],value[i]);
            int in = i+v[n].weight;
            if(in <= W){
                tmp[in] = max(tmp[in],value[i]+v[n].value);
            }
        }
        value.swap(tmp);
    }

    int retw=0;
    int retv=0;
    repeat(i,len(value)){
        if(value[i] > retv){
            retv = value[i];
            retw = i;
        }
    }
    cout << retv << endl << retw << endl;
    return 0;
}

桁 DP¶

N 以下の整数で、かつ桁ごとや bit ごとに処理していけば良いような DP を考える。

例: 桁の合計が 3 の倍数である数
例: 立っている bit が n 個である数

例えば 10 進数で計算しているとき、とある桁が動く範囲は 0 ~ 9 もしくは 0 ~ N のある桁の数になる (N を超えないようにするため) なので一般に以下のような DP を考える。

DP[何桁目まで決めた][すでに絶対 N より小さいか] = 通り数

DP[i+1][true] <- dp[i][1] から i 桁目を自由に決める
DP[i+1][false] <- dp[i][false] から i 桁目を N と同じ
DP[i+1][true] <- dp[i][false] から i 桁目を 0 ~ (N の i 桁目 - 1)

https://drken1215.hatenablog.com/entry/2019/02/04/013700

以下が上から bit を決めていく DP

// https://atcoder.jp/contests/abc117/tasks/abc117_d
ll n, k; cin >> n >> k;
auto a = input_vector<ll>(n);

int MAX_BIT = 62;
// dp[i][smaller]
vector<vector<ll>> dp(MAX_BIT+1, vector<ll>(2, -1));
dp[0][0] = 0;
for (int i = 0; i < MAX_BIT; i++) {
    ll fil = 1ll << (MAX_BIT - i - 1);
    int select0 = 0, select1 =  0;
    for (ll ia : a) {
        if (ia & fil) {
            select0++;
        } else {
            select1++;
        }
    }
    // smaller -> smaller
    if (dp[i][1] != -1) {
        chmax(dp[i+1][1], dp[i][1] + fil * select0);
        chmax(dp[i+1][1], dp[i][1] + fil * select1);
    }
    if (dp[i][0] != -1) {
        // !smaller -> !smaller
        if (k & fil) {
            chmax(dp[i+1][0], dp[i][0] + fil * select1);
        } else {
            chmax(dp[i+1][0], dp[i][0] + fil * select0);
        }
        // !smaller -> smaller
        if (k & fil) {
            chmax(dp[i+1][1], dp[i][0] + fil * select0);
        }
    }
}

cout << *max_element(all(dp.back())) << endl;

ヒストグラム内の最大長方形のサイズ¶

#define repeat(i,n) for(int i=0;i<static_cast<int>(n);i++)
typedef vector<int> vi;
typedef vector<vi> vvi;

// see http://algorithms.blog55.fc2.com/blog-entry-132.html
//  max rectangle area in histogram.
template<typename T>
T max_rectangle_size_in_histogram(const vector<T>& v){
    // where,height. in this function,use index as where.
    typedef pair<int,T> piT;
    stack<piT> sta;
    T ret = 0;
    for(int i=0;i<v.size();i++){
        if(sta.empty() or sta.top().second < v[i]){
            sta.push(make_pair(i,v[i]));
        }else if(sta.top().second == v[i]){
            // pass
        }else if(sta.top().second > v[i]){
            piT t;
            while(not sta.empty() and sta.top().second >= v[i]){
                t = sta.top();
                ret = max(ret,t.second*(i-t.first));
                sta.pop();
            }
            sta.push(make_pair(t.first,v[i]));
        }
    }
    int total_width = v.size();
    while(not sta.empty()){
        piT t = sta.top();
        sta.pop();
        ret = max(ret,t.second*(total_width-t.first));
    }
    return ret;
}

// square area 
int solve(const vvi &matrix){
    int h = matrix.size();
    int w = matrix[0].size();
    vvi up(h,vi(w));

    {
        for(int i=0;i<h;i++){
            for(int j=0;j<w;j++){
                if(i == 0){
                    up[i][j] = matrix[i][j] == 1;
                }else {
                    if(matrix[i][j] == 0){
                        up[i][j] = 0;
                    }else{
                        up[i][j] = up[i-1][j]+1;
                    }
                }
            }
        }
    }

    int ret = 0;
    for(int y=0;y<h;y++){
        ret = max(ret,max_rectangle_size_in_histogram(up[y]));
    }
    return ret;
}

なんか長方形の和をもとめるやつ¶

#include <vector>
#include <iostream>

using namespace std;

// ret[i][j] = sum of all field in [0,i) x [0,j).
template<typename T>
vector<vector<T>> calc_sum(const vector<vector<T>>& field){
    int h = field.size();
    int w = field[0].size();
    vector<vector<T>> ret(h+1,vector<T>(w+1));
    for(int y=1;y<=h;y++){
        for(int x=1;x<=w;x++){
            ret[y][x] = ret[y][x-1] + field[y-1][x-1];
        }
    }
    for(int x=0;x<=w;x++){
        for(int y=1;y<=h;y++){
            ret[y][x] = ret[y][x] + ret[y-1][x];
        }
    }
    return ret;
}

// calc [sy,gy] x [sx,gx]
template<typename T>
T calc(const vector<vector<T>>& v,int sy,int sx,int gy,int gx){
    // for close section
    gy++;gx++;
    return v[gy][gx] - v[sy][gx] - v[gy][sx] + v[sy][sx];
}

int main(){
    vector<vector<int>> field = {{0,1,2},
                                 {3,4,5},
                                 {6,7,8}};
    auto s = calc_sum(field);
    // calc 3+5+6 + 6+7+8
    cout << calc(s,1,0,2,2) << endl;

    return 0;
}

// 座標圧縮
int main(){
    int n,m;cin >> n >> m;

    vector<int> index_of_x(n),index_of_y(n);
    for(int i=0;i<n;i++){
        cin >> index_of_x[i] >> index_of_y[i];
    }
    vector<int> value_of_x = index_of_x,
                value_of_y = index_of_y;

    sort(all(value_of_x));
    sort(all(value_of_y));
    value_of_x.erase(unique(all(value_of_x)),value_of_x.end());
    value_of_y.erase(unique(all(value_of_y)),value_of_y.end());
    for(int i=0;i<n;i++){
        index_of_x[i] = lower_bound(all(value_of_x),index_of_x[i]) - value_of_x.begin();
        index_of_y[i] = lower_bound(all(value_of_y),index_of_y[i]) - value_of_y.begin();
    }
    vector<vector<int>> p(value_of_y.size(),vector<int>(value_of_x.size()));
    for(int i=0;i<n;i++){
        p[index_of_y[i]][index_of_x[i]] += 1;
    }
    auto v = calc_sum(p);

    for(int i=0;i<m;i++){
        int x1,y1,x2,y2;
        cin >> x1 >> y1 >> x2 >> y2;

        int sx = lower_bound(all(value_of_x),x1) - value_of_x.begin();
        int gx = upper_bound(all(value_of_x),x2) - value_of_x.begin();
        int sy = lower_bound(all(value_of_y),y1) - value_of_y.begin();
        int gy = upper_bound(all(value_of_y),y2) - value_of_y.begin();

        int ans = v[gy][gx] - v[sy][gx] - v[gy][sx] + v[sy][sx];
        cout << ans << endl;
    }

    return 0;
}

累積和¶

一次元¶

// 累積和を計算する
// - ret[k] = 0 <= i < k を満たす v[i] の和
// - ret[a] - ret[b] = b <= i < a を満たす v[i] の部分和
// ret[0] = init
// ret[1] = v[0] + init
// ret[2] = v[1] + (v[0] + init)
// ruisekiwa({1, 2, 3}) = {0, 1, 3, 6}
template<typename T, typename R = T>
std::vector<R> ruisekiwa(const std::vector<T>& v, std::function<R(const T&, const R&)> op = std::plus<>(), R init = R {})
{
    std::vector<R> ret;
    ret.reserve(v.size() + 1);
    ret.push_back(init);
    for (size_t i = 1; i <= v.size(); i++) {
        ret.push_back(op(v[i-1], ret.back()));
    }
    return ret;
}
// 累積和へのアクセス
// x = [l, r) を満たす v[x] の和を返す
template<typename R>
R bubunwa(const std::vector<R>& v, int r, int l = 0)
{
    assert(r >= l);
    return v[r] - v[l];
}

二次元¶

// 2 次元の累積和
// - ret[a][b] = 0 <= y < a && 0 <= b < b を満たす v[x][y] の和
//  (bubunwa2 参照)
// ref: https://atcoder.jp/contests/abc366/submissions/56581627, 鉄則本 p60
template<typename T>
std::vector<std::vector<T>> ruisekiwa2(const std::vector<std::vector<T>>& v, std::function<T(const T&, const T&)> op = std::plus<>(), T init = T {})
{
#if DEBUG
    assert(v.size() > 0);
    size_t min_col = v[0].size(), max_col = v[0].size();
    for (const auto& c : v) {
        min_col = min(min_col, c.size());
        max_col = max(max_col, c.size());
    }
    assert(min_col == max_col && min_col > 0);
#endif
    size_t row = v.size();
    size_t col = v[0].size();
    std::vector<std::vector<T>> ret(row + 1, std::vector<T>(col + 1, init));
    // 横の累積和
    for (size_t r = 1; r <= row; r++)
        for (size_t c = 1; c <= col; c++)
            ret[r][c] = ret[r][c-1] + v[r-1][c-1];

    // 縦の累積和
    for (size_t r = 1; r <= row; r++)
        for (size_t c = 1; c <= col; c++)
            ret[r][c] += ret[r-1][c];

    return ret;
}

// 累積和へのアクセス
// y = [fy, ty) かつ y = [fx, tx) を満たす v[y][x] の和を返す
template<typename R>
R bubunwa2(const std::vector<std::vector<R>>& v, int ty, int tx, int fy, int fx)
{
    assert(tx >= fx);
    assert(ty >= fy);
    return v[ty][tx] - v[fy][tx] - v[ty][fx] + v[fy][fx];
}

ダブリング¶

// one に対し、 x で表される binary_op を n 回適用する
template<typename T>
T apply_doubling(const T& x, const uint64_t n, const std::function<T(T, T)>& binary_op, const T& one)
{
    if (n == 0) return one;
    T ret = apply_doubling(binary_op(x, x), n/2, binary_op, one);
    if (n % 2 == 1) ret = binary_op(ret, x);
    return ret;
}

// trans で表されるような遷移を n 回適用する
// https://atcoder.jp/contests/abc167/tasks/abc167_d
std::vector<int> apply_trans_vector(const uint64_t n, const std::vector<int>& trans)
{
    std::vector<int> init(trans.size());
    iota(init.begin(), init.end(), 0);
    return apply_doubling<vector<int>>(trans, n, [n = trans.size()](auto x, auto y) {
        vector<int> ret(n);
        for (size_t i = 0; i < n; i++)
            ret[i] = x[y[i]];
        return ret;
    }, init);
}

// n 回 trans で遷移した時の場所とコストを計算する
// trans[0].first = 行き先、 trans[0].second = コスト
// https://atcoder.jp/contests/abc370/tasks/abc370_f
template<typename Cost>
std::vector<pair<int, Cost>> apply_trans_vector2(const uint64_t n, std::vector<pair<int, Cost>> trans)
{
    std::vector<pair<int, Cost>> init(trans.size());
    for (size_t i = 0; i < init.size(); i++) {
        init[i].first = i;
        init[i].second = 0;
    }

    return apply_doubling<std::vector<pair<int, Cost>>>(trans, n, [n = trans.size()](const auto& f, const auto& x) {
        std::vector<pair<int, Cost>> ret(n);
        for (size_t i = 0; i < n; i++) {
            ret[i].first = f[x[i].first].first;
            // ここの式がややこしいが、 現在の値 x を f で移すイメージがわかりやすい
            // 現在値 x[i].first から生えているエッジが f[x[i].first] なのでそのコストを現在の値に追加する
            ret[i].second = x[i].second + f[x[i].first].second;
        }
        return ret;
    }, init);
}

// x を n 乗する
template<typename T, typename R>
constexpr T powi2(T x, R n, T one = 1) {
    return apply_doubling(x, n, std::multiplies<T>(), one);
}

グラフ¶

構成要素¶

隣接行列を使うか、vector<Edge>みたいのを使うかの二択。場合によってはNode みたいなのも使うかも。隣接行列を使うとメモリとか探すのとか重い。

struct Edge{
    int to,cost;
    Edge(int to,int cost) : to(to),cost(cost) {};
};

ベルマンフォード¶

ある点から他の点への最短経路を計算する。負の閉路があっても対応可能で、負の閉路の検出も可能 O(NE)

template<typename T>
struct Edge {
    int to;
    T value;
    Edge(int to, const T& v) : to(to), value(v) {}
};

template<typename T>
struct BellmanFordResult {
    bool has_negative_cycle;
    vector<T> dist;
    vector<int> prev;
    const T INF = std::numeric_limits<T>::max() / 2;
};

// 始点 s からの全てのノードへの最短距離を計算する
// result.has_negative_cycle: 負の閉路が存在するかどうか
// result.dist: それぞれのノードへの最短距離、ただし到達不可能な時は result.INF, 負の閉路の影響で無限に負になるときは -result.INF
// result.prev: 最短経路を実現するための経路木
template<typename T>
BellmanFordResult<T> bellman_ford(const vector<vector<Edge<T>>>& graph, int s) {
    BellmanFordResult<T> result;
    result.has_negative_cycle = false;
    result.dist = vector<T>(graph.size(), result.INF);
    result.dist[0] = 0;
    result.prev = vector<int>(graph.size(), -1);

    for (size_t i = 0; i < graph.size(); i++) {
        for(size_t j = 0; j < graph.size(); j++) {
            for(size_t k = 0; k < graph[j].size(); k++) {
                const auto& e = graph[j][k];
                if(result.dist[e.to] > result.dist[j] + e.value){
                    result.dist[e.to] = result.dist[j] + e.value;
                    result.prev[e.to] = j;
                    if(i == graph.size()-1){
                        result.dist[e.to] = -result.INF;
                        result.has_negative_cycle = true;
                    }
                }
            }
        }
    }
    return result;
}

ダイクストラ¶

負の経路があったらダメ

#include <iostream>
#include <vector>
#include <queue>
#include <utility>
#include <algorithm>
#include <functional>

using namespace std;

typedef pair<int,int> pii;

struct Edge{
    int from,to,cost;
    Edge(int from,int to,int cost)
        : from(from),to(to),cost(cost) {};
};

int dijkstra(const int start,const int goal,
             const vector<vector<Edge> > &graph){
    typedef pair<int,int> pii;
    int N = graph.size();
    vector<char> visited(N,false);
    // cost where
    priority_queue<pii,vector<pii>,greater<pii> > que;
    que.push(make_pair(0,start));
    while(not que.empty()){
        int cost,where;
        cost = que.top().first;
        where = que.top().second;
        que.pop();
        if(visited[where]) continue;
        if(where == goal) return cost;
        visited[where] = true;
        for(int j=0;j<(int)graph[where].size();j++){
            que.push(make_pair(graph[where][j].cost+cost,graph[where][j].to));
        }
    }
    return -1;
}


int main(){
    int n,m;
    cin >> n >> m;
    vector<vector<Edge> > V(m);

    for(int i=0;i<n;i++){
        int a,b,cost;
        cin >> a >> b >> cost;
        V[a].push_back(Edge(a,b,cost));
        V[b].push_back(Edge(b,a,cost));
    }

    int ret = -1;
    int p,q;
    cin >> p >> q;
    vector<char> visited(m,false);
    //                 cost,where
    priority_queue<pii,vector<pii>, greater<pii> > Q;
    Q.push(make_pair(0,p));
    while(!Q.empty()){
        int cost,where;
        cost = Q.top().first;
        where = Q.top().second;
        Q.pop();
        if(visited[where]) continue;
        if(where == q){
             ret = cost;
             break;
        }
        visited[where] = true;
        for(int j=0;j<(int)V[where].size();j++){
            Q.push(make_pair(V[where][j].cost+cost,V[where][j].to));
        }
    }
    // 到達不能なときは-1
    cout << ret << endl;
    return 0;
}

直径¶

木の直径 (最も遠い２点の距離) を求める O(n)

struct DiameterAnswer {
    int s, t; // 直径になる２点
    ll distance; // 直径
};
DiameterAnswer tree_diameter(const vector<vector<Edge>>& graph)
{
    ll INF = 1ll << 60;
    assert(graph.size() >= 1);
    // 適当なノードからスタートし、一番遠い点 s を求める
    // その後その点から一番遠い点 t を求める。
    // s -- t が木の直径

    // ノード番号、距離
    std::function<pair<int, ll>(int, int)> find_farthest_node = [&](int c, int from) {
        auto ret = make_pair(c, 0ll);
        for (const Edge& e : graph[c]) {
            if (e.to == from) continue;
            auto p = find_farthest_node(e.to, c);
            p.second += e.cost;
            if (p.second > ret.second) {
                ret = p;
            }
        }
        dump(ret);
        return ret;
    };
    int s = find_farthest_node(0, -1).first;
    auto p = find_farthest_node(s, -1);
    int t = p.first;
    auto dist = p.second;

    return {s, t, dist};
}

ワーシャルフロイド¶

すべてのノードに対してすべてのノードへの距離を求める。負の経路があってもOK。もし負の閉路があったらiからiはマイナスになる。 $O(|V|^{3})$

// 隣接行列の形でグラフを受け取って、 全点間距離を計算する。 O(頂点の数 ^3)
// 渡す際に graph[0][0] = 0 としなければならず、また多重辺があるなら最小の辺を行列の値としておく。
// 負の辺があっても良く、負の閉路があった場合は graph[i][i] < 0 になる。
template<typename T>
void warshall_floyd(std::vector<std::vector<T>>& graph)
{
#ifdef DEBUG
    // N x N の行列でなければならない
    for (size_t i = 0; i < graph.size(); i++) {
        assert(graph[i].size() == graph.size());
    }
    // graph[x][x] = 0 でなければならない
    for (size_t i = 0; i < graph.size(); i++) {
        assert(graph[i][i] == 0);
    }
#endif
    const size_t n = graph.size();
    for(size_t k = 0; k < n; k++){
        for(size_t i = 0; i < n; i++){
            for(size_t j = 0; j < n; j++){
                graph[i][j] = std::min(graph[i][j], graph[i][k] + graph[k][j]);
            }
        }
    }
}

最小全域木¶

クラスカル法¶

クラスカル法によって最小全域木を構築する、大体 O(MlogM) (M = グラフの辺の数) 辺をコストによってソートしておいて、もしその辺を追加したら連結成分が更新される時だけ追加するプリム法より早めかもしれない。

template<typename T>
struct MinimumSpanningTreeResult {
    bool success; // そもそも全域木が構築できたかどうか
    vector<vector<Edge<T>>> graph; // 最小全域木
    T total_cost; // 最小全域木全体のコスト
};

// クラスカル法によって最小全域木を構築する、大体 O(MlogM) (M = グラフの辺の数)
// 辺をコストによってソートしておいて、もしその辺を追加したら連結成分が更新される時だけ追加する
// プリム法より早めかもしれない。
// Ref: https://atcoder.jp/contests/typical90/tasks/typical90_aw
template<typename T>
MinimumSpanningTreeResult<T> minimum_spanning_tree(const vector<vector<Edge<T>>>& graph)
{
    // UnionFind を使う、ありものを使ってもいいがライブラリ間の依存を減らすために内部で定義することにする
    vector<int> uf(graph.size(), -1); // uf の値は負ならその集合のサイズに -1 をかけたもの、非負なら親のインデックス
    // 親のインデックスを返す
    function<int(int)> root = [&uf, &root](int x) {
        if (uf[x] < 0) return x;
        return uf[x] = root(uf[x]);
    };
    // x と y を併合
    auto unite = [&uf, &root](int x, int y) {
        x = root(x);
        y = root(y);
        if (x == y) return;
        // 集合 x に 集合 y を併合する、 もし y の方が大きいなら x, y を入れ替え
        if (uf[y] < uf[x]) swap(x, y);
        uf[x] += uf[y];
        uf[y] = x;
    };

    // Edge 構造体に辺の元の情報がないので追加しておく
    vector<pair<int, Edge<T>>> sorted_edges;
    for (int i = 0; i < static_cast<int>(graph.size()); i++) {
        for (const auto& e : graph[i]) {
            sorted_edges.push_back({i, e});
        }
    }
    // コストでソート
    std::sort(sorted_edges.begin(), sorted_edges.end(), [](const auto& e1, const auto& e2) {
        return e1.second.value < e2.second.value;
    });
    T total_cost = 0;
    vector<vector<Edge<T>>> minimum_spanning_tree(graph.size());
    for (const auto& e : sorted_edges) {
        int from = e.first;
        int to = e.second.to;
        const auto& cost = e.second.value;
        if (root(from) == root(to)) continue; // すでに同じ連結成分に存在
        unite(from, to);
        total_cost += cost;
        minimum_spanning_tree[from].push_back(e.second);
    }

    int largest_group_size = -*min_element(uf.begin(), uf.end());
    if (largest_group_size != static_cast<int>(graph.size())) {
        // 全域木が構築できていない
        return MinimumSpanningTreeResult<T>{ false, {}, 0 };
    }
    return MinimumSpanningTreeResult<T> { true, minimum_spanning_tree, total_cost };
}

プリム法¶

$O(N^2log(N))$ だと思う。最小コストを求めるコードが以下。

struct Edge{
    int to, cost;
    Edge(int to_,int cost_)
        : to(to_),cost(cost_) {}
};


vector<vector<Edge>> minimum_spanning_tree(const vector<vector<Edge>> &graph){
    typedef tuple<int,int,int> cft;  // cost,from,to
    int N = graph.size();
    vector<vector<Edge>> ret(N);
    vector<char> used(N,false);
    priority_queue<cft,vector<cft>,greater<cft>> que;
    que.push(make_tuple(0,-1,0));
    while(not que.empty()){
        int cost = get<0>(que.top());
        int from = get<1>(que.top());
        int to = get<2>(que.top());
        que.pop();
        if(used[to]) continue;
        used[to] = true;
        // ignore first.
        if(from != -1){
            ret[from].push_back(Edge(to,cost));
        }
        for(const Edge& e : graph[to]){
            que.push(make_tuple(e.cost,to,e.to));
        }
    }
    return ret;
}

閉路の検出¶

// グラフにある閉路を返す。 ({0, 1, 2, 3, 0}) のように最初と最後の値が同じものを返す
// DAG (Directed acyclic graph) かどうかの判定にも利用できる
// https://drken1215.hatenablog.com/entry/2023/05/20/200517
// TODO: 未検証: is_directed_graph が false、つまり無向グラフにこの関数が呼ばれた時、 同じ辺を二度使うようなループは許さない
template<typename T>
optional<vector<int>> find_cycle(const vector<vector<Edge<T>>>& graph, bool is_directed_graph = true)
{
    if (!is_directed_graph) {
        // 無向グラフの時、 多重辺が存在したときにはそれを使ってすぐにループを作ることができる。
        set<pair<int, int>> already_appears;
        for (int c = 0; c < static_cast<int>(graph.size()); c++) {
            for (const auto& e : graph[c]) {
                if (already_appears.contains({c, e.to})) {
                    return vector<int> { c, e.to, c };
                }
                already_appears.insert({c, e.to});
            }
        }
    }

    // 自己ループがあったらそれを返す
    for (int c = 0; c < static_cast<int>(graph.size()); c++) {
        for (const auto& e : graph[c]) {
            if (e.to == c) {
                return vector<int> { c, c };
            }
        }
    }

    vector<bool> seen(graph.size()); // dfs を呼び出したか
    vector<bool> finished(graph.size()); // dfs を最後まで完了したか
    vector<int> history;
    std::function<int(int, int)> dfs = [&](int k, int p) -> int {
        seen[k] = true;
        history.push_back(k);
        for (const auto& e : graph[k]) {
            int v = e.to;
            // 無向グラフの時、すぐに逆に戻るようなもの (同じ辺を二度使う) はスキップ
            if (!is_directed_graph && v == p) continue;
            if (finished[v]) continue; // 操作済み
            if (seen[v] && !finished[v]) {
                // この呼び出しは v の中で発生しているので v を始点終点とするループがある
                history.push_back(v);
                return v;
            }
            int pos = dfs(v, k);
            if (pos >= 0) return pos;
        }
        finished[k] = true;
        history.pop_back();
        return -1;
    };

    for (int c = 0; c < static_cast<int>(graph.size()); c++) {
        if (seen[c]) continue;
        history.clear();
        int pos = dfs(c, -1);
        if (pos >= 0) {
            vector<int> ret;
            assert(history.size() >= 2);
            // history の最後はループの最初・最後の値。ループ開始前の値をここで削除する
            while (!history.empty()) {
                ret.push_back(history.back());
                history.pop_back();
                if (ret.size() >= 2 && ret.back() == pos) break; // ループの先頭までたどり着いた
            }
            // 逆順に追加したので戻す
            reverse(ret.begin(), ret.end());
            return ret;
        }
    }
    return nullopt;
}

最大流¶

Dinic 法と FordFulkerson 法がある、 FordFulkerson は計算量がフローの総量に依存する。 Dinic 法はノードの数とエッジの数に依存する。特に二部マッチングの時、 FordFulkerson は O(n^2)、 Dinic 法は O(n^(3/2)) になる。また、最大流最小カット定理より、最大流と最小カットは一致する。

template<typename Capacity>
class MaxFlow {
public:
    virtual ~MaxFlow() = default;
    virtual int add_edge(int from, int to, Capacity cap) = 0;
    virtual Capacity flow(int from, int to) = 0;
};

// フローを流す処理は ...
template<typename Capacity>
struct Dinic : public MaxFlow<Capacity> {
    static_assert(std::is_integral_v<Capacity>, "Capacity must be integral");

    explicit Dinic(size_t number_of_node)
    : graph(number_of_node)
    {
    }

    struct Edge {
        int to; // エッジの行き先
        int rev; // graph[to][rev] が残余エッジ
        Capacity flow; // 現在流れている量
        Capacity cap; // 辺の容量
        bool is_rev_edge; // 残余グラフ用に追加したエッジかどうか、フィルター用
    };

    friend ostream& operator<<(ostream& os, const Edge& e) noexcept
    {
        return os << "{" << e.to << ", " << e.flow << "/" << e.cap << "}";
    }

    // エッジを追加し、そのインデックスを返す。残余グラフがあるため返り値が連続とは限らない
    int add_edge(int from, int to, Capacity cap) override
    {
        assert(from < static_cast<int>(graph.size()));
        assert(to < static_cast<int>(graph.size()));
        graph[from].push_back(Edge { to, static_cast<int>(graph[to].size()), 0, cap, false });
        graph[to].push_back(Edge { from, static_cast<int>(graph[from].size())-1, cap, cap, true });

        return graph[from].size() - 1;
    }

    // from から to まで最大までフローを流し、流量を返す
    Capacity flow(int from, int to) override
    {
        assert(from < static_cast<int>(graph.size()));
        assert(to < static_cast<int>(graph.size()));
        Capacity total_flow = 0;

        // Dinic 法では from から to の最短経路 DAG (閉路のない有向グラフ) を BFS で作成し、
        // DFS でその上の増加パスに流せるだけ流す、を繰り返す
        while (true) {
            auto level = levels(from);
            if (level[to] < 0) break; // そもそも to まで到達不可能なら終了する
            vector<int> iter(graph.size());
            // DAG 上で流せるだけ流す
            while (true) {
                Capacity f = good_path(from, to, std::numeric_limits<Capacity>::max(), iter, level);
                if (f == 0) break;
                total_flow += f;
            }
        }
        return total_flow;
    }

    // from から伸びているエッジを返す、デフォルトで残余グラフはフィルターアウトする
    std::vector<Edge> get_edges(int from, bool filter_rev_edge = true)
    {
        assert(from < static_cast<int>(graph.size()));

        std::vector<Edge> ret;
        for (const auto& e : graph[from]) {
            if (!filter_rev_edge || !e.is_rev_edge)
                ret.push_back(e);
        }
        return ret;
    }

    // from から伸びているエッジを返す、index は add_edge が返した値
    Edge get_edge(int from, int index)
    {
        return graph[from][index];
    }

    size_t size() { return graph.size(); }
    bool resize(size_t n)
    {
        assert(n >= graph.size());
        graph.resize(n);
    }
private:
    // from から始まる BFS により DAG を生成する。
    // 戻り値は v[i] = from から到達可能なら from からの距離、そうでないなら -1
    vector<int> levels(int from) {
        vector<int> level(graph.size(),-1);
        level[from] = 0;
        queue<int> que; que.push(from);
        while(!que.empty()){
            int v = que.front(); que.pop();
            for (auto& e : graph[v]) {
                // まだフローが流せて、かつまだ行っていないノードを追加
                if(e.cap - e.flow > 0 and level[e.to] == -1){
                    level[e.to] = level[v]+1;
                    que.push(e.to);
                }
            }
        }
        return level;
    }

    // levels で生成された DAG 上を DFS してフローを流す
    Capacity good_path(int from, int to, Capacity flow, vector<int>& iter, const vector<int>& level){
        if(from == to) return flow;
        for(int &i = iter[from]; i < static_cast<int>(graph[from].size()); i++){
            auto& e = graph[from][i];
            if(e.cap - e.flow > 0 && level[from] < level[e.to]){
                Capacity nf = good_path(e.to, to, min(flow, e.cap - e.flow), iter, level);
                if(nf > 0){
                    e.flow += nf;
                    graph[e.to][e.rev].flow -= nf;
                    return nf;
                }
            }
        }
        return 0;
    }
    std::vector<std::vector<Edge>> graph;
};

FordFulkerson は以下。あまり使わないかも。

// フローを流す処理は O(FM) (F = フローの総量、 M = エッジの数)
template<typename Capacity>
struct FordFulkerson : public MaxFlow<Capacity> {
    static_assert(std::is_integral_v<Capacity>, "Capacity must be integral");

    explicit FordFulkerson(size_t number_of_node)
    : graph(number_of_node)
    {

    }

    struct Edge {
        int to; // エッジの行き先
        int rev; // graph[to][rev] が残余エッジ
        Capacity flow; // 現在流れている量
        Capacity cap; // 辺の容量
        bool is_rev_edge; // 残余グラフ用に追加したエッジかどうか、フィルター用
    };

    friend ostream& operator<<(ostream& os, const Edge& e) noexcept
    {
        return os << "{" << e.to << ", " << e.flow << "/" << e.cap << "}";
    }

    // エッジを追加し、そのインデックスを返す。残余グラフがあるため返り値が連続とは限らない
    int add_edge(int from, int to, Capacity cap) override
    {
        assert(from < static_cast<int>(graph.size()));
        assert(to < static_cast<int>(graph.size()));
        graph[from].push_back(Edge { to, static_cast<int>(graph[to].size()), 0, cap, false });
        graph[to].push_back(Edge { from, static_cast<int>(graph[from].size())-1, cap, cap, true });

        return graph[from].size() - 1;
    }

    // from から to まで最大までフローを流し、流量を返す
    Capacity flow(int from, int to) override
    {
        assert(from < static_cast<int>(graph.size()));
        assert(to < static_cast<int>(graph.size()));
        Capacity total_flow = 0;

        // フォードフルカーソンアルゴリズムは繰り返し dfs を行い、残余グラフ含めたグラフにフローを流す
        // 最悪ケースは 1 ずつ流量が増えていく、 dfs は O(エッジの数) であり、この dfs を総量分行うことになる。
        while (true) {
            std::vector<bool> visited(graph.size(), false);
            Capacity f = dfs(from, to, std::numeric_limits<Capacity>::max(), visited);
            if (f == 0) break; // 新しくフローが流せなかった
            total_flow += f;
        }

        return total_flow;
    }

    // from から伸びているエッジを返す、デフォルトで残余グラフはフィルターアウトする
    std::vector<Edge> get_edges(int from, bool filter_rev_edge = true)
    {
        assert(from < static_cast<int>(graph.size()));

        std::vector<Edge> ret;
        for (const auto& e : graph[from]) {
            if (!filter_rev_edge || !e.is_rev_edge)
                ret.push_back(e);
        }
        return ret;
    }

    // from から伸びているエッジを返す、index は add_edge が返した値
    Edge get_edge(int from, int index)
    {
        return graph[from][index];
    }

private:
    // from から to まで dfs を行い、その途中の最小の容量を f として保存する
    int dfs(int from, int to, Capacity flow, vector<bool>& visited)
    {
        if (from == to) return flow;
        visited[from] = true;
        for (auto& v : graph[from]) {
            if (v.flow == v.cap) continue; // すでにマックスまで流している
            if (visited[v.to]) continue; // 訪問済み、ループ
            Capacity nf = dfs(v.to, to, min(flow, v.cap - v.flow), visited);
            if (nf >= 1) {
                // 残余グラフの flow を増やしつつ、 通常グラフの flow は減らす
                v.flow += nf;
                graph[v.to][v.rev].flow -= nf;
                return nf;
            }
        }
        return 0; // どの辺にも流せなかった
    }
    std::vector<std::vector<Edge>> graph;
};

例

int main(){
    // https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_6_A&lang=ja
    int v, e; cin >> v >> e;
    int source = 0;
    int sink = v-1;

    Dinic<int> maxflow(v);
    for (int i = 0; i < e; i++) {
        int u, v, c; cin >> u >> v >> c;
        maxflow.add_edge(u, v, c);
    }
    cout << maxflow.flow(source, sink) << endl;
    return 0;
}

燃やす埋める (Project Selection)¶

変数 x1, x2 ... xn があり、それぞれ 0 か 1 のどちらかを選択できる。 xi を 0 にしたときのコストを Fi, 1 にしたときのコストを Ti とする。さらに xi = 1 かつ xj = 0 の時追加のコスト Cij がかかる。 (俗に 0/1 を燃やす・埋めるにした問題になぞらえて燃やす・埋めると呼ばれているらしい) https://drken1215.hatenablog.com/entry/2023/11/24/034300

上記の問題を最小カット (グラフの頂点を S 側、 T 側にわけ、 S -> T の辺の重さの最小を求める問題) に帰着する。

// Ref:
// - https://github.com/drken1215/algorithm/blob/master/GraphNetworkFlow/two_variable_submodular_optimization.cpp
// - https://drken1215.hatenablog.com/entry/2023/11/24/034300
// 最小カットに帰着するため以下のようにグラフを構築する。
// - source を含むカットに含まれる => 選択する (true, 1)
// - destination を含むカットに含まれる => 選択しない (false, 0)
// 内部で　Dinic 法を使うため、 O(N^2*ルールの数) になる。
// TODO: 最小カットの値を求めるだけではなく、カット (どのノードがどちらのグループに含まれるのか) を計算する方法を検討する
template<typename Cost>
struct ProjectSelection {
    explicit ProjectSelection(size_t number_of_node, Cost infinity = std::numeric_limits<Cost>::max() / 1000)
    :
    source(static_cast<int>(number_of_node)),
    destination((static_cast<int>(number_of_node) + 1)),
    flow(number_of_node + 2),
    offset(0),
    infinity(infinity)
    {
    }
    // i 番目の要素を選択したときのコストを true_cost, しない時を false_cost
    // ※ 問題が「利得を得る」であった場合は利得に -1 をかけたものを入力し、最後の答えを *-1 する。
    // この関数は選択したときの「損失」であるため
    void add_single_cost(const int i, const Cost true_cost, const Cost false_cost)
    {
        // アイデアとして、
        // - source 側にノードを含める = true にするならば true_cost がカットに含まれる
        // - destination 側にノードを含める = false にするならば false_cost がカットに含まれる
        // ようにグラフを構築する
        assert(0 <= i && i < static_cast<int>(flow.size()));
        if (true_cost >= 0 && false_cost >= 0) {
            // source から i に false_cost を生やすことで、 もし i を true にした場合は false_cost の辺は含まれず、
            // 逆に true_cost の辺がカットに含まれるようになる。
            // (source は true 側であるが、 false_cost の辺を生やすので直感とは逆になるかもしれない)
            flow.add_edge(source, i, false_cost);
            flow.add_edge(i, destination, true_cost);
        } else if (false_cost >= true_cost) {
            // フローを流す際にグラフに負の辺があると困るので、
            // 最終的な答えに +true_cost し、 false にした場合は (false_cost - true_cost) が増えることにする
            flow.add_edge(source, i, false_cost - true_cost);
            offset += true_cost;
        } else {
            // 上に同じ
            flow.add_edge(i, destination, true_cost - false_cost);
            offset += false_cost;
        }
    }
    // どちらかが true, もう片方が false の時のコストを追加する (いわゆる埋める・燃やす)
    void add_true_false_cost(const int true_i, const int false_i, const Cost cost)
    {
        assert(0 <= true_i && true_i < static_cast<int>(flow.size()));
        assert(0 <= false_i && false_i < static_cast<int>(flow.size()));
        assert(cost >= 0);
        // true となるノードから false となるノードに辺を張ることで、カットにその辺が含まれるようになる
        flow.add_edge(true_i, false_i, cost);
    }
    // true_i を選択し、 false_i を選択しないことを許さないようにする
    void ban_true_false(const int true_i, const int false_i)
    {
        // 無限のコストを追加する
        return add_true_false_cost(true_i, false_i, infinity);
    }
    // 両方が true の時の利得を追加する
    // TODO: 未検証
    void add_true_true_profit(const int true_i, const int true_j, const Cost profit)
    {
        // source から i に、 i から j にそれぞれ profit の辺をはり、オフセットに -profit する
        // true, true => 辺はどちらも source 側にあるので -profit
        // true, false => source から i の辺はカウントされないが、 i から j の辺はカット間のあるので、 0
        // false, true => source から i の辺はカウントされ、 i から j の辺はカウントされない (destination 側から source 側への辺なので) 合計 0
        // false, false => source から i の辺はカウントされ、 i から j の辺は両方 destination 側なのでカウントされない、合計 0
        assert(profit >= 0);
        assert(0 <= true_i && true_i < static_cast<int>(flow.size()));
        assert(0 <= true_j && true_j < static_cast<int>(flow.size()));

        flow.add_edge(source, true_i, profit);
        flow.add_edge(true_i, true_j, profit);
        offset += (-profit);
    }
    // 両方が false の時の利得を追加する
    // TODO: 未検証
    void add_false_false_profit(const int false_i, const int false_j, const Cost profit)
    {
        // i から destination に、 j から i にそれぞれ profit の辺をはり、オフセットに -profit する
        // true, true => i から destination の辺はカウントされ、 j から i はカウントされない。合計 0
        // true, false => i から destination の辺はカウントされ、 j から i はカウントされない (destination 側から source 側の辺なので)。合計 0
        // false, true => i から destination の辺はカウントされず、 j から i はカウントされる。 合計 0
        // false, false => i から destination の辺も j から i の辺もカウントされない -profit
        assert(profit >= 0);
        assert(0 <= false_i && false_i < static_cast<int>(flow.size()));
        assert(0 <= false_j && false_j < static_cast<int>(flow.size()));
        flow.add_edge(false_i, destination, profit);
        flow.add_edge(false_j, false_i, profit);
        offset += (-profit);
    }
    // true_is に含まれるものがすべて true の時の利得を追加する
    // TODO: 未検証
    void add_all_true_profit(const vector<int>& true_is, const Cost profit)
    {
        assert(profit >= 0);
        assert(all_of(all(true_is), [&](int i){ 0 <= i && i < static_cast<int>(flow.size());}));
        // 補助頂点 y を追加し、 source から y に profit の辺を張り、
        // y から true になるノード全てに♾️の辺をはる。またオフセットに -profit を設定
        // => y が source 側に存在する時のみ全体のコストが -profit になるが、一方で true_is の一つでも false だと無限のコストが発生する
        const int y = flow.size();
        flow.resize(flow.size() + 1);
        flow.add_edge(source, y, profit);
        for (int i : true_is) {
            flow.add_edge(y, i, infinity);
        }
        offset += (-profit);
    }
    // false_is に含まれるものがすべて false の時の利得を追加する
    // TODO: 未検証
    void add_all_false_profit(const vector<int>& false_is, const Cost profit)
    {
        assert(profit >= 0);
        assert(all_of(all(false_is), [&](int i){ 0 <= i && i < static_cast<int>(flow.size());}));
        // 補助頂点 y を追加し、 y から destination に profit の辺を張り、
        // false になるノード全てから y に♾️の辺をはる。またオフセットに -profit を設定
        // => y が destination 側に存在する時のみ全体のコストが -profit になるが、一方で false_is の一つでも true だと無限のコストが発生する
        const int y = flow.size();
        flow.resize(flow.size() + 1);
        flow.add_edge(y, destination, profit);
        for (int i : false_is) {
            flow.add_edge(i, y, infinity);
        }
        offset += (-profit);
    }
    // より一般的な関係コストの追加
    // TODO: 未検証
    void add_cost_general(int i, int j, Cost true_true_cost, Cost true_false_cost, Cost false_true_cost, Cost false_false_cost)
    {
        assert(0 <= i && i < static_cast<int>(flow.size()));
        assert(0 <= j && j < static_cast<int>(flow.size()));
        assert(false_true_cost + true_false_cost >= true_true_cost + false_false_cost);
        offset += true_true_cost;
        add_single_cost(i, false_false_cost - false_true_cost, 0);
        add_single_cost(j, false_true_cost - true_true_cost, 0);
        add_true_false_cost(i, j, false_false_cost + true_false_cost - true_true_cost - false_false_cost);
    }
    Cost solve()
    {
        return flow.flow(source, destination) + offset;
    }
private:
    size_t number_of_node;
    int source;
    int destination;
    Dinic<Cost> flow;
    Cost offset;
    Cost infinity;
};

以下例

// https://atcoder.jp/contests/typical90/tasks/typical90_an
int main([[maybe_unused]] int argc, [[maybe_unused]] char** argv)
{
    // 小数点以下10桁表示にする場合は以下をコメントアウト
    // cout << fixed << setprecision(10);
    int n; ll w; cin >> n >> w;
    ProjectSelection<ll> project_selection(n);
    auto a = input_vector<ll>(n);

    // 選択する = 入る
    // 選択しない = 入らない
    vector<vector<int>> keys(n);
    for (int i = 0; i < n; i++) {
        int k; cin >> k;
        for (int j = 0; j < k; j++) {
            int c; cin >> c;
            c--;
            keys[i].push_back(c);
        }
    }

    for (int i = 0; i < n; i++) {
        project_selection.add_single_cost(i, w - a[i], 0);
        // i を選択しない場合、 その家にある鍵を取れないので入れなくなる
        for (int k : keys[i]) {
            project_selection.ban_true_false(k, i);
        }
    }

    auto ret = -project_selection.solve();
    cout << ret << endl;
    return 0;
}

最小費用流¶

Primal-Dual法による。

typedef pair<int,int> pii;

struct Edge{
    int to,cap,cost,rev;
    Edge(int to,int cap,int cost,int rev)
        : to(to),cap(cap),cost(cost),rev(rev) {};
};

void add_edge(vector<vector<Edge> > &E,int from,int to,int cap,int cost){
    E[from].push_back(Edge(to,cap,cost,E[to].size()));
    E[to].push_back(Edge(from,0,-cost,E[from].size()-1));
}

// s -> t (flow f)
//  if cant, return -1.
int min_cost_flow(vector<vector<Edge> > E,int s,int t,int f){
    const int INF = 1 << 30;
    int ret = 0;
    // potential
    vector<int> h(E.size());
    vector<int> prevv(E.size());
    vector<int> preve(E.size());

    while(f > 0){
        vector<int> dist(E.size(),INF);
        dist[s] = 0;
        priority_queue<pii,vector<pii>,greater<pii> > que;
        que.push(make_pair(0,s));
        while(!que.empty()){
            pii p = que.top();
            que.pop();
            int pf = p.first,ps = p.second;
            if(dist[ps] < pf) continue;
            for(int i=0;i<E[ps].size();i++){
                Edge &e = E[ps][i];
                if(e.cap > 0 and dist[e.to] > dist[ps] + e.cost + h[ps] - h[e.to]){
                    dist[e.to] = dist[ps] + e.cost + h[ps] - h[e.to];
                    prevv[e.to] = ps;
                    preve[e.to] = i;
                    que.push(make_pair(dist[e.to],e.to));
                }
            }
        }
        if(dist[t] == INF){
            return -1;
        }
        for(int v=0;v<E.size();v++){
            h[v] += dist[v];
        }
        int d = f;
        for(int v=t;v!=s;v=prevv[v]){
            d = min(d,E[prevv[v]][preve[v]].cap);
        }
        f -= d;
        ret += d * h[t];
        for(int v=t;v!=s;v=prevv[v]){
            Edge &e = E[prevv[v]][preve[v]];
            e.cap -= d;
            E[v][e.rev].cap += d;
        }
    }
    return ret;
}

int main(){
    while(true){
        int N,M;
        cin >> N >> M;
        if(N == 0 and M == 0) break;
        vector<pair<int,int> > men;
        vector<pair<int,int> > houses;
        for(int h=0;h<N;h++){
            for(int w=0;w<M;w++){
                char x;
                cin >> x;
                if(x == 'H') houses.push_back(make_pair(h,w));
                else if(x == 'm') men.push_back(make_pair(h,w));
            }
        }

        vector<vector<Edge> > E(men.size()+houses.size()+2);
        int s = men.size()+houses.size();
        int t = s+1;

        for(int i=0;i<men.size();i++){
            add_edge(E,s,i,1,0);
            for(int j=0;j<houses.size();j++){
                int dist = abs(men[i].first - houses[j].first)
                    + abs(men[i].second - houses[j].second);
                add_edge(E,i,men.size()+j,1,dist);
            }
        }

        for(int i=0;i<houses.size();i++){
            add_edge(E,men.size()+i,t,1,0);
        }

        cout << min_cost_flow(E,s,t,men.size()) << endl;
    }
    return 0;
}

無向グラフにおける最小カット¶

Nagamochi-Ibaraki/Storer-Wagnerの方法によってO(V^3)で計算できる。

struct Edge{
    int to,cap,rev;
    Edge(int _to,int _cap,int _rev) : to(_to),cap(_cap),rev(_rev) {};
};

void add_edge(vector<vector<Edge> >& E,int from,int to,int cap){
    E[from].push_back(Edge(to,cap,E[to].size()));
    E[to].push_back(Edge(from,0,E[from].size()-1));
}

// Nagamochi-Ibaraki/Stoer-Wagner
//  http://www.prefield.com/algorithm/graph/minimum_cut.html
int minimum_cut_of_undirected_graph(const vector<vector<Edge> > &graph){
    const int INF = 1 << 30;
    int n = graph.size();
    vector<vector<int> > adj(n,vector<int>(n));
    for(int i=0;i<n;i++){
        for(size_t j=0;j<graph[i].size();j++){
            Edge e = graph[i][j];
            adj[i][e.to] += e.cap;
        }
    }
    vector<int> h(n);
    for(int i=0;i<n;i++){
        h[i] = i;
    }
    int cut = INF;
    for(int m = n;m > 1;m--){
        vector<int> ws(m,0);
        int u,v,w;
        for(int k=0;k<m;k++){
            u = v;
            v = max_element(ws.begin(),ws.end()) - ws.begin();
            w = ws[v];
            ws[v] = -1;
            for(int i=0;i<m;i++){
                if(ws[i] >= 0) {
                    ws[i] += adj[h[v]][h[i]];
                }
            }
        }
        for(int i=0;i<m;i++){
            adj[h[i]][h[u]] += adj[h[i]][h[v]];
            adj[h[u]][h[i]] += adj[h[v]][h[i]];
        }
        h.erase(h.begin()+v);
        cut = min(cut,w);
    }
    return cut;
}

int main(){
    int N,M;
    while(scanf("%d %d\n",&N,&M) != EOF){
        vector<vector<Edge> > graph(N);
        for(int i=0;i<M;i++){
            int a,b,c;
            scanf("%d %d %d\n",&a,&b,&c);
            add_edge(graph,a,b,c);
            add_edge(graph,b,a,c);
        }
        printf("%d\n",minimum_cut_of_undirected_graph(graph));
    }
    return 0;

強連結成分分解¶

template<typename T>
vector<vector<int>> strongly_connected_components(const vector<vector<Edge<T>>>& graph){
    const int n = graph.size();
    vector<int> back_number;

    {
        vector<bool> used(n);
        function<void(int)> check_back_number = [&](int v) {
            used[v] = true;
            for(const auto& e : graph[v]){
                if(!used[e.to]){
                    check_back_number(e.to);
                }
            }
            back_number.push_back(v);
        };
        for(int i = 0; i < n; i++){
            if(!used[i]){
                check_back_number(i);
            }
        }
    }

    vector<vector<int>> scc;

    {
        vector<vector<Edge<T>>> reversed_graph(n);
        for(int from = 0; from < n; from++){
            for(const auto& e : graph[from]){
                reversed_graph[e.to].push_back(Edge(from, e.value));
            }
        }
        vector<char> used(n);
        function<void(int, vector<int>&)> collect_nodes = [&](int v, vector<int>& s) {
            used[v] = true;
            s.push_back(v);
            for(const auto& e : reversed_graph[v]){
                if (!used[e.to]) {
                    collect_nodes(e.to,s);
                }
            }
        };

        reverse(back_number.begin(), back_number.end());
        for(int k : back_number){
            if(!used[k]){
                scc.emplace_back(0);
                collect_nodes(k, scc.back());
            }
        }
    }
    return scc;
}


// SCC した後のグラフを計算。
template<typename T>
vector<vector<Edge<T>>> construct_scced_graph(const vector<vector<Edge<T>>>& graph)
{
    auto scc = strongly_connected_components(graph);
    dump(scc);
    vector<vector<Edge<T>>> scc_graph(scc.size());

    // 元のノード番号から新しいノード番号へのマップを作成
    map<int, int> to_scc_node_id;
    for (int i = 0; i < scc.size(); i++) {
        for (int c : scc[i]) {
            to_scc_node_id[c] = i;
        }
    }

    // SCC 後のグラフにエッジをはる
    for (int i = 0; i < scc.size(); i++) {
        for (int c : scc[i]) {
            for (const auto& e : graph[c]) {
                // 自己ループになる時はスキップ
                if (to_scc_node_id[e.to] == i) continue;
                // TODO: 多重辺をとりのぞくか、あるいはどのように取りのぞくかカスタマイズできるようにする
                scc_graph[i].emplace_back(to_scc_node_id[e.to], e.value);
            }
        }
    }
    return scc_graph;
}

2 SAT¶

namespace SAT{
    // variable is integer .
    struct Literal{
        bool is_not;
        int var;
        Literal(bool _is_not,int _var)
            : is_not(_is_not),var(_var){}
        Literal not_() const{
            return Literal(not is_not,var);
        }
    };
    Literal make_yes(int var){
        return Literal(false,var);
    }
    Literal make_no(int var){
        return Literal(true,var);
    }
    using Clause = std::vector<Literal>;
    bool solve_2SAT(const int number_of_variable,
                    const vector<Clause>& problem,
                    vector<char>& output){
        // variable is [0..N).
        // "not n" will be assigned to (n+N).
        auto index = [number_of_variable](const bool& is_not,const int& var){
            return number_of_variable*is_not+var;
        };
        auto to_int = [number_of_variable,&index](const Literal& lit){
            return index(lit.is_not,lit.var);
        };

        vector<vector<Edge>> graph(number_of_variable*2);
        for(const Clause& clause : problem){
            // clause = {Literal,Literal}
            // a \/ b -> ~a -> b /\ ~b -> a
            for(int i=0;i<2;i++){
                Literal c = clause[i].not_(),
                        d = clause[not i];
                graph[to_int(c)].push_back(Edge(to_int(d),1));
            }
        }
        vector<vector<int>> scc = strongly_connected_components(graph);
        // where[i] = which component holds i
        vector<int> where(number_of_variable*2);
        for(int i=0;i<scc.size();i++){
            for(int node : scc[i]){
                where[node] = i;
            }
        }
        for(int var=0;var<number_of_variable;var++){
            if(where[index(false,var)] == where[index(true,var)]){
                return false;
            }
        }
        for(int var=0;var<number_of_variable;var++){
            output[var] = where[index(false,var)] < where[index(true,var)];
        }
        return true;
    }
}

// Codeforces 228E
signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    int n,m;
    cin >> n >> m;
    vector<SAT::Clause> problem;
    for(int i=0;i<m;i++){
        int a,b,c;
        cin >> a >> b >> c;
        a--;b--;
        if(c){
            // already
            // ~(A/\~B) /\ ~(~A/\B)
            //  = (~A \/ B) /\ (A \/ ~B)
            problem.push_back({SAT::make_no(a),SAT::make_yes(b)});
            problem.push_back({SAT::make_yes(a),SAT::make_no(b)});
        }else{
            // ~(A/\B) /\ ~(~A/\~B)
            //  = (~A \/ ~B) /\ (A \/ B)
            problem.push_back({SAT::make_no(a),SAT::make_no(b)});
            problem.push_back({SAT::make_yes(a),SAT::make_yes(b)});
        }
    }
    vector<char> output(n);
    bool is_solvable = SAT::solve_2SAT(n,problem,output);
    if(is_solvable){
        cout << count(all(output),true) << endl;
        for(int i=0;i<n;i++){
            if(output[i]){
                cout << i+1 << " ";
            }
        }
        cout << endl;
    }else{
        cout << "Impossible" << endl;
    }

    return 0;
}

橋の列挙¶

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

#define all(c) c.begin(),c.end()
#define repeat(i,n) for(int i=0;i<static_cast<int>(n);i++)

struct Edge{
    int to,cost;
    Edge(int to_,int cost_)
        :to(to_),cost(cost_){
    }
};

int find_bridge_sub(int cur,int from,int number,
                    const vector<vector<Edge>>& graph,
                    vector<int> &pre_order,
                    vector<int> &low,
                    vector<vector<Edge> > &bridges){
    pre_order[cur] = low[cur] = number;
    number++;
    for(const Edge& e : graph[cur]){
        // not visited
        if(pre_order[e.to] == -1){
            number = find_bridge_sub(e.to,cur,number,graph,pre_order,low,bridges);
            if(pre_order[e.to] == low[e.to]){
                bridges[cur].push_back(e);
            }
        }
        if(from != e.to){ //ignore parent.
            low[cur] = min(low[cur],low[e.to]);
        }
    }
    return number;
}

// assume it is connected graph.
vector<vector<Edge>> find_bridge(const vector<vector<Edge>>& graph){
    int N = graph.size();
    vector<int> pre_order(N,-1),lowest_number_of_connected(N,N+1);
    vector<vector<Edge>> ret(N);
    find_bridge_sub(0,-1,0,graph,pre_order,lowest_number_of_connected,ret);

    // optional: make sure from < e.to
    //  You can use following if graph is bidirectional.
    // vector<vector<Edge>> swapped(N);
    // repeat(f,graph.size()){
    //     for(const Edge& e: ret[f]){
    //         if(f < e.to){
    //             swapped[f].push_back(e);
    //         }else{
    //             swapped[e.to].push_back(Edge(f,e.cost));
    //         }
    //     }
    // }
    // ret = swapped;
    return ret;
}


int main(){
    int V_SIZE,E_SIZE;
    cin >> V_SIZE >> E_SIZE;
    vector<vector<Edge>> graph(V_SIZE);
    repeat(i,E_SIZE){
        int s,t;
        cin >> s >> t;
        graph[s].push_back(Edge(t,1));
        graph[t].push_back(Edge(s,1));
    }

    vector<vector<Edge>> bridges = find_bridge(graph);
    repeat(f,graph.size()){
        vector<Edge> &es = bridges[f];
        sort(all(es),[](Edge lhs,Edge rhs){return lhs.to < rhs.to;});
        for(const Edge& e : es){
            cout << f << " " << e.to << endl;
        }
    }
    return 0;
}

Lowest Common Ancestor¶

木において、根から最も遠い、u,vの共通の祖先をLCAと呼ぶ。

// 木上で共通の祖先を高速に求める
struct LowestCommonAncestor
{
    // ノードごとの子供を vector で渡す、循環が存在しないこと
    LowestCommonAncestor(const std::vector<std::vector<int>> &children, const int root) :
        children(children), root(root) {
#ifdef DEBUG
        // 循環が存在しないこと、 ルートを子供に持つ要素がないことをチェックする
        for (int p = 0; p < static_cast<int>(children.size()); p++) {
            for (int child : children[p]) {
                assert(0 <= child && child < static_cast<int>(children.size()));
                assert(child != root);
                assert(find(children[child].begin(), children[child].end(), p) == children[child].end());
            }
        }
#endif
        assert(0 <= root && root < static_cast<int>(children.size()));
        calc_parent_and_depth();
        calc_parent_pow2();
    }

    template<typename T>
    LowestCommonAncestor(const std::vector<std::vector<Edge<T>>> &tree, const int root)
    : LowestCommonAncestor(generate_children_from_edge_list(tree, root), root)
    {
    }


    // u, v の共通の祖先を求める
    int solve(int u, int v) {
        // make sure depth(u) > depth(v).
        if (depth[u] < depth[v]) swap(u,v);
        for (size_t k = 0; k < parent_pow2.size(); k++){
            if(((depth[u] - depth[v]) >> k) & 1){
                u = parent_pow2[k][u];
            }
        }

        if (u == v) return u;
        for (int k = static_cast<int>(parent_pow2.size()) - 1; k >=0; k--) {
            if (parent_pow2[k][u] != parent_pow2[k][v]) {
                u = parent_pow2[k][u];
                v = parent_pow2[k][v];
            }
        }
        return parent_pow2[0][u];
    }

    // u, v の距離を求める
    // 3 点の合計の距離は (distance(a[0], a[1] + distance(a[1], a[2]) + distance(a[2], a[0])) / 2
    // https://atcoder.jp/contests/typical90/tasks/typical90_ai
    int distance(int u, int v)
    {
        int lca = solve(u, v);
        // root から u の距離　+ root から v の距離 - 2 * 共通祖先への距離
        return depth[u] + depth[v] - 2 * depth[lca];
    }

private:
    void calc_parent_and_depth() {
        parent = vector<int>(children.size(), -1);
        depth = vector<int>(children.size(), -1);
        sub_calc_parent_and_depth(root, -1, 0);
    }
    void sub_calc_parent_and_depth(int cur, int par, int dep){
        parent[cur] = par;
        depth[cur] = dep;
        for (int child : children[cur]) {
            if (child != par) {
                sub_calc_parent_and_depth(child,cur,dep+1);
            }
        }
    }
    void calc_parent_pow2() {
        // parent_pow2[k][i] = 2^k parent of node i.
        parent_pow2 = vector<vector<int>>(ceil(log(children.size())/log(2)),
                                        vector<int>(children.size(),-1));
        parent_pow2[0] = parent;
        for (size_t k = 0; k + 1 < parent_pow2.size(); k++) {
            for (size_t v = 0; v < children.size(); v++) {
                if (parent_pow2[k][v] >= 0) {
                    parent_pow2[k+1][v] = parent_pow2[k][parent_pow2[k][v]];
                }
            }
        }
    }

    template<typename T>
    static std::vector<std::vector<int>> generate_children_from_edge_list(const std::vector<std::vector<Edge<T>>> &tree, const int root)
    {
        vector<vector<int>> children(tree.size());
        function<void(int, int)> dfs = [&](int c, int p) {
            for (int child : tree[c]) {
                if (child == p) continue;
                children[c].push_back(child);
                dfs(child, c);
            }
        };
        dfs(root, -1);
        return children;
    }

    std::vector<std::vector<int>> children;
    int root;
    // if root,parent is -1.
    vector<int> parent;
    vector<int> depth;
    vector<vector<int>> parent_pow2;
};

Auxiliary Tree¶

// 木のいくつかの頂点を選んでそれらを含む最小の部分木を作成する
template<typename T>
struct AuxiliaryTree {
    explicit AuxiliaryTree(const vector<vector<Edge<T>>>& tree, int root = 0)
        : lca(tree, root), node_index_to_dfs_order(tree.size())
    {
        // グラフの頂点を全て巡回するような DFS の行き掛け順をメモする
        vector<int> dfs_order;
        function<void(int, int)> dfs = [&](int c, int p) {
            dfs_order.push_back(c);
            for (auto [child, _] : tree[c]) {
                if (child == p) continue;
                dfs(child, c);
            }
        };
        dfs(root, -1);
        assert(dfs_order.size() == tree.size());

        for (int i = 0; i < dfs_order.size(); i++) {
            node_index_to_dfs_order[dfs_order[i]] = i;
        }
    }

    // DFS の行き掛け順で与えられたインデックスをソートする
    vector<int> sort_by_dfs_order(std::vector<int> node_index)
    {
        std::sort(node_index.begin(), node_index.end(), [&](int x, int y) {
            return node_index_to_dfs_order[x] < node_index_to_dfs_order[y];
        });
        return node_index;
    }

    // 与えられたノードをつなげるのに必要なエッジの数を計算
    // https://atcoder.jp/contests/typical90/tasks/typical90_ai
    int number_of_edges_to_connect_nodes(const std::vector<int>& node_index)
    {
        auto sorted = sort_by_dfs_order(node_index);
        int ret = 0;
        for (size_t i = 0; i < sorted.size(); i++) {
            ret += lca.distance(sorted[i], sorted[(i + 1) % sorted.size()]);
        }
        ret /= 2; // 辺ごとに ２ 回カウントするので /2
        return ret;
    }

    // 与えられたノードの位置関係を接続するために必要な中間ノードを返す (与えられたノードも含む)
    // Auxiliary Tree に含まれるようなノードを返す
    // TODO: 未検証、エッジも貼るようにする
    vector<int> compress_tree(const std::vector<int>& node_index)
    {
        auto sorted = sort_by_dfs_order(node_index);
        size_t orig_n = sorted.size();
        for (size_t i = 0; i < orig_n - 1; i++) {
            sorted.push_back(lca.solve(sorted[i], sorted[i+1]));
        }
        return sorted;
    }

private:
    LowestCommonAncestor lca;
    vector<int> node_index_to_dfs_order;
};

トポロジカルソート¶

// 有向であってループがないグラフ (DAG) を入力とする
// もし DAG ではないものが入力されたら nullopt を返す
template<typename T>
std::optional<vector<int>> topological_sort(const std::vector<std::vector<Edge<T>>>& graph)
{
    if (find_cycle(graph)) {
        return std::nullopt;
    }
    // 実装としては入る辺が 0 になったものを順番にリストに突っ込んでいく
    vector<int> ret;
    ret.reserve(graph.size());
    vector<int> number_of_in_edges(graph.size());
    for (const auto& edges: graph) {
        for (const auto& e : edges) {
            number_of_in_edges[e.to] += 1;
        }
    }

    queue<int> que;
    for (size_t i = 0; i < number_of_in_edges.size(); i++) {
        if (number_of_in_edges[i] == 0) {
            que.push(static_cast<int>(i));
        }
    }

    while (!que.empty()) {
        int cur = que.front();
        que.pop();
        ret.push_back(cur);
        for (const auto& e : graph[cur]) {
            number_of_in_edges[e.to] -= 1;
            assert(number_of_in_edges[e.to] >= 0);
            if (number_of_in_edges[e.to] == 0) {
                que.push(e.to);
            }
        }
    }
    assert(ret.size() == graph.size());
    return ret;
}

最小道被覆問題¶

// グラフをいくつかの道に分割する、最小道被覆問題を解く
// https://kyopro.hateblo.jp/entry/2018/06/04/000659
// https://atcoder.jp/contests/abc374/tasks/abc374_g
template<typename T>
int minimum_path_cover_on_dag(const vector<vector<Edge<T>>>& graph)
{
#if DEBUG
    assert(find_cycle(graph) == nullopt);
#endif
    // 二部マッチング問題に帰着する
    Dinic<int> maxflow(2 * graph.size() + 2);
    const int source = 2 * graph.size();
    const int sink = 2 * graph.size() + 1;
    // source => 0 ~ n-1 に辺をはる
    for (int i = 0; i < graph.size(); i++)
        maxflow.add_edge(source, i, 1);
    // n ~ 2 * n -1 => sink に辺をはる
    for (int i = 0; i < graph.size(); i++)
        maxflow.add_edge(graph.size() + i, sink, 1);

    // 元のグラフで u => v の辺があるならば、 u => n + v に辺をはる
    for (int i = 0; i < graph.size(); i++) {
        for (const auto& e : graph[i]) {
            maxflow.add_edge(i, graph.size() + e.to, 1);
        }
    }
    // この時フローが流れている辺が元のグラフ上の道になる
    int ret = graph.size() - maxflow.flow(source, sink);
    return ret;
}

整数論¶

最大公約数,最小公倍数¶

ユークリッドの互除法を使う。intをlong longに置換してもいい。 O(logn)

最大公約数¶

int gcd(int a,int b){
    return b==0 ? a : gcd(b,a%b);
}

最小公倍数¶

int lcm(int a,int b){
    return a*b / gcd(a,b);
}

mod¶

long long に入らないような答えのときに mod が登場する。

modの計算式について¶

\[\begin{split}\begin{aligned} a \equiv c & \pmod m \\ b \equiv d & \pmod m\end{aligned}\end{split}\]

の条件下では以下の式が成り立つ。

\[\begin{split}\begin{aligned} a+b \equiv c+d & \pmod m \\ a-b \equiv c-d & \pmod m \\ a \times b \equiv c \times d & \pmod m\end{aligned}\end{split}\]

さらに、mが素数の場合、以下の関係が成り立つ。

\[\begin{split}\begin{aligned} a ^ m \equiv a \pmod m \\ a ^ {m-1} \equiv 1 \pmod m \\ a ^ {m-2} \equiv \frac{1}{a} \pmod m\end{aligned}\end{split}\]

つまり、 $a$ で割ることと、$a^{m-2}$ を掛けることは同じである。これは、 $C(10000,5000) \pmod p$ といった式を計算する際、次の冪乗の演算と組みあわせて用いる。

Mod を扱うクラス¶

template<typename T, typename R>
constexpr T powi(T x, R n, T one = 1) {
    if (n == 0) return one;
    T ret = powi(x * x, n / 2, one);
    if (n % 2 == 1) ret *= x;
    return ret;
}

template <uint_fast64_t MOD>
struct mod_int {
    mod_int() noexcept : m_value(0) { }
    constexpr mod_int(uint_fast64_t x) noexcept : m_value((x + MOD) % MOD)  { } // NOLINT
    constexpr mod_int(long long x) noexcept : m_value((x + MOD) % MOD) { } // NOLINT
    constexpr mod_int(int x) noexcept : m_value((x + MOD) % MOD) { }  // NOLINT
    friend constexpr mod_int operator+(mod_int lhs, const mod_int& rhs) noexcept { lhs += rhs; return lhs; }
    friend constexpr mod_int operator-(mod_int lhs, const mod_int& rhs) noexcept { lhs -= rhs; return lhs; }
    friend constexpr mod_int operator*(mod_int lhs, const mod_int& rhs) noexcept { lhs *= rhs; return lhs; }
    friend constexpr mod_int operator/(mod_int lhs, const mod_int& rhs) noexcept { lhs /= rhs; return lhs; }
    friend constexpr bool operator==(const mod_int& lhs, const mod_int& rhs) noexcept { return lhs.m_value == rhs.m_value; }
    friend constexpr bool operator!=(const mod_int& lhs, const mod_int& rhs) noexcept { return lhs.m_value != rhs.m_value; }
    constexpr mod_int& operator+=(const mod_int& rhs) noexcept { m_value += rhs.m_value; m_value %= MOD; return *this; }
    constexpr mod_int& operator*=(const mod_int& rhs) noexcept { m_value *= rhs.m_value; m_value %= MOD; return *this; }
    constexpr mod_int& operator-=(const mod_int& rhs) noexcept { m_value += (MOD - rhs.m_value); m_value %= MOD; return *this; }
    constexpr mod_int& operator/=(mod_int rhs) noexcept {
        uint_fast64_t exp = MOD - 2;
        while (exp > 0) {
            if (exp % 2 == 1)
                *this *= rhs;
            rhs *= rhs;
            exp /= 2;
        }
        return *this;
    }
    [[nodiscard]] mod_int inv() const noexcept { return mod_int(1) /= *this; }
    [[nodiscard]] uint_fast64_t value() const { return m_value; }
private:
    uint_fast64_t m_value;
};
template <uint_fast64_t MOD> std::ostream& operator<<(std::ostream& os, const mod_int<MOD>& v) { return os << v.value(); }
template <uint_fast64_t MOD> std::istream& operator>>(std::istream& is, mod_int<MOD>& v) { return is >> v.value(); }

constexpr ll MOD = 1000000007;
using md = mod_int<MOD>;

素数¶

エラトステネスの篩¶

#include <iostream>
#include <vector>
using namespace std;

//上限より余裕を取ること。
vector<bool> sieve(const int M){
    vector<bool> isPrime(M);
    for(int i=2;i<M;i++) isPrime[i] = true;
    for(int i=2;i*i<M;i++){
        if(not isPrime[i]) continue;
        for(int j=i*i;j<M;j+=i){
            isPrime[j] = false;
        }
    }
    return isPrime;
}

// リストがほしいときはこっち
vector<int> sieve(const int M){
    vector<bool> isPrime(M);
    vector<int> primes;
    for(int i=2;i<M;i++) isPrime[i] = true;
    for(int i=2;i*i<M;i++){
        if(not isPrime[i]) continue;
        for(int j=i*i;j<M;j+=i){
            isPrime[j] = false;
        }
    }
    for(int i=2;i<M;i++){
        if(isPrime[i]) primes.push_back(i);
    }
    return primes;
}

from math import *

def sieve(N):
    primes = set()
    for i in range(2,N):
        primes.add(i)

    for i in range(2,ceil(sqrt(N))):
        if i in primes:
            for j in range(i*i,N,i):
                primes.discard(j)

    return primes

素数のリストが欲しかったら、適当に突っ込むこと。実際には $O(n \log \log n)$ だけれど、大体 $O(n)$ だと思っていい。

素因数分解¶

// primes is primes under sqrt(N).
vector<int> prime_decomposition(int N,const vector<int>& primes){
    vector<int> ret;
    for(int p : primes){
        while(N % p == 0){
            N /= p;
            ret.push_back(p);
        }
    }
    if(N != 1) ret.push_back(N);
    return ret;
}

コンビネーション¶

くみあわせ。

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;
typedef long long ll;

//いい感じのやつ (n=61まで大丈夫)
//  modにしても良い。O(N)
ll nCr(int n,int r){
    if(n < r) return 0;
    r = min(r,n-r);
    ll ret = 1;
    for(int i=0;i<r;i++){
        ret *= n-i;
        ret /= i+1;
    }
    return ret;
}


// http://www.ravco.jp/cat/view.php?cat_id=6179
// 区別するn 個のものから，繰り返し用いることを許して，r 個取り出して作った組
// <-> 区別しないr 個のボールを，区別するn 個の箱に配る場合の総数
ll nHr(int n, int r){
    return nCr(r+n-1,r);
}

//パスカルの三角形 (n=66までオーバーフローしない)
//たくさん必要になるときはこっちのほうがいい。
ll combi_tri(int n,int r){
    int N = n+1;
    vector<vector<ll> > memo(N,vector<ll>(N,0));
    memo[0][0] = 1;
    for(int i=1;i<N;i++){
        memo[i][0] = memo[i-1][0];
        for(int j=1;j<N;j++){
            memo[i][j] = memo[i-1][j-1] + memo[i-1][j];
        }
    }
    return memo[n][r];
}

// testing...
template <typename BidirectionalIterator>
bool next_combination(BidirectionalIterator first1,
                      BidirectionalIterator last1,
                      BidirectionalIterator first2,
                      BidirectionalIterator last2){
    if((first1 == last1) or (first2 == last2)){
        return false;
    }
    BidirectionalIterator m1 = last1;
    BidirectionalIterator m2 = last2;
    --m2;
    while(--m1 != first1 && not (*m1 < *m2)){
    }
    bool result = (m1 == first1) and not (*first1 < *m2);
    if(not result){
        while(first2 != m2 and not(*m1 < *first2)){
            ++first2;
        }
        first1 = m1;
        iter_swap(first1,first2);
        ++first1;++first2;
    }
    if((first1 != last1) and (first2 != last2)){
        m1 = last1; m2 = first2;
        while((m1 != first1) and (m2 != last2)){
            iter_swap(--m1,m2);
            ++m2;
        }
        reverse(first1,m1);
        reverse(first1,last1);
        reverse(m2,last2);
        reverse(first2,last2);
    }
    return not result;
}

template<typename BidirectionalIterator>
bool next_combination(BidirectionalIterator first,
                      BidirectionalIterator middle,
                      BidirectionalIterator last){
    return next_combination(first,middle,middle,last);
}

template<class BidirectionalIterator>
inline bool prev_combination(BidirectionalIterator first,
                             BidirectionalIterator middle,
                             BidirectionalIterator last){
  return next_combination(middle,last,first,middle);
}

template<typename T>
ostream& operator<<(ostream& os,const vector<T>& vec){
    os << "[";
    for(const auto& v : vec){
        os << v << ",";
    }
    os << "]";
    return os;
}

typedef unsigned long long ull;

void print_bit(ull S, int n=64){
    for(int i=n-1; i>=0; i--){
        if(S>>i & 1) std::cout << 1;
        else std::cout << 0;
    }
    std::cout << std::endl;
}

void subset_combination(int n, int k){
    ull S = (1ULL << k) - 1ULL;
    ull E = ~((1ULL << n) - 1ULL);
    while(!(S & E)){
        print_bit(S, n);
        ull smallest = S & -S;
        ull ripple = S + smallest;
        ull nsmallest = ripple & -ripple;
        S = ripple | (((nsmallest / smallest) >> 1) - 1);
    }
}


int main(){
    // poyoから2つえらぶ
    vector<int> poyo = {1,2,3,4,5};
    int N = poyo.size(),R = 3;
    do{
        vector<int> c(poyo.begin(),poyo.begin()+R);
        cerr << c  << endl;
    }while(next_combination(poyo.begin(),poyo.begin()+R,poyo.end()));

    return 0;
}

// mod 上の　nCr が大量に必要になるときは n!, 1/r!, 1/(n-r)! を先に求めておく
// https://atcoder.jp/contests/abc042/submissions/57187098
int main() {
    vector<md> exp(n+1);
    exp[0] = 1;
    for (int i = 1; i < static_cast<int>(exp.size()); i++) {
        exp[i] = exp[i-1] * i;
    }
    vector<md> dexp(exp.size());
    for (int i = 0; i < static_cast<int>(exp.size()); i++) {
        dexp[i] = 1 / exp[i];
    }
    // n 個から r 個を選ぶ組み合わせ数、 r の方が大きい時は 0
    auto nCr = [&exp, &dexp](ll n, ll r) -> md {
        if (n < r) return 0;
        assert(0 <= n && n < static_cast<ll>(exp.size()));
        assert(0 <= r && r < static_cast<ll>(exp.size()));
        return exp[n] * dexp[r] * dexp[n-r];
    };
    // 区別するn 個のものから，繰り返し用いることを許して，r 個取り出して作った組
    // <-> 区別しないr 個のボールを，区別するn 個の箱に配る場合の総数
    auto nHr = [&nCr](ll n, ll r) -> md {
        return nCr(r+n-1, r);
    };
}

確率的なアレ¶

テスト中。

// n< 341550071728321 , ok.
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

typedef long long ll;

// will not overflow even if mod mod is too big..
ll mod_mul(ll a,ll b,ll mod){
    if(b == 0) return 0;
    ll res = mod_mul((a+a)%mod,b/2,mod);
    if(b % 2 == 1) res = (res + a)%mod;
    return res;
}

// use mod_mul if mod is too big.
ll mod_pow(ll x,ll n,ll mod){
    if(n==0) return 1;
    ll ret = mod_pow(x * x % mod,n/2,mod);
    if(n%2== 1)  ret = ret * x % mod;
    return ret;
}

// return probably prime.
//  if n < 341550071728321 return correct answer.
bool miller_rabin(ll n){
    if(n == 2) return true;
    if(n % 2 == 0 or n <= 1) return false;
    ll s=0,d=n-1;
    static const ll a[] = {2,3,5,7,11,13,17,n+1};

    while(d % 2 == 0){
        s++;d/=2;
    }
    for(int i=0;a[i]<n;i++){
        ll x = mod_pow(a[i],d,n);
        if(x != 1){
            ll r;
            for(r=0;r<s;r++){
                if(x == n-1) break;
                x = mod_mul(x,x,n);
            }
            if(r == s) return false;
        }
    }
    return true;
}

ll gcd(ll a,ll b){
    return b==0 ? a : gcd(b,a%b);
}

inline ll random(ll x,ll c,ll m){
    return (mod_mul(x,x,m)+c)%m;
}

vector<bool> sieve(const int M){
    vector<bool> isPrime(M);
    for(int i=2;i<M;i++) isPrime[i] = true;
    for(int i=2;i*i < M;i++){
        if(not isPrime[i]) continue;
        for(int j=i*i;j<M;j+=i){
            isPrime[j] = false;
        }
    }
    return isPrime;
}

// primes must contatin 2.
vector<ll> pollard_rho(ll n,const vector<ll>& primes,bool precheck=true){
    if(n == 0 or n == 1) return vector<ll>(1,n);
    vector<ll> ret;
    static const int constants[] = {1,51,73,0};
    if(precheck){
        for(size_t i=0;i<primes.size();i++){
            while(n % primes[i] == 0){
                n /= primes[i];
                ret.push_back(primes[i]);
            }
            if(n == 1) return ret;
        }
    }
    if(miller_rabin(n)){
        ret.push_back(n);
        return ret;
    }
    for(int i=0;constants[i];i++){
        ll x = 2,y = 2,d = 1;
        while(d == 1){
            x = random(x,constants[i],n);
            y = random(random(y,constants[i],n),constants[i],n);
            d = gcd(abs(x-y),n);
        }
        if(d == n) continue;

        vector<ll> dp = pollard_rho(d,primes,false);
        vector<ll> ndp = pollard_rho(n/d,primes,false);
        for(size_t j=0;j<dp.size();j++){
            ret.push_back(dp[j]);
        }
        for(size_t j=0;j<ndp.size();j++){
            ret.push_back(ndp[j]);
        }
        return ret;
    }
    return ret;
}


int main(){
    int N = 100;
    vector<bool> isprime = sieve(5);
    vector<ll> primes;
    for(int i=0;i<isprime.size();i++){
        if(isprime[i]){
            primes.push_back(i);
        }
    }
    for(int i=0;i<N;i++){
        //cout << i << " " << miller_rabin(i) << " " << endl;
        vector<ll> ps = pollard_rho(i,primes,true);
        for(int j=0;j<ps.size();j++){
            cout << ps[j] << " ";
        }
        cout << endl;
    }
    return 0;
}

乱数¶

XORShiftをつかったらうれしいかもしれない。

unsigned long xor128(){
    static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;
    unsigned long t;
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}

基数変換¶

vector<int> int_to_digits(int base,int N){
    vector<int> rev_ret;
    while(N != 0){
        rev_ret.push_back(N % base);
        N /= base;
    }
    reverse(all(rev_ret));
    return rev_ret;
}

int digits_to_int(int base,const vector<int>& digits){
    int ret = 0;
    for(int i=0;i<(int)digits.size();i++){
        ret *= base;
        ret += digits[i];
    }
    return ret;
}

二次方程式¶

// solve ax^2 + bx + c = 0.
vector<double> solve_quadratic_equation(double a,double b,double c){
    vector<double> ret;
    if(abs(a-EPS) < 0){
        ret.push_back(-c/b);
    }else{
        double d = b*b-4*a*c;
        if(d >= 0){
            if(b >= 0){
                ret.push_back(-2*c/(b+sqrt(d)));
                ret.push_back((-b-sqrt(d))/(2*a));
            }else{
                ret.push_back(-2*c/(b-sqrt(d)));
                ret.push_back((-b+sqrt(d))/(2*a));
            }
            if(abs(ret[0]-ret[1]) < EPS){
                ret.pop_back();
            }
        }
    }
    return ret;
}

nと互いに素な数の個数¶

EulerのTotient関数による。

// nと互いに素な数の個数
// n (1-1/p1) ... (1-1/pn)
template<typename T>
T totient(T n){
    if(n == 0) return 0;
    T ans = n;
    for(T x=2;x*x <= n;x++){
        if(n % x == 0){
            ans -= ans/x;
            while(n % x == 0) n /= x;
        }
    }
    if(n > 1){
        ans -= ans/n;
    }
    return ans;
}

幾何¶

#include <complex>
using namespace std;

const double EPS = 1e-9;
typedef complex<double> point;
typedef vector<point> vertex;

namespace std {
    bool operator<(const point &lhs, const point &rhs){
        if(real(lhs) == real(rhs)){
            return imag(lhs) < imag(rhs);
        }else{
            return real(lhs) < real(rhs);
        }
    }
}

// ベクタの長さ
double vector_length(point a){
    return abs(a);
}

// 二点間距離
double point_distance(point a,point b){
    return abs(a-b);
}

// 単位ベクトル
point unit_vector(point a){
    return a / abs(a);
}

// 法線ベクトル
pair<point,point> normal_vector(point a){
    point n1 = a * point(0,1);
    point n2 = a * point(0,-1);
    return make_pair(n1,n2);
}

// 点が一緒かどうか
bool point_eq(point a,point b){
    return abs(a-b) < EPS;
}
// 内積 (dot product) : a・b = |a||b|cosΘ
double dot(point a,point b){
    return real(conj(a)*b);
}

// 外積 (cross product) : |a×b| = |a||b|sinΘ
double cross(point a,point b){
    return imag(conj(a)*b);
}

// rotate by theta (-pi,+pi)
// arg will return angle of point
point rotate(point p,double theta){
    double cos_theta = cos(theta),
           sin_theta = sin(theta);
    return point(p.real()*cos_theta - p.imag()*sin_theta,
                 p.real()*sin_theta + p.imag()*cos_theta);
}
// 二つの円の共通接線 (たぶんおかしい)
//  http://www.ftext.org/text/subsubsection/1427
//  - 離れている → 4本
//  - 外接している → 3本
//  - 交わっている → 2本
//  - 内接している → 1本
//  - 含む → 0本
vector<pair<point,point> > common_tangent_of_two_circles(point c1,double r1,point c2,double r2){
    // shift c2 by c1.
    point nc = c2 - c1;
    // c.img == 0
    double ar = arg(nc);
    point c = rotate(nc,-ar);

    vector<pair<point,point> > ret;
    vector<double> xs;
    for(int sign=-1;sign<=1;sign+=2){
        xs.push_back((r1*r1 + sign*r1*r2)/c.real());
    }
    for(double x1 : xs){
        double y2 = r1*r1 - x1*x1;
        if(y2 < EPS) continue; // maybe two circle is crossed
        for(int sign=-1;sign<=1;sign+=2){
            double y1 = sign*sqrt(y2);
            // x1*x + y1*y = r1*r1 is such line
            point h = point(x1,y1);
            point t = point(x1+1,(-x1*(x1+1)+r1*r1)/y1); // tekito
            if(abs(y1) < EPS){
                t = point(x1,y1);
            }

            ret.push_back(make_pair(h,t));
        }
    }
    transform(ret.begin(),ret.end(),ret.begin(),[ar,c1](pair<point,point> p){
            return make_pair(rotate(p.first,ar)+c1,
                             rotate(p.second,ar)+c1);
        });
    return ret;
}


// a1,a2を端点とする線分(la)とb1,b2を端点(lb)とする線分の交差判定
bool is_intersected_linesegment(point a1,point a2,point b1,point b2){
    if(max(a1.real(),a2.real()) + EPS < min(b1.real(),b2.real())) return false;
    if(max(b1.real(),b2.real()) + EPS < min(a1.real(),a2.real())) return false;
    if(max(a1.imag(),a2.imag()) + EPS < min(b1.imag(),b2.imag())) return false;
    if(max(b1.imag(),b2.imag()) + EPS < min(a1.imag(),a2.imag())) return false;
    return (cross(a2-a1,b1-a1)*cross(a2-a1,b2-a1) < EPS) and
           (cross(b2-b1,a1-b1)*cross(b2-b1,a2-b1) < EPS);

}

// 与えられた4点が正方形を成すかどうかを判定する
// ただし，ほぼ同じ点が含まれたときにはfalseとする
//  4辺および対角線が1:1:1:1:sqrt(2):sqrt(2)ならばyes
bool is_square(vector<point> ps){
    // assert ps.size() == 4
    sort(ps.begin(),ps.end());
    vector<double> dist;
    for(int i=0;i<4;i++){
        for(int j=i+1;j<4;j++){
            // there is almost same point
            if(abs(ps[i] - ps[j]) < EPS) return false;
            dist.push_back(abs(ps[i]-ps[j]));
        }
    }
    sort(dist.begin(),dist.end());
    for(int i=0;i<4;i++){
        if(abs(dist[0] - dist[i]) > EPS) return false;
    }
    for(int i=4;i<6;i++){
        if(abs(sqrt(2)*dist[1]-dist[i]) > EPS) return false;
    }
    return true;
}


// a1,a2を端点とする線分(la)とb1,b2を端点とする線分(lb)の交点計算
point intersection_point_linesegment(point a1,point a2,point b1,point b2) {
    if(a1 == b1 or a1 == b2) return a1;
    if(a2 == b1 or a2 == b2) return a2;
    point b = b2-b1;
    double d1 = abs(cross(b, a1-b1));
    double d2 = abs(cross(b, a2-b1));
    double t = d1 / (d1 + d2);
    return a1 + (a2-a1) * t;
}

// 線分同士の最短距離
double dist_linesegment_and_linesegment(point a1,point a2,point b1,point b2){
    if(is_intersected_linesegment(a1,a2,b1,b2)){
        return 0;
    }
    return min(min(dist_linesegment_and_point(a1,a2,b1),
               dist_linesegment_and_point(a1,a2,b2)),
           min(dist_linesegment_and_point(b1,b2,a1),
               dist_linesegment_and_point(b1,b2,a2)));
}


// 2直線の直交判定 : a⊥b <=> dot(a, b) = 0
// みけんしょう
bool is_orthogonal(point a1,point a2,point b1,point b2) {
    return dot(a1-a2,b1-b2) < EPS;
}

// 2直線の平行判定 : a//b <=> cross(a, b) = 0
bool is_parallel(point a1,point a2,point b1,point b2) {
    return abs(cross(a2-a1,b2-b1)) < EPS;
}

// a1,a2を通る直線とb1,b2を通る直線の交差判定
bool is_intersected_line(point a1,point a2,point b1,point b2) {
    return not is_parallel(a1,a2,b1,b2);
}

// a1,a2を通る直線とb1,b2を通る直線の交点計算
point intersection_line(point a1,point a2,point b1,point b2) {
    point a = a2 - a1,b = b2 - b1;
    return a1 + a * cross(b, b1-a1) / cross(b, a);
}

// 直線と点との距離
double dist_line_and_point(point a1,point a2,point b){
    return abs(cross(a2-a1,b-a1)) / abs(a2-a1);
}

// 線分と点との距離
double dist_linesegment_and_point(point a1,point a2,point b){
    if(dot(a2-a1,b-a1) < EPS) return abs(b-a1);
    if(dot(a1-a2,b-a2) < EPS) return abs(b-a2);
    return dist_line_and_point(a1,a2,b);
}

// 直線と点の最短距離を実現する直線の点(すいせんの足)(みけんしょう)
point nearest_point_line_and_point(point a1,point a2,point b){
    return a1 + (a2-a1) * dot((a2-a1),(b-a1)) / norm(a2-a1);
}

// 線分と点の最短距離を実現する線分嬢の点(みけんしょう)
point nearest_point_linesegment_and_point(point a1,point a2,point b){
    if(dot(a2-a1,b-a1) < EPS) return a1;
    if(dot(a1-a2,b-a2) < EPS) return a2;
    return nearest_point_line_and_point(a1,a2,b);
}
// 円と線分の交差判定
bool is_cross_linesegment_and_circle(point c,double r,point a1,point a2){
    return (dist_linesegment_and_point(a1,a2,c) < r+EPS and
             (r < abs(c-a1) + EPS or r < abs(c-a2) + EPS));
}

// 点の進行方向
int ccw(point a,point b,point c){
    b -= a;c -= a;
    if(cross(b,c) > 0) return +1;    // counter clockwise
    if(cross(b,c) < 0) return -1;    // clockwise
    if(dot(b,c) < 0) return +2;      // c -- a -- b
    if(norm(b) < norm(c)) return -2; // a -- b -- c
    return 0;
}

// 点が真に多角形(凸?)の中にはいっているか
bool is_inner_point_vertex(const vector<point> &ps,point a){
    int cc = ccw(ps[0],ps[1],a);
    if(not(cc == 1 or cc == -1)) return false;
    for(size_t i=0;i<ps.size();i++){
        if(cc != ccw(ps[i],ps[(i+1)%ps.size()],a)) return false;
    }
    return true;
}

// 点が辺上、もしくは内部にある。(未検証)
bool is_inner_point_vertex_or_line(const vector<point> &ps,point a){
    for(size_t i=0;i<ps.size();i++){
        if(dist_linesegment_and_point(ps[i],ps[(i+1)%ps.size()],a) < EPS){
            return true;
        }
    }
    return is_inner_point_vertex(ps,a);
}

// 凸包 (UVA 109)
vector<point> convex_hull(vector<point> ps){
    int n = ps.size();
    int k = 0;
    sort(ps.begin(),ps.end());
    vector<point> ch(2*n);
    for(int i=0;i<n;ch[k++] = ps[i++]){
        while(k >= 2 and ccw(ch[k-2],ch[k-1],ps[i]) <= 0) --k;
    }
    for(int i=n-2,t=k+1;i>=0;ch[k++]=ps[i--]){
        while(k >= t and ccw(ch[k-2],ch[k-1],ps[i]) <= 0) --k;
    }
    ch.resize(k-1);
    return ch;
}

// remember,pts make convex.
// (http://judge.u-aizu.ac.jp/onlinejudge/cdescription.jsp?cid=ACAC002&pid=C)
double convex_diameter(const vector<point>& pts){
    const int n = pts.size();
    int is=0,js=0; // initial antipode.
    for(int i=1;i<n;i++){
        if(pts[i].imag() > pts[is].imag()) is = i;
        if(pts[i].imag() < pts[js].imag()) js = i;
    }
    double maxd = abs(pts[is]-pts[js]);
    int i,j,maxi,maxj;
    i = maxi = is;
    j = maxj = js;
    do{
        if(cross(pts[(i+1)%n]-pts[i],
                 pts[(j+1)%n]-pts[j]) >= 0){
            j = (j+1)%n;
        }else{
            i = (i+1)%n;
        }
        if(abs(pts[i]-pts[j]) > maxd){
            maxd = abs(pts[i]-pts[j]);
            maxi = i;maxj = j;
        }
    } while(not(i == is and j == js));
    // pts[maxi],pts[maxj] is pair of max diff.
    return maxd;
}

// 円と円の交点(2点ある前提)
vector<point> circles_point(point c1,double r1,point c2,double r2){
    double d = abs(c1-c2);
    double s = (r1+r2+d) / 2;
    double S = sqrt(s*(s-r1)*(s-r2)*(s-d));
    double h = 2 * S / d;
    point v = (c2-c1) / (abs(c2-c1));

    double m = sqrt(r1*r1 - h*h);

    vector<point> ret;
    ret.push_back(c1 + m*v+h*v*point(0,1));
    ret.push_back(c1 + m*v-h*v*point(0,1));
    return ret;
}

// clockwiseだと負
double triangle_area(point a,point b,point c){
    return cross(b-a,c-a)/2;
}

// clockwiseだと負
double vertex_area(vector<point> v){
    double ret = 0;
    for(int i=1;i<v.size()-1;i++){
        ret += triangle_area(v[0],v[i],v[i+1]);
    }
    return ret;
}

// 線分上の格子点をカウントする (端点も含む)
// ※ ただし端点も格子点であることを仮定
long long count_grid_point_on_line(const point& a, const point& b)
{
    assert(abs(a.imag() - round(a.imag())) < EPS && abs(a.real() - round(a.real())) < EPS);
    assert(abs(b.imag() - round(b.imag())) < EPS && abs(b.real() - round(b.real())) < EPS);
    return gcd(
        abs(static_cast<long long>(a.real()) - static_cast<long long>(b.real())),
        abs(static_cast<long long>(a.imag()) - static_cast<long long>(b.imag()))) + 1;
}

// エリアの中にある格子点 (x, y が両方とも整数である) を数える
// ピックの定理より、 面積 = *内部*の格子点の数 + 1/2 * 辺上の格子点 - 1 が成り立つ
// Ref: https://atcoder.jp/contests/typical90/tasks/typical90_ao
// この関数は 「内部の点 (端点や辺上の点を含まない)」と「端点を含む辺上の点」をペアにして返す
pair<long long, long long> count_grid_points_in_area(const vector<point>& ps)
{
    // ps の点は全て格子点である必要がある
    // assert(ps.size() >= 3);
    assert(all_of(ps.begin(), ps.end(), [](const point& p) {
        return abs(p.imag() - round(p.imag())) < EPS &&
            abs(p.real() - round(p.real())) < EPS;
    }));
    // FIXME: ここで単に point を使って面積を求めると誤差で死ぬケースが存在する
    // なので整数型で以降の計算を行う
    // 1. 面積を求める
    long long area2 = 0; // 面積の二倍を保持
    for (size_t i = 0; i < ps.size(); i++) {
        long long ax = ps[(i + 0) % ps.size()].real(), ay = ps[(i + 0) % ps.size()].imag();
        long long bx = ps[(i + 1) % ps.size()].real(), by = ps[(i + 1) % ps.size()].imag();
        area2 += (bx - ax) * (by + ay);
    }
    area2 = abs(area2); // 見る順番によってはマイナスになる可能性がある
    // 2. 辺上の点を求める
    long long grid_point_on_lines = 0;
    for (size_t i = 0; i < ps.size(); i++) {
        const auto& a = ps[i];
        const auto& b = ps[(i+1)%ps.size()];
        grid_point_on_lines += static_cast<long long>(count_grid_point_on_line(a, b));
    }
    // 重複してカウントされている端点分マイナスしておく
    grid_point_on_lines -= ps.size();

    // 面積の二倍を使って計算する場合
    // 2 * 面積 = 2 * 内部の格子点の数 + 辺上の格子点 - 2
    // => 内部の格子点の数 = (2 * 面積 - 辺上の格子点 + 2) / 2
    return { (area2 - grid_point_on_lines + 2) / 2, grid_point_on_lines };
}

ゲーム¶

Nim¶

いくつかのコインの山がある。この中からプレイヤーは山を一つ選び、1個以上の任意の数のコインを取る。最後のコインを取ったプレイヤーが勝ちである。この問題に対しては以下のことが知られている。すべての山のxorをとったとき、それが0であるとき、後攻の勝ち、それ以外のときは先攻の勝ち。これは次のことから言える。

全ての山が 0 の時、 xor は 0 である
現在の状態が xor 0 の時、次の操作で xor を 0 にすることはできない
現在の状態が xor 0 以外の時、次の操作で xor を 0 にすることができる

#include <iostream>
#include <vector>

using namespace std;

int main(){
    int xor_sum = 0;
    vector<int> coins = {1,2,4,5,6};
    for(size_t i=0;i<coins.size();i++){
        xor_sum = xor_sum ^ coins[i];
    }
    if(xor_sum != 0) cout << "First Player Win" << endl;
    else cout << "Second Player Win" << endl;
}

Grundy 数¶

Nim を拡張し、コインの山ではなくなんらかのゲームがいくつかあり、その一つを選び操作をする。例えば「A1, A2 , ... An が与えられ、どれかを選んでそれ自身ではないその約数に変更する。最初に操作が行えなくなった法の負け」 - https://atcoder.jp/contests/abc368/tasks/abc368_f

このような問題を解く際には一つのゲームの盤面を Grundy 数を使ってエンコードする。この Grundy 数は以下の性質を持つ。

その状態からどこにも遷移できない (負けの状態) のなら 0
遷移できる状態の Grundy 数を x1, x2, ... xn としたとき、それに含まれない最小の非負整数

いくつかあるゲームの Grundy 数を全て求めたらその xor 和が 0 なら初手が負ける、それ以外なら初手が勝てる。

そのほかの例: https://atcoder.jp/contests/typical90/tasks/typical90_ae

実装面では多くの場合は grundy 数を再帰的な関数か、 DP などで計算することになる。

// https://atcoder.jp/contests/abc368/tasks/abc368_f
// A -> A そのものではないその約数に変更する
int n; cin >> n;
vector<ll> a(n);
for (int i = 0; i < n; i++) cin >> a[i];
vector<int> grundy(*max_element(all(a)) + 1);
grundy[0] = grundy[1] = 0; // 遷移できない
for (int i = 2; i < grundy.size(); i++) {
    // 遷移できる先を計算する
    // 少なくともそれ自身で破ることで 1 (grundy 数 0) に遷移できる
    set<int> transition = { 0 };
    for (int c = 2; c*c <= i; c++) {
        // i が c で割れるので c に遷移できるし、 i/c にも遷移できる
        if (i % c == 0) {
            transition.insert(grundy[c]);
            transition.insert(grundy[i/c]);
        }
    }
    // 遷移先に現れない最小の値を計算する
    int r = 0;
    bool filled = false;
    for (auto t : transition) {
        if (t != r) {
            grundy[i] = r;
            filled = true;
            break;
        }
        r++;
    }
    if (!filled) {
        grundy[i] = r;
    }
}

// grundy 数の xor をとって、もし 0 なら初手が負ける
int grundy_xor_sum = 0;
for (int i = 0; i < n; i++) {
    grundy_xor_sum ^= grundy[a[i]];
}
bool first_lose = grundy_xor_sum == 0;

if (first_lose) {
    cout << "Bruno" << endl;
} else {
    cout << "Anna" << endl;
}

ゲーム木探索¶

TODO

なんか最適化とか¶

焼き鈍し¶

http://shindannin.hatenadiary.com/entry/20121224/1356364040

#include <chrono>
#include <iostream>
#include <algorithm>
using namespace std;
#define debug(x) #x << "=" << (x)
#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl

unsigned long xor128(){
    static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;
    unsigned long t;
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    const chrono::milliseconds annealing_time(2000);
    const auto start_time = chrono::system_clock::now();
    const auto end_time = start_time + annealing_time;
    const int precision_for_force_next = 1000000;
    auto current_time = start_time;

    const auto score = [](double x){
        return max(-pow(-x+100,2)+100,-pow(-x+333,2)+333);
    };

    auto current_state = 0.0;     // initial state
    auto current_score = score(current_state);
    auto best_state = current_state;
    auto best_score = current_score;

    // chrono::system_clock should NOT be called many times...
    while((current_time=chrono::system_clock::now()) < end_time){
        const auto progress = 1.0*(current_time-start_time)/(end_time-start_time);
        const bool force_next = progress*precision_for_force_next
                              < xor128()%precision_for_force_next;

        // [0..100] -> [-50..50] balance is IMPORTANT.
        const auto next_state = current_state + xor128()%101-50;
        const auto next_score = score(next_state);
        if(force_next or (next_score > current_score)){
            current_state = next_state;
            current_score = next_score;
            dump(current_state);
            dump(current_score);
        }else{
            // undo change.
        }

        if(current_score > best_score){
            best_score = current_score;
            best_state = current_state;
            dump(best_score);
        }
    }
    dump(best_score);
    dump(best_state);

    return 0;
}

いろんなデータ¶

階乗¶

$N$	$N!$
0	1
1	1
2	2
3	6
4	24
5	120
6	720
7	5040
8	40320
9	362880
10	3628800

数単位変換¶

N	日本語	英語
$10^{0}$	一	one
$10^{1}$	十	ten
$10^{2}$	百	hundred
$10^{3}$	千	thousand
$10^{4}$	万	ten thousand
$10^{5}$	十万	hundred thousand
$10^{6}$	百万	million
$10^{7}$	千万	ten million
$10^{8}$	億	hundred million
$10^{9}$	十億	billion
$10^{10}$	百億	ten billion
$10^{11}$	千億	hundred billion

bit¶

N	$2^N$	備考
0	1
1	2	boolの大きさ
2	4
3	8
4	16
7	128	charの最大値+1
8	256	unsigned charの最大値+1
16	65,534
24	16,777,216
31	2,147,483,648	intの最大値+1(about 2*10^9)
32	4,294,967,296	unsigned intの最大値+1
52	4,503,599,627,370,496	doubleのprecision
63	9,223,372,036,854,775,808	long longの最大値+1
64	18,446,744,073,709,551,616	unsigned long longの最大値+1(about 10^19)

最低限の設定ファイル¶

vim用¶

最低限のもの。ホームにおく。

set fileencoding=utf-8
set nocompatible
set t_Co=256
set ambiwidth=double

syntax on
filetype plugin indent on

set nobackup
set noswapfile

set completeopt=menuone
set wildmode=list:longest

set smartindent
set autoindent
set tabstop=4
set softtabstop=4
set shiftwidth=4

set smarttab
set expandtab

set incsearch ignorecase hlsearch
set showmatch
set wildmenu

set listchars=tab:>-,trail:-
set list

set backspace=indent,eol,start

アスキーコード¶

Char	Dec	Oct	Hex
(nul)	0	0000	0x00
(soh)	1	0001	0x01
(stx)	2	0002	0x02
(etx)	3	0003	0x03
(eot)	4	0004	0x04
(enq)	5	0005	0x05
(ack)	6	0006	0x06
(bel)	7	0007	0x07
(bs)	8	0010	0x08
(ht)	9	0011	0x09
(nl)	10	0012	0x0a
(vt)	11	0013	0x0b
(np)	12	0014	0x0c
(cr)	13	0015	0x0d
(so)	14	0016	0x0e
(si)	15	0017	0x0f
(dle)	16	0020	0x10
(dc1)	17	0021	0x11
(dc2)	18	0022	0x12
(dc3)	19	0023	0x13
(dc4)	20	0024	0x14
(nak)	21	0025	0x15
(syn)	22	0026	0x16
(etb)	23	0027	0x17
(can)	24	0030	0x18
(em)	25	0031	0x19
(sub)	26	0032	0x1a
(esc)	27	0033	0x1b
(fs)	28	0034	0x1c
(gs)	29	0035	0x1d
(rs)	30	0036	0x1e
(us)	31	0037	0x1f
(sp)	32	0040	0x20
!	33	0041	0x21
"	34	0042	0x22
#	35	0043	0x23
$	36	0044	0x24
%	37	0045	0x25
&	38	0046	0x26
'	39	0047	0x27
(	40	0050	0x28
)	41	0051	0x29
	42	0052	0x2a
	43	0053	0x2b
,	44	0054	0x2c
	45	0055	0x2d
.	46	0056	0x2e
/	47	0057	0x2f
0	48	0060	0x30
1	49	0061	0x31
2	50	0062	0x32
3	51	0063	0x33
4	52	0064	0x34
5	53	0065	0x35
6	54	0066	0x36
7	55	0067	0x37
8	56	0070	0x38
9	57	0071	0x39
:	58	0072	0x3a
;	59	0073	0x3b
<	60	0074	0x3c
=	61	0075	0x3d
>	62	0076	0x3e
?	63	0077	0x3f
@	64	0100	0x40
A	65	0101	0x41
B	66	0102	0x42
C	67	0103	0x43
D	68	0104	0x44
E	69	0105	0x45
F	70	0106	0x46
G	71	0107	0x47
H	72	0110	0x48
I	73	0111	0x49
J	74	0112	0x4a
K	75	0113	0x4b
L	76	0114	0x4c
M	77	0115	0x4d
N	78	0116	0x4e
O	79	0117	0x4f
P	80	0120	0x50
Q	81	0121	0x51
R	82	0122	0x52
S	83	0123	0x53
T	84	0124	0x54
U	85	0125	0x55
V	86	0126	0x56
W	87	0127	0x57
X	88	0130	0x58
Y	89	0131	0x59
Z	90	0132	0x5a
[	91	0133	0x5b
	92	0134	0x5c
]	93	0135	0x5d
^	94	0136	0x5e
_	95	0137	0x5f
`	96	0140	0x60
a	97	0141	0x61
b	98	0142	0x62
c	99	0143	0x63
d	100	0144	0x64
e	101	0145	0x65
f	102	0146	0x66
g	103	0147	0x67
h	104	0150	0x68
i	105	0151	0x69
j	106	0152	0x6a
k	107	0153	0x6b
l	108	0154	0x6c
m	109	0155	0x6d
n	110	0156	0x6e
o	111	0157	0x6f
p	112	0160	0x70
q	113	0161	0x71
r	114	0162	0x72
s	115	0163	0x73
t	116	0164	0x74
u	117	0165	0x75
v	118	0166	0x76
w	119	0167	0x77
x	120	0170	0x78
y	121	0171	0x79
z	122	0172	0x7a
{	123	0173	0x7b
	124	0174	0x7c
}	125	0175	0x7d
~	126	0176	0x7e
(del)	127	0177	0x7f

N	日本語	英語
\(10^{0}\)	一	one
\(10^{1}\)	十	ten
\(10^{2}\)	百	hundred
\(10^{3}\)	千	thousand
\(10^{4}\)	万	ten thousand
\(10^{5}\)	十万	hundred thousand
\(10^{6}\)	百万	million
\(10^{7}\)	千万	ten million
\(10^{8}\)	億	hundred million
\(10^{9}\)	十億	billion
\(10^{10}\)	百億	ten billion
\(10^{11}\)	千億	hundred billion