ABC 186 E - Throne
一次合同式についての教育的な良問だったので、考えたことを書き残しておく。
問題
考えたこと
この問題は言い換えると、初期値がsの状態からkを何度か足すことで、nの倍数を作れるか、という問題になる。
つまり、(kx+s)%n==0となる最小のxを求める問題（解が存在しない場合も有り得る）であり、この式を更に変形すると、kx+s≡0 (mod n) => k*x≡-s (mod n) という一次合同式を解く問題に帰着できる。
一次合同式 a*x≡b (mod m) の解を求める手順は以下の通り。
- 整数 a, b, mの最大公約数gを求めて、gでそれぞれの値を割った値に置き換える
- これはただの前処理で、gで割る前と後の合同式は同値
- 次に、a, mが互いに素であるかを確認するためにa, mの最大公約数を求める
- このとき、その最大公約数が1でないなら、解は存在しない（これは乗法逆元が定まらないため）
- 解が存在するなら、mod mにおけるaの乗法逆元a'とすると、その解は、a'*bとなる
注意点
ここに丁寧に記されている通り、乗法逆元の存在する条件は、aとmが互いに素であること。
以下の2点は正であるのに誤解しがちなところので注意する。
上記を踏まえると逆元が存在する場合には、与えられるa, b, mによってはaとmが互いに素ではあるが、mは素数ではないという場合が有り得る。
さらに整理すると、合成数である可能性もあるMが与えられ解として計算結果のmod Mを求めよという問題で計算結果がWAとなるのは、計算過程に除法があるならばmod Mの乗法逆元を求める必要があり、割る数とMが互いに素でない場合（フェルマーの小定理等では）その逆元を正しく求められないからである。
この問題もその一例で、(m-2)乗で逆元が求まらない場合がある。結局この場合の求め方がわからなかったので、以下の解答コードではACLのmodintを使いお茶を濁している。（誰か教えてください）
解答コード
#include <bits/stdc++.h> using namespace std; #include <atcoder/all> using namespace atcoder; int solve(int a, int b, int m) { int g = gcd(gcd(a, b), m); a /= g, b /= g, m /= g; if (gcd(a, m) != 1) return -1; modint::set_mod(m); modint ma = a; modint mb = b; return (mb * ma.inv()).val(); } signed main() { int t; cin >> t; while (t--) { int n, s, k; cin >> n >> s >> k; cout << solve(k, -s, n) << endl; } }
個人的なメモ
gcdと__gcdの挙動は異なる
※この問題で__gcdを使って実装するとバグります。（__gcdはパラメータに負数を含む場合に結果が負数になる可能性があるため）
__gcd() と gcd(), 実質一緒だと思ったら負数に対する結果が違うらしく嵌まりかけた
— si💊 (@iiljj) 2020年5月2日
実験の結果，戻り値は
・__gcd() → 0でなく絶対値の小さい引数と同符号
・gcd() → 常に非負整数
どうしてこうなった
__gcd は rotate の実装のために置かれてる内部用の関数で、我々が使うことを意図して置かれてるものではないので
— えびちゃん (@rsk0315_h4x) 2020年5月2日
結論gcdを使う。
- ACL modintを呼ぶ前に
modint::set_modを呼ぶのを忘れない
参考
合成数の可能性があるMが与えられて、計算結果にmod Mを求める必要があり、かつ、計算過程に逆元を求める必要がある問題
DFSによるオイラーツアー順序を求める実装
高難度帯の部分問題として出てくるので実装を置いておく。
実装
以下の実装は、訪れた順に頂点を採番していく。
/// オイラーツアー struct EulerianTrail{ vector<int> order; EulerianTrail(vector<vector<int>>& G, int root=0){ int n=G.size(); order.resize(n); int cur=0; function<void(int,int)> dfs=[&](int v, int p){ order[v]=cur; cur++; for(auto nv: G[v]){ if(nv==p) continue; dfs(nv, v); } }; dfs(root, -1); } };
トポロジカルソートの実装をライブラリ化した
過去問埋めしてて、トポロジカルソートで解く問題に当たって面白かった↓。
良い機会なので、改めて実装についてまとめておこうと思う。
前提
トポロジカルソートとは、閉路無し有向グラフ(DAG)において、どの頂点もその出力辺の先の頂点より前にくるように並べること。
ようするに、そのグラフの頂点を「実行順序があるタスク」と見た時、先に実行すべきタスクの順に並べることを言う。
- 各頂点の入次数をカウントしておく
- 閉路無し有向グラフなら必ず入次数0の頂点があるので、queue に push する
- 以下の処理を queue が空になるまで繰り返す
- queue に入っている頂点を順に取り出し（pop）、今見ている頂点から移動可能な頂点の入次数を -1 する
- 入次数が0になった頂点は都度 queue に push する
こうして queue に push された頂点の順序がトポロジカル順序になる
検証用問題
実装
/** * トポロジカルソート * 頂点数をV、辺数をEとして O(V+E) */ vector<int> topological_sort(const vector<vector<int>>& G){ const int n=G.size(); vector<int> in_cnt(n); // 各頂点の入次数 // 各頂点の入次数をカウント for(int v=0; v<n; v++){ for(auto nv: G[v]){ in_cnt[nv]++; } } // 入次数が0の頂点を始点としてqueueにpush（同時にresにもpush） vector<int> res; queue<int> q; for(int v=0; v<n; v++){ if(in_cnt[v]==0){ q.push(v); res.push_back(v); } } // queueに入っている頂点から移動可能な頂点の入次数をデクリメントし、入次数が0になったものをqueueにpush while(!q.empty()){ int v=q.front(); q.pop(); for(auto nv: G[v]){ in_cnt[nv]--; if(in_cnt[nv]==0){ q.push(nv); res.push_back(nv); } } } // reverse(res.begin(), res.end()); // 入次数0の頂点を先頭にしたい場合 return res; }
個人的なメモ
頂点 i => j (i < j)という様に辺が与えられて、2点間から計算されるスコアの最高値を求めよ、みたいな問題（上に載せたABC 188 E Peddler がまさにそう）はトポロジカルソートで解く。
この手の問題を dfs で解こうとすると大概TLEになる。（たまにやらかす）それはそうで、1つ頂点が複数の親を持つ可能性があるので、各2点間におけるスコアを求めるには少なくとも入次数0である全ての頂点を始点として dfs する必要がある。つまり頂点数をV、辺数をEとすると計算量は、O( V * (V + E) )になる。
与えられるグラフが木であるなら dfs で問題ないが、上の様な条件では木にならないことに注意する。（木なら各頂点の親の個数は1個以下になる。また辺の数が V-1 個になるので、まず辺数の制約を確認する）
エラトステネスの篩の実装
エラトステネスの篩と、それを利用した素因数の高速列挙の実装をライブラリ化した。（どちらかというと後者の実装を残す目的で作った）
前提
高速素因数分解
- エラトステネスの篩でnまでの整数を素数判定する過程で、各値の最小の素因数をテーブル
min_prime_factorに保存しておく x / min_prime_factor[x]でxの約数が得られるので、値xの最小の素因数を確認すると同時に、xをその約数に上書きする- 上記をxが1になるまでの間繰り返すことで素因数を列挙できる
- ループ回数は高々log x回なので、O(log x)
検証用問題
実装
/// エラトステネスの篩 struct Eratosthenes { vector<bool> is_prime; /// nまでの素数一覧 vector<int> primes; /// 整数xの最小の素因数 vector<int> min_prime_factor; Eratosthenes(const int n) { min_prime_factor.resize(n + 1); /// nまでの素数一覧 is_prime.assign(n + 1, true); is_prime[0] = is_prime[1] = false; min_prime_factor[0] = min_prime_factor[1] = 1; for (int i = 2; i <= n; i++) { if (is_prime[i]) { for (int j = 2 * i; j <= n; j += i) { is_prime[j] = false; if (min_prime_factor[j] == 0) min_prime_factor[j] = i; } primes.push_back(i); min_prime_factor[i] = i; } } } /** * xの素因数の種類を取得 * @param x n以下の数値を指定する必要がある */ set<int> get_prime_factors(int x) { set<int> s; while (x > 1) { s.insert(min_prime_factor[x]); x /= min_prime_factor[x]; } return s; } /** * xの素因数とその数を取得 * @param x n以下の数値を指定する必要がある */ map<int, int> get_prime_factor_and_cnt(int x) { map<int, int> mp; while (x > 1) { mp[min_prime_factor[x]]++; x /= min_prime_factor[x]; } return mp; } };
Trie木の実装をライブラリ化した
前提
- Trie木は、複数の単語（文字列）を登録でき、ある文字列（またはそのprefix）が登録済みであるかを高速に検索できるデータ構造
- 基本的な機能はシンプルに、
insert/searchの2つだけ - 実態は有向木
- 基本的な動作は、単語のprefixが同じならNodeを共有し、異なる文字が現れた際に新規にNodeを追加していく
- prefixを共有しながら枝分かれしていくイメージ
- 各Nodeは3つの基本情報を持つ
next: 次のNodeへのポインタ- 扱いたい文字種の数
char_size(例: a-zなら26)だけ、ポインタを持つ必要がある
- 扱いたい文字種の数
commom: そのノードがいくつの単語で共有されているか- この情報を持たせることで、そのNode以降にいくつの単語が登録済みかがわかる
- 故に、Node(0) を確認すればTrie木に登録済みの単語の総数がわかる
accept: そのノードがある単語の末尾文字であるか- 単語検索の際には
acceptを確認し、prefix検索の際には確認しない
- 単語検索の際には
- 今Trie木が管理している文字数を知りたい場合は、単にNodeの数をカウントすればよい（Node(0)は文字ではないので、正確には
size - 1になる）
Trie木のわかりやすい解説
検証用問題
※ 以下の実装を使うには、メソッドのシグネチャとtemplate部分を変更が必要
実装
/** * Trie-Tree * * @param char_size Trie木で扱う文字種数 * @param base Trieで扱う文字種のうちの先頭 */ template<int char_size, int base> class Trie{ private: struct Node{ int common; // この頂点を共有する単語数 vector<int> next; // 次の頂点id vector<int> accept; // この頂点を終端とする単語idを保持する Node(): common(0), next(vector<int>(char_size, -1)) {}; }; vector<Node> nodes; public: Trie(): nodes(vector<Node>(1)) {}; void insert(const string& s, int string_id=0){ nodes[0].common++; int last=0; for(int i=0; i<s.size(); i++){ int char_id=(int)(s[i]-base); if(nodes[last].next[char_id]==-1){ nodes.push_back(Node()); nodes[last].next[char_id]=nodes.size()-1; last=nodes.size()-1; }else{ last=nodes[last].next[char_id]; } nodes[last].common++; } nodes[last].accept.push_back(string_id); } bool search(const string& s){ int last=0; for(int i=0; i<s.size(); i++){ int char_id=(int)(s[i]-base); if(nodes[last].next[char_id]==-1) return false; else last=nodes[last].next[char_id]; } return nodes[last].accept.size()>0 ? true : false; } bool search_prefix(const string& s){ int last=0; for(int i=0; i<s.size(); i++){ int char_id=(int)(s[i]-base); if(nodes[last].next[char_id]==-1) return false; else last=nodes[last].next[char_id]; } return true; } int count_words(){ return nodes[0].common; } int size(){ return nodes.size(); } };
三分探索の実装をライブラリ化した
二分探索と比べて、区間の内分点の取り方やら、下に凸か上に凸かやら、地味に考えることが多いので汎用性を目指してライブラリ化してみた。
前提
- 三分探索とは、ある区間[l, r]において、極値がただ一つだけ存在するとき、その極値の近似値を区間の長さをnとして、およそ O(log n) で求める手法
- 二分探索が単調グラフにおける最大値、最小値を求めるものに対して、三分探索は凸関数における最大値、最小値を求めるためのもの
- 三分探索の実装では、以下のことを行っている
三分探索のわかりやすい解説
実装
l, rに、対象の区間を指定する- 第3引数
fに、対象の関数を与える - 上に凸な関数なら、第4引数
is_downward_convexに false を指定する
/// 三分探索 double ternary_search(double l, double r, const function<double(double)> f, const bool is_downward_convex=true){ int loop=500; if(is_downward_convex){ // 下に凸 while(loop--){ auto mid_l=(l*2+r)/3, mid_r=(l+2*r)/3; if(f(mid_l)>=f(mid_r)) l=mid_l; else r=mid_r; } }else{ // 上に凸 while(loop--){ auto mid_l=(l*2+r)/3, mid_r=(l+2*r)/3; if(f(mid_l)<=f(mid_r)) l=mid_l; else r=mid_r; } } return l; }
検証用問題
呼び出し元
上記問題を例に、以下に実際の使い方を載せる。
double p; /* ternary_search */ signed main() { cin>>p; auto f=[&](double x){ return x+p/pow(2,x/1.5); }; auto ans_x=ternary_search(0,1e18,f); printf("%.15lf\n", f(ans_x)); }
浮動小数を文字列として受け取り、整数型に変換する実装（個人的なメモ）
これは何
最近コンテストで、入力値が小数第4位程度まで与えられ、その入力値を素直にdoubleで受け取って計算すると丸め誤差で WA となるような問題が散見されるようになった。
（直近の問題だと ABC 191 D Circle Lattice Points がそう）
この手の問題の対策はいくつかあると思う。
- 多倍長浮動小数型がある言語ならそれを使う
long double等、longより精度の高い型 + EPSを使う- 入力値を10p 倍して整数に変換して計算した後、10p で割って最終的な解を求める
以下の実装は、3における入力値を10pするためのもので、丸め誤差が生じないように文字列として入力値を受け取り、それを long long 型で返すものである。（毎回実装考える時間がもったいないので、実装を残しておく）
実装
/** * @param s 文字列化した小数値 * @param p sの数値を 10^p 乗倍する * * 10^p乗倍した結果、小数点以下が0になることと、その値がlong long型に収まる必要がある */ long long string_double_to_long(const string & s, const int p){ if(s.find('.')==string::npos){ long long p=1; while(p--) p*=10; return stoll(s)*p; }else{ int n=s.size(); for(int i=0; i<n; i++){ if(s[i]=='.'){ int tail_len=n-i-1; int loop_cnt=p-tail_len; auto t=s.substr(0,i)+s.substr(i+1,n); while(loop_cnt--) t.push_back('0'); return stoll(t); } } } }
