AtCoderで青になりました

　いつも記事を読んでくれてありがとうございます。
　先日（2020/2/22）AtCoderで開かれたABC156（AtCoder Beginner Contest 156 - AtCoder）でレート1608になり、いわゆる青コーダーになりました。ということで、本ポストでは私が競プロのコンテスト中に使うライブラリについて解説します。全部説明するのは大変ですから、コンテスト中に役に立っている気がするものだけ解説します。それらについても、特に工夫した点とかコンテスト中にどうやって使ってるかの紹介に留めます。詳しい解説は詳しく書いてる人の記事見てください。良いのがあったらリンク貼っときます。
　紹介したコードは自由に使ってもらって構いません。ソースコードは私のGitHubに公開されていますが、動作は保証していません。ソースコード内に問題へのリンクがあれば、基本的にその問題はACしてるって意味です。本ポスト内で紹介してるものについては動作確認済み。C++14でコンパイルしてください。
github.com

Segment Tree

　通称「セグ木」。数列に対する操作を要素数のlog時間くらいの計算量でやってくれて便利。概要はこの辺見てください。
tsutaj.hatenablog.com
beet-aizu.hatenablog.com

　私がコンテスト中によく使うのは「RMQ」「RSM」「RAQ」とかです。区間最小・最大、区間和、区間に加算など。工夫している点は特になし。それぞれのコードは
comprolib/rmqsegtree.cpp at master · wakuwinmail/comprolib · GitHub
comprolib/rsqsegtree.cpp at master · wakuwinmail/comprolib · GitHub
comprolib/raqlazysegtree.cpp at master · wakuwinmail/comprolib · GitHub
区間加算については「遅延セグ木」とかでググってください。

tsutaj.hatenablog.com
beet-aizu.hatenablog.com

　遅延セグ木が必要になるのは区間に対するクエリを処理しないといけない時だと思ってます。遅延セグ木に関しては行列積だったりと、妙なクエリとか代数的構造の処理を要求する問題に対応したい時のためにC++のtemplate機能を使ってます。
comprolib/lazysegtree.cpp at master · wakuwinmail/comprolib · GitHub
以下、気になる人はコードと並べて読んでみてください。（今後の更新のせいで）行数ズレたらメンゴ。投稿時点では合ってると思う。

7行目

template <typename T,typename E>

これはT：セグ木で管理する値の型、E：クエリで使う値の型、を柔軟に受け取れるやつです。これをやっておくと整数小数にコード変更無しで対応出来たり、行列を扱う時に楽です。すげー

15-18行目

F cal;//function for merge
G upd;//function for update
H ecal;//function for evaluate
P rcal;//function for range calculate

配列の要素をどうやって取り出すかとか、演算をどんな風に行うかを柔軟に指定できる。演算を変更してもコード全体を見直す必要がない。std::functionを使って実現してます。ただ、どの関数に何を指定すればいいのか使い始めてすぐの時はわからなくなってたから、実際に使ってみて慣れた方が良い。関数呼び出しが遅いと思うかもですが、定数倍で落ちたことないし多分大丈夫。

25,53行目

explicit LazySegmentTree(//配列で初期化する場合
explicit LazySegmentTree(//特定の要素で初期化する場合

空の配列から始めて一々updateするのはなんか遅そうなのと、めんどくさいことが多いので初期値をstd::vectorで受け取れるようにしました。それはそうとして、全部0とか特定の値で埋めたいこともあるので指定した値で埋められるようにもしてあります。どっちの場合も使うので両方あった方がコンテスト中に苦労しないと思います。
セグ木で工夫したのはこのくらい。実装とかは上に貼った記事見てください。ほぼ丸パクリしてるので。

Euler tour

　通称「オイラーツアー」。そこまで頻出な気はしてないけど、セグ木と組み合わせて使うことが多い。グラフ理論における「木」を列に変換するアルゴリズムです。変換された後の列において、部分木が連続した区間になります。すると、部分木に対するクエリがセグ木で処理できます。すごいですね！ただ、これは頂点に対するクエリしかできなくて、辺の集合を扱いたい場合はHL分解とかいうのが必要になってくるんですけど、持ってないので紹介しません。（調べたら辺属性のオイラーツアーっていうのがあるらしいんですけど知らないから🐛）良い感じの記事見つからなかったんでちょっとだけ詳しめに解説します。
comprolib/euler_tour.cpp at master · wakuwinmail/comprolib · GitHub

7行目

using Graph=std::vector<std::vector<int>>;//無向グラフ

いつも使ってるグラフ用のテンプレ。これ便利。

9行目

struct EulerTour{

構造体にしてます。オーバーヘッドとかで遅くなるかもしれないけど、各種操作が関数呼び出しでできるのは直接配列にアクセスするより直感的で良い。コンテスト中は考察に集中したい派なので添え字のめんどくさい所とかはライブラリに任せたい。

28-33行目

EulerTour(Graph g,int root=0):g(g){//initialize
    in.resize(g.size());
    out.resize(g.size());
    ord=0;
    dfs(root,-1);
}

コンストラクタです。基本的にdfsするだけ。

17-26行目

void dfs(int v,int p){//0-indexed
    in[v]=ord;
    ++ord;
    for(int c:g[v]){
        if(c==p)continue;
        dfs(c,v);
    }
    out[v]=ord;
    ++ord;
}

頂点に入った時と出た時にそれぞれ別の配列にordの値を記録する。ordの値は頂点に入る時、出るときに増える。これは配列にした時のindexを指してます。他はただのdfsでOK。

35-37行目

std::pair<int,int> interval(int n){
    return std::make_pair(in[n],out[n]);
}

dfsした時に調整してあるからこれだけでOK。閉区間になってるところだけ注意（セグ木は半開で揃えてある）。別に半開で返してもいいだけど、意味的に閉区間の方が合ってそうだからそうしてる。

169-170行目

auto t=et.interval(v);
st.update(t.first,t.second+1,1);

セグ木と一緒に使う時はこんな感じで終点を+1して使う。
どっかに図入りの良い資料あった気がするんだよな。知ってる人いたらコメントださい。

強連結成分分解

　有向グラフをDAGにできる。DAGだと何かと嬉しいので便利。ただ、強連結成分分解（は長いですから、SCCと表記します）やれって言われることそんなに多くないような・・・最初から入力をDAGにしてくれてる人もよくいるらしいと聞く。アルゴリズムの概要とか正当性の証明は貼った記事見て。（SlideShareの方にはLCAとかHL分解も載ってて良い感じ。本ポストではどっちも解説しませんので。）
mathtrain.jp

www.slideshare.net

8行目

using Graph=std::vector<std::vector<int>>;//無向グラフ

無向グラフ用のテンプレかと思いきや、有向グラフでもOK。[元の頂点][先の頂点]だけ更新するだけ。無向は両方向更新する。

10行目

struct StronglyConnectedComponents{

構造体にしておくとdfsするのにグローバル変数使ったり引数を長くしたりしなくて済む。速度面でどうなってるのかは知らないけど、まあ大したことないでしょ。関数呼び出しのオーバヘッドの方がどうせボトルネックなんじゃないですか？（知らんけど）（詳しい人いたら教えてください）

15行目

Graph RES;//強連結成分に分けた後のグラフ

これ大事だよ。

53-60行目

RES.resize(n);
for(int i = 0; i < n; ++i){
    for(auto v: G[i]){
        int x=gp[i],y=gp[v];
        if(x==y)continue;
        RES[x].push_back(y);
    }
}

RESはこれだけで作れる。Gは元のグラフ、gpは強連結成分の番号。元のグラフで辺があって、強連結成分が異なるなら辺を張る。多重辺が嫌な人は適当に頑張ってください。そこまでめんどくさくはないです。std::setで管理すればオーダーもlog倍くらいで済む。

63-69行目

void dfs(int pos){//帰りがけ順の記録
    visited[pos]=true;
    for(auto i:G[pos]){
        if(!visited[i])dfs(i);
    }
    porder.emplace_back(pos);
}

帰りがけ順っていうのはこういうことです。AOJのライブラリ問題に木の探索する問題があったはず。ちなみにこのSCCはAOJのSCCって問題を通してます。実装はらずひるどめもをパクってたりそうじゃなかったりするので見といてください。

ei1333.github.io

Union Find(disjont set)

素集合データ構造とかいうやつ。集合の各要素が高々1つのグループに入っている状態を管理できる。UnionFind木を使うと、グループの併合・異なる要素が同一グループに属するかの判定が要素数に対してだいたいlog時間で出来る。便利。chokudaiスライドが有名かつわかりやすいと思います。

www.slideshare.net

　私のライブラリには3つのunionfindと名の付くソースコードがありますね。どれも基本的な機能はUnionFindです。素のunionfin.cppはchokudaiスライドの丸写しですから、解説は省略します。以下、chokudaiスライド相当の知識は前提にします。知らなかったら読んでね。

pertial persistent

comprolib/pertial_persistent_unionfind.cpp at master · wakuwinmail/comprolib · GitHub

「pertial persistent」は部分永続なUnionFindです。永続データ構造の詳しい解説はしません。ググってでてきたのを見ましょう*1*2*3
出来ることは、それぞれの併合クエリを処理する前の状態を参照すること。ただし、部分永続なので過去に戻って操作をやり直すことは出来ない。永続が必要ならstd::vectorじゃなくて永続配列を使ってUnionFinを実装すればそれで終わり。

10行目

std::vector<int> time;

これが重要です。通常のUnionFindとの違いはこれだけ。配列の中身は、「最後に親が変更された時刻」です。chokudaiスライドで紹介されているテクニックのうち、経路圧縮を用いない場合を考えると、親は高々1回しか変更されないことが分かりますね？それが変更された時刻を記録します。

19行目

time[i]=std::numeric_limits<int>::max();

初期値は十分に大きい数にしておく。これは、親が変更されていないことを表してます。

24-27行目

int root(int t,int x){
    if(time[x]>t)return x;
    else return root(t,par[x]);
}

引数のtの意味は、「時刻tにおける状態を参照する」という意味。親は高々1回しか変更されないので変更された時刻以前であれば自身が親、それ以外の場合は通常通り探索を行う。

38-42行目

if(size[x]<size[y]){
    par[x]=y;
    time[x]=t;
    size[y]+=size[x];
}

unite内の処理です。chokudaiスライドではrankになってるけど、sizeで比較しても同様の効果が得られます。証明気になる人は各自考えてみてください。確かiwiさんのブログ記事であったよなーって思って「マージテク」とかで検索かけたらはてなグループ日記がサービス終了してた。気になる人は「データ構造をマージする一般的なテク」とかが関係あるので調べてみてね。*4
timeの更新に関して、sizeの小さい方が大きい方に吸収されるので併合のタイミングで書き換えておきます。

50-60行目

int whenMerged(int t,int x,int y){//0-indexed
    if(!same(t,x,y))return -2;
    int ng=-1,ok=t,mid;
    mid=(ng+ok)/2;
    while(ok-ng!=1){
        mid=(ng+ok)/2;
        if(same(mid,x,y))ok=mid;
        else ng=mid;
    }
    return ok;
}

読んで関数名のごとく。時刻tまでの操作において、xとyが同じグループに属した時刻を調べる。当然、単調性があるので二分探索ができる。実際にsameクエリに投げて確認する。計算量はO(log T log N)くらいで済む。Tは時刻、Nは要素数。時刻に関しては座圧っぽいことをやれば高々Nに抑えられるのでO((log N)^2)くらいになると思う。今回verifyに使った問題ではTがNまでしかないので工夫する必要は無かったので何もしてないです。

Potentialized

comprolib/potentialized_unionfind.cpp at master · wakuwinmail/comprolib · GitHub

　ポテンシャル付きUnionFind、重み付きとも呼ばれる。要素間の差を管理できる。UnionFindっぽくないと思いきやガッツリ応用になっている。2要素の差がたくさん与えられたとき、それを併合操作だと思うと同じ集合に属する2要素間の差は原理的に計算できるはず。これを要素数に対してlog時間くらいで計算できます。
qiita.com


7行目

template<typename T>

ポテンシャルの値の型はテンプレート機能で柔軟に受け取ります。+と-が定義されていればOK。ただし、交換法則は成り立ってないと困る。まあ問題を解く上では性質の良いものしか与えられないことが多いからそんなに気にしないです。

13行目

std::vector<T> diff;

diffには、「rootとのポテンシャル差」が格納されています。

23-32行目

int root(int x){
    if(par[x]==x){
        return x;
    }
    else{
        int r=root(par[x]);
        diff[x]+=diff[par[x]];
        return par[x]=r;
    }
}

普通のUnionFindと異なる点は、再帰的に親を探索する際にdiffの値を更新していく部分です。自身の親のdiffは再帰が終わった時点で確定してるので正しくdiffが求まります。経路圧縮をしているので何度もrootを実行してもdiffの値は増えません（再帰が1回行われるんですが、根になっているノードのdiffは0なので値は変化しません）。

34-41行目

T potential(int x){
    root(x);
    return diff[x];
}

T potential_diff(int x,int y){
    return potential(y)-potential(x);
}

ポテンシャルの差は、各要素のポテンシャルをrootを実行することで計算した後、その差を返すという処理になっています。これだと、xとyが異なる集合に属する場合に正しくないんですけど、C++14ではstd::optionalみたいなやつが使えないんですよね。後で書きますが関数呼び出しの時点でチェックしてます。

47-53行目

void unite(int x,int y,T w){//0-indexed,potential(y)-potential(x)=w
        w+=potential(x);
        w-=potential(y);

        x=root(x);
        y=root(y);
        if(x==y)return;

コメントに書いてる通り、yのポテンシャルがxのポテンシャルよりもw大きいというクエリを処理します。wに足したり引いたりしてるのは、UnionFindの都合上xとyじゃなくてそれぞれの根を接続するので重み調整です。調整を行ったので以降、x,yは通常のUnionFindと同様にx,yの根に置き換えて処理を行います。この時点でxとyが等しい場合、ポテンシャルを上書きしてしまうことになるので処理を終えます。この部分で整合性チェックとかする設計でもいいと思います。

55-59行目

if(size[x]<size[y]){
    par[x]=y;
    size[y]+=size[x];
    diff[x]=-w;
}

sizeが小さい方を大きい方の下に繋げる。部分永続のと同じ。diff[x]=-wはクエリの通り。ここで直感的になるように調節しておいたのだよ。

85-86行目

if(pu.same(x,y))std::cout<<pu.potential_diff(x,y)<<std::endl;
else std::cout<<"UNKNOWN"<<std::endl;

こんな感じでpotential_diffを使う時はsameかチェックしないと変なことが起こる。

UnionFindシリーズはこんな感じです。AtCoderでは良く出てるイメージがある。グラフのある部分が連結かどうかの判定にも使えるので、とりあえず貼るみたいな場面は多いと思います。

Rolling hash

　通称「ロリハ」、正式っぽい名称は「ラビン-カープ文字列検索アルゴリズム」らしい長くてウザいですね。前計算O(|S|)してhashを計算しておくことで文字列の比較（一致するかどうか）を定数時間で行える。連続した部分文字列の比較についても定数時間で行える。え、やばくねこれ？

www.slideshare.net

9行目

struct HashString{

例のごとく、という感じで構造体にしています。ハッシュの計算をコンストラクタで、取得をメンバ関数で行えるのでまあ分かりやすさ重視です。

11-12行目

using hashtype=long long;
using hashcomptype=std::tuple<hashtype,hashtype,hashtype>;

いちいち書いてたら長すぎてコードが読めなくなったのでやってます。hashcomptypeなんですけど、とりあえずhashtypeが3つ入ったtupleになっています。ハッシュの衝突とかに関してはもう少し後で書きますが、3つで比較してると結構遅いのでTLが危なくなったら2つで比較すると良いかも。この辺は柔軟に対応できるように設計したいんだけど、今まで困ったことがないのでこれ使っています。

14-19行目

std::vector<hashtype> table1;
...
std::vector<hashtype> powtable1;
..

tableは計算したhashの値を格納するための配列、powtableは基数の累乗を格納するための配列。累乗の方ってコンパイル時計算できると思うんだけど、よく知らないです。

21-24行目

hashtype mod=1000000007ll;
hashtype base1=1919ll;
hashtype base2=114514ll;
hashtype base3=334334ll;

hashを作るための数字たち。modはなるべく大きめの素数、ただし掛け算したときにlong longに入りきらないと面倒なのでそのくらいの大きさまでで。baseは基数とか呼ばれる数字で、任意。ただしこっちも掛け算したときにlong longに入りきる程度が良い。基数1とか0とか流石に良くないはず。実行時にランダムで決めてもOK。Codeforcesとかhackがある時は固定だともしかしたら危ないのかも。3つとって時間が間に合うならhackは気にするほどじゃないとは思う。

28行目

table1.resize(s.size()+1);

サイズを1つ多く取ってます。s[i]に対応する値がtable[i+1]に入ってます。これはhashの取得の時に半開区間で取れて便利にするための工夫です。この後もそのつもりで見てください。

42-56行目

for(int i=0;i < s.size(); ++i){
    table1[i+1]=table1[i]*base1+s[i];
    ...
    powtable1[i+1]=powtable1[i]*base1;
    ...
    table1[i+1]=table1[i+1]%mod;
    ...
    powtable1[i+1]=powtable1[i+1]%mod;
    ...
}

ここはhashを計算してmod取ってるところです。hashの計算方法はwikipediaの該当箇所が詳しいです。実用上ではどんな計算してるのか気にする必要ないです。何をしているかというと、1つ前のインデックスまでのhash値にbaseを掛けたものに、今のインデックスに対応する文字列中の文字の値を足します。こうして文字列の値を反映しつつ、連続した区間の値を取り出せるようになってるわけです。
ラビン-カープ文字列検索アルゴリズム - Wikipedia

59-65行目

bool operator == (HashString s){
    if(this->table1.size()!=s.table1.size())return false;
    if(this->table1.back()!=s.table1.back())return false;
    ...
    return true;//maybe :)
}

HashString同士の比較用です。実際には一部のhashを取得してそれを比較することの方が多いと思うのであんまり使わないかも。普通に比較するのと計算量変わらんし。やってることは文字列長を比較した後に文字列全体のhashを比較。

67行目

hashcomptype get_hash(int l,int r){//[l,r)でのhash

これがメインの関数。よく使う。半開区間だと色々とやりやすいことが多いです。コンテスト中はどうせ添え字気にしたりとかしないんで閉区間で取ってもいいっちゃいいです。

70-72行目

hashtype hash1=table1[r]-(powtable1[r-l]*table1[l]%mod);
    while(hash1<0)hash1+=mod;
    hash1=hash1%mod;

これが大事です。右端の値（s[l-1]までのhash値）から、左端までの値（s[l-1]までのhash値）にbaseのr-l乗を掛けた値を引きます。 table \left[ r \right]=table[l ] \times pow(base,r-l)+(sの[l,r)に関するhash)であることを考えれば分かりますね。

83行目

return std::make_tuple(hash1,hash2,hash3);

tupleにして返します。これは==演算子で比較できるので比較してみて、一致してれば文字列も一致と考えます。

　hashと聞けばすぐに思い浮かぶ衝突についてですが、競プロの世界ではパフォーマンス重視なんでhashが一致したときに元の文字列を確認して判断とかしません（諸説あり）。そもそも、自分のロリハは元の文字列を保存してません。つまり、ある程度の確率で衝突はするけど、しないことを期待します。どの程度の確率で衝突が起こるのか？とか、ハッシュを衝突させる方法*5とかは詳しい資料があります。気になる人は読んでみて、それでも良いって人はロリハを使いましょう。
snuke.hatenablog.com

qiita.com

Mod Pでの演算

数え上げ系の問題ではおなじみの「Mod P」で答えを求めるやつです。以下、Pは素数とします。合成数が持たない性質とかも利用してます。分かる範囲で注釈しておきます。

mod付累乗

comprolib/modpow.cpp at master · wakuwinmail/comprolib · GitHub
2の1000000005乗のような計算をMod Pで計算できます。計算時間は累乗の指数のlog程度。繰り返し自乗法とか呼ばれてるテクニックらしい。Pythonの組み込み関数のpowではデフォルトで行われているらしく、速い（らしい）。C++ではそういうのはできないので自分で実装するしかないです。
satanic0258.hatenablog.com
qiita.com


4行目

long long mod=1000000007LL;

modの法はグローバルかどこかに宣言しておくと便利だと思います。どうせ変更しないことがほとんどでしょ。他のライブラリでmodの宣言が省略されてる場合があるんですけど、基本はグローバルにこれがあるものと思ってください。

6行目

template<typename T>

たぶん64bit整数しか使わないから要らないかも。unsignedを使ってもめんどくさくないようになってます。

7行目

constexpr T modpow(T a,int n){//(a^n)%MOD

constexprを付けておくとコンパイル時に値が決まってる範囲で勝手に計算をしておいてくれるらしいです。付け得。

9-15行目

while(n>0){
    if((n&1)!=0){//n%2==1
        ret=ret*a%mod;
    }
    a=a*a%mod;
    n=n/2;
}

ループの度にnを2で割っていき、aは2乗していく。これにより、k回目（0-indexed）のループ開始時にはaは、最初のaの[(2のk乗)乗]になっています。nが奇数の時だけ答えにaを掛けているのは、nを2進数として表した時に1の桁だけ使うイメージ。再帰でやる実相もあるらしいけど、たぶんこっちの方が速いと思う。

積の逆元

comprolib/modinv.cpp at master · wakuwinmail/comprolib · GitHub

まず、「フェルマーの小定理」を知らない人がいたら検索してください。 x^{p-1} \equiv 1(mod p)という式の左辺を x \times x^{p-2} と分解してみると、xにx^{p-2} を掛けた値はmod pにおいて1になることが分かります。xに掛けると1になる数のことを積の逆元といいます。注意点として、フェルマーの小定理は法が素数の時だけ成り立ちます。なので、素数modじゃないときは使わないでね。あんまりハマったことはないですが、xがpの倍数であるときにはxの逆元は定義できません。形式的に計算はできますが、結果は0であり意味がありません。
mathtrain.jp


って、思ってたんですけどいま確認したら自分のライブラリは拡張ユークリッドの互除法とか呼ばれるテクニックを使うやつでした*6。これは、xとmが互いに素な場合に x \times y \equiv 1(mod m)であるyが計算できるそうです。簡単にやってることを書いておくと、 x \times y - a \times mb = 1であるようなyを求めています。右辺を1にするためにxとmの最大公約数が1でないといけないので、このような制約が付きます。素数じゃなくても使えるので便利ですね。こっちの方はよくわからないのでコメントなしです。勉強したらまた何か書きます。

組み合わせ（コンビネーション）

comprolib/modcomb.cpp at master · wakuwinmail/comprolib · GitHub
組み合わせとは、「n個の中からm個を選ぶ通り数（ただし、それぞれの物は区別しない）」です。以下、この組み合わせの通り数を \tbinom{n}{m}と表記します。これは、中学校・高校で習った記憶がある人もいるかもしれません。これを、前計算O(n log n)かけることで各組み合わせについて定数時間で求めます。 \tbinom{n}{m} = \frac{n!}{(n-k)!k!}ですので、積の逆元を利用して実際にこれを計算します。計算に必要な数は、「nまでの階乗」とそれらの逆元です。逆元計算の計算量ですけど、拡張ユークリッドの互除法の方はよくわからないですが、フェルマーの小定理と繰り返し自乗法を併せて用いるとlog nにはなるので上はこれで抑えられてます。

1-2行目

std::vector<long long> fac;
std::vector<long long> ifac;

階乗とその逆元を入れるようの配列です。配列名はfactorialから取ってるんですけど、なるべく発音しないでください。構造体にすればよかったのに、なぜか関数で書いてるせいで配列がグローバルにあります。は？？？（青になった報告のブログなんで修正はしないけど、そのうち書き直します。）

7行目

T modcomb(T a, T b){

本体。実際に組み合わせを計算する用の関数。中身は定義に基づいて計算をするだけ。

14行目

void combinit(int maxn){

前計算用の関数。使いそうな最大のnを渡す。

15-16行目

fac.resize(maxn+1);
ifac.resize(maxn+1);

nの階乗について、fac[n]で取り出せるようにしたいので配列サイズは+1で取る*7。

17-22行目

fac[0] = 1;
ifac[0] = 1;
for(long long i = 0; i<maxn; i++){
    fac[i+1] = fac[i]*(i+1) % mod;
    ifac[i+1] = ifac[i]*modpow(i+1, mod-2) % mod;
}

実際に計算している部分です。0の階乗を1として、順番にmaxn+1まで計算します。上で紹介したように、フェルマーの小定理と繰り返し自乗法で逆元を計算してます。ループ変数をintで取るとオーバーフローするのでコンテスト中に実装すると絶対にミスる。必携ライブラリ間違いなしですね。

modint

comprolib/modint.cpp at master · wakuwinmail/comprolib · GitHub

これまでに紹介した関数たちを過去の存在にする最強の構造体です。普通の整数と同じように演算できる上に、計算結果は常にmod pで得られる。抽象化した遅延セグ木と組み合わせるなどの使い方が便利。
noshi91.hatenablog.com

4行目

template<std::uint64_t Mod> struct ModInt{

modをテンプレートで受け取るようにしてある。コンストラクタで受け取っても変わらんと思った人、後でわかる、この便利さ。

5行目

using u64=uint64_t;

内部で持つ値は符号無し64bitということにしてある。工夫すれば32bitで持って高速化みたいなことが出来そうだけど、まあそんなに急ぐ必要はない。先に言っておくと、高速化の工夫はあまりしてないです。剰余が重いのでなるべく剰余を取らないようにするとか、逆元計算とか、色々工夫できる点はありそう。必要になったらまたやります。

12-13行目

constexpr u64 &value(){return a;}
constexpr const u64 &value()const{return a;}

いわゆるgetメソッドみたいなもの。内部で持つ値はprivate指定してあるので、関数で取り出す。constにしてるやつとしてないやつが両方あるのは、なんかこうしないとうまくいかないから。四則演算のoperatorを定義しているんですが、そこで右辺値をconstとして受け取ったりする関係です。詳しい説明はできないので、知ってる人がいたら教えてください。

16-22行目

constexpr ModInt &operator+=(const ModInt rhs){
    a+=rhs.value();
    if(a>=Mod){
        a-=Mod;
    }
    return *this;
}

64bit整数で値を管理してるんで、とりあえずオーバーフローとか気にせずに和を取る。よくある1e9+7を法としてる場合には32bit整数で持ってもギリギリオーバーフローはしない。普段なら剰余を取って済ませるところですが、ifと引き算で済ませる。これが本当に速いのかは知りませんけど、内部がこうなってることで実用上困ることはない。正しく演算できてれば問題ないって話です。
引き算もほとんど同じ。ifを計算前に評価するのに気を付ける（符号無し整数で管理しているため）。

32-35行目

constexpr ModInt &operator*=(const ModInt rhs){
    a=a*rhs.value()%Mod;
    return *this;
}

積は諦めて剰余取ってます。割り算して引いた方が速いとかあるのかな？知りません。

37-47行目

constexpr ModInt &operator/=(ModInt rhs){
    u64 n=Mod-2;
    while(n>0){
    if((n&1)!=0){//n%2==1
        *this*=rhs;
    }
    rhs*=rhs;
    n=n/2;
    }
    return *this;
}

割り算は、フェルマーの小定理から逆元を実際に計算してます。サボりまくって遅くなってそうだけど、これしか分からず。速い実装知ってる人がいたら教えてください。

49-63行目

constexpr ModInt operator+(const ModInt rhs){
    return ModInt(*this)+=rhs;
}
...

この辺、C++に詳しくないので説明うまくできないんですけど、2項演算の結果として計算結果を返すことで1+2+3みたいなのが正しく計算できるみたいです。四則演算はぜんぶこの形。

64-78行目

constexpr ModInt &operator++(){
    ++a;
    if(a>=Mod){
        a-=Mod;
    }
    return *this;
}
...

インクリメント・デクリメントの処理。1足す1引くの処理を書いてるだけ。この辺はあると使う、というよりは使いたくなることが多そうだから先手必勝しておいた。

80-83行目

constexpr ModInt operator- (){
    if(a==0)return ModInt(0);
    else return ModInt(Mod-a);
}

単項の-演算子。これは間違いなく使う。日常的に使ってる組み込み型の偉大さの片鱗を感じるライブラリになったと思う。

88行目

using Mi=ModInt<1000000007ull>;

これがやりたかったんですよ。usingを使うとtemplateまで含めてalias（というのか知らんけど）を作れるので、コンストラクタで渡すよりも書きやすい。だからtemplateを使う必要があったんですね*8。

90-92行目

Mi operator"" _mi(unsigned long long n){
    return Mi(n);
}

これを定義しておくと、定数のmodintが114514_miみたいな形で書けるようになる。わざわざコンストラクトの形にしなくて良いので楽。

他にも、cinとcoutに対応させるやつとか累乗を計算するやつとかを作ることもできます。色々カスタマイズして最強のmodintを創り上げよう！！！

行列

comprolib/matrix.cpp at master · wakuwinmail/comprolib · GitHub
名前の通り、行列を扱う構造体。一般に半環を要素に持てば行列としての演算も半環になるらしいです。ちなみに、このライブラリは半環では不味くて、環を乗せないと壊れる可能性があります。それはさておき、環とか半環ってなんじゃ？？？ってなっている人にざっくり説明しておくと、環は「足し算引き算かけ算」ができる代数的構造で、半環は「足し算かけ算」ができる代数的構造です。和については交換法則、積については分配法則・結合法則が成り立っているものを考えることが多い。詳しくは代数学を学んでください。私もよく知りません。
このライブラリに関してはまだまだ整備中で、自分でも解説できないところがあるので、解説が無い部分は各自勉強してみてください。読者への課題です（まるなげ）。

19行目

template<typename T>

いつもの。modintを乗せられるのがかなり便利。

21行目

using scalarvalue=T;

scalarvalueの方がわかりやすい。

28-31行目

constexpr Matrix(std::vector<std::vector<scalarvalue>> data):data(data){
    col=data.size();
    row=data[0].size();
}

2次元配列を受け取ってそのままコピーするコンストラクト。行と列の数を知ってると何かと便利なので覚えておく。

33-35行目

constexpr Matrix(int col,int row,scalarvalue init=0):row(row),col(col){
    data.assign(col,std::vector<scalarvalue>(row,init));
}

大きさと要素の初期値を受け取る。最初は全部0で埋まってる時とかには使う。

37-59行目

constexpr Matrix &operator+=(const Matrix rhs){
    for(int i = 0; i < col; ++i){
        for(int j = 0; j < row; ++j){
            data[i][j]+=rhs[i][j];
        }
    }
}
constexpr Matrix &operator-=(const Matrix rhs){
...
constexpr Matrix &operator*=(const scalarvalue rhs){//スカラー倍
...

和は行列サイズ分だけ時間を掛けて足していくだけ。+=が定義されている必要がある。スカラー倍も同じく。
問題は差の方。半環は引き算ができないって言ったんですが、環ならできます。乗せてるscalarvalueに-=が定義されていない時に引き算をすると壊れます。
念の為、具体例を出しておきましょう。
a+b=min(a,b)という演算に関して、a-bとか意味が分からないですね。-bをbに関する和の逆元とみなすと、和の単位元が∞なのでb+(-b)=min(b,-b)=∞である必要があります。もちろん、このようなbは（実数上には）存在しないのでminを取る演算の逆元は定義できません。よって、引き算も定義できないことになります。このような時に行列の引き算をすると壊れる、というわけです。まともな頭してればそんな関数呼び出ししないと思うから大丈夫。

61-71行目

constexpr Matrix operator*(Matrix rhs){
    Matrix<scalarvalue> ret(col,rhs.get_row(),0);
    for(int i = 0; i < col; ++i){
        for(int k = 0; k < rhs.get_row(); ++k){
            for(int j = 0; j < row; ++j){
                ret[i][j]+=data[i][k]*rhs[k][j];
            }
        }
    }
    return ret;
}

行列積です。愚直に3乗のアルゴリズム。色々調べると、シュトラッセンのアルゴリズム*9というものがあって、だいたいO(n^2.807)くらいで計算できるらしい。めんどくさいのでやってません。代わりって程効果があるか知りませんが、ループ順を工夫してます*10。vectorというのは連続したメモリ領域に確保されているらしい*11ので、近い要素を行ったり来たりするのは速いです。それを利用できているかもしれないし、できていないかもしれません。

77-92行目

constexpr Matrix pow(long long n){
    assert(is_square());//正方形じゃないと計算が定義できない
    Matrix<scalarvalue> ret(col,row,0);
    Matrix<scalarvalue> a(*this);
    for(int i = 0; i < col; ++i){
        ret[i][i]=1;//単位行列を作る
    }
    while(n>0){
        if((n&1)!=0){//n%2==1
            ret=ret*a;
        }
        a=a*a;
        n=n/2;
    }
    return ret;
}

これは何回か紹介してる繰り返し自乗法を使って高速化してます。計算量はO(N^3 log n)くらい。

行列は、3乗オーダーの計算がよく出てくるので2次元vectorじゃなくて1次元に展開したvectorか配列を使う方が良いかもしれないです。TLEギリギリになったら考えるかも...


ここまでお付き合い頂きありがとうございました。5回くらい読み直してますが、誤字脱字や不正確な部分があればご指摘お願いします。はてなブログのコメントかTwitter(@wakuwinmail)までお願いします。余談ですが、このブログ書いてる間にAtCoderのレートは水色になりました。

それでは、今回の記事で見つけたライブラリの穴を埋めなきゃいけないんでこの辺で失礼します。今後も応援よろしくお願いします。

*1:https://www.slideshare.net/qnighy/re-15443018

*2:http://noshi91.hatenablog.com/entry/2019/02/04/175100

*3:http://www.kmonos.net/pub/Presen/PFDS.pdf

*4:https://tomorinao.blogspot.com/2018/10/blog-post_39.html

*5:http://hos.ac/blog/#blog0003

*6:https://drken1215.hatenablog.com/entry/2018/06/08/210000

*7:https://twitter.com/beet_aizu/status/1146400552844050432?s=20

*8:https://dic.nicovideo.jp/a/%E3%83%A1%E3%82%AC%E3%83%88%E3%83%B3%E3%82%B3%E3%82%A4%E3%83%B3

*9:https://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A5%E3%83%88%E3%83%A9%E3%83%83%E3%82%BB%E3%83%B3%E3%81%AE%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0

*10:https://qiita.com/Syo_pr/items/92b3cf7d7fc5dab4a3a7

*11:http://crane-house.sakura.ne.jp/program/vector.html