ARC 118 A - Tax Included Price 解説
問題
消費税T%のとき課税後の整数値として現れないN番目のものを出力せよ。
制約
解説
ある自然数Aがあって、T*A の結果、100の位へ繰り上げがあるなら、(100+T)*A/100-1 は登場しない整数値である。
なぜなら、T*A で100の位への繰り上げが生じるならば、 の差が2になるためである。
また、100の位への繰り上げが発生する A と、その次に繰り上げが発生する A' の差は、高々100なので（T=1の場合、例えば、A=100, A'=200 と差が100になる）、登場しない整数値のN番目を愚直に求めようとするならば計算量は O(N9*100) になる。
A=1 から始めて、 の値が始めて0になった時の値を A の値を B、1～B の範囲で登場しなかった整数値の個数を cnt とすると、B+1～2B、2B+1～3B、... における個数も同様に cnt 個になる。（この B、cnt を求める計算量は高々 O(100) で十分高速）
cnt が求まっているなら cnt 倍数番目の登場しない整数値を O(1) で求めることができ、例えば 2*cnt 番目は、 と求めることができる。よって、この問題を高速に解くためには、N/cnt 番目の登場しない整数値を求め、そこから更に、N%cnt 番目の登場しない整数値を求めればよい。
N%cnt 番目の登場しない整数値が と表せるとすると、求める解は
となる。
実装
#include <bits/stdc++.h> using namespace std; using ll = long long; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int t,n; cin>>t>>n; int b=1,cnt=1; int sum=0; // 100の位への桁上がり検査のため、t*bのmod100の値を持つ while(t*b%100){ sum+=t; if(sum>=100){ cnt++; sum-=100; } b++; } int d=0; if(b>1) { d=n/cnt; n%=cnt; } int c=0; sum=0; if(n) while(n--){ while(sum<100){ sum+=t; c++; } sum-=100; } cout<<((100+t)*((ll)d*b+c))/100-1<<endl; }
感想
区間内にある条件を満たす要素の個数を高速に数え上げる典型手法に、ある周期ではその要素数が同じになる → 周期内、周期の余り部分で別々に数え上げるというものがあり、その良い事例だった。
ARC 118 C - Coprime Set 解説
問題
以下の条件を満たす長さNの数列Aを構築せよ。
- 数列のどの2つの要素も互いに素ではない
解説
まず、どの2つの要素も互いに素ではなく、一方で全体のgcd は1である、という点に着眼し、数列の要素を3つのグループに分けることを考える。
- グループ1：2*3の倍数でかつ5の倍数ではない
- グループ2：3*5の倍数でかつ2の倍数ではない
- グループ3：5*2の倍数でかつ3の倍数ではない
このように3つのグループを作ると、どのグループの要素も他のグループの要素と1以外の約数を持ち、かつ全体の gcd は1にできる。
それぞれのグループの要素の最大個数をカウントしてみると以下の様になる。
- グループ1：10000/6-10000/30=1666-333=1333
- グループ2：10000/15-10000/30=666-333=333
- グループ3：10000/10-10000/30=1000-333=667
- 合計：2333
合計が最大ケースに対して167個足りないとわかる。（初見時はここで詰まった。）
ここで構築した数列Aを観察してみると、グループ1～3の要素が1つずつ既に存在しているなら、2*3*5の倍数をAに追加しても問題ないことに気付く。（全体の gcd が1であることは保たれるし、2*3*5の倍数の要素は他のどの要素とも1以外の約数を持つ）
グループ4を2*3*5の倍数の要素とすれば、その最大個数は333個となるので、これで2500個の最大ケースにも対応することができた。
（素数7を増やすことでグループを増やしたくなるが、この場合は逆に作れる要素数が減少してしまう。）
実装
#include <bits/stdc++.h> using namespace std; #define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i) template<class t> using vc=vector<t>; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int n; cin>>n; vc<int> a; a.push_back(6); a.push_back(15); a.push_back(10); int now=2, x=6; while(now*x<=10000){ if(now%5==0){ now++; continue; } a.push_back(now*x); now++; } now=2,x=15; while(now*x<=10000){ if(now%2==0){ now++; continue; } a.push_back(now*x); now++; } now=2, x=10; while(now*x<=10000){ if(now%3==0){ now++; continue; } a.push_back(now*x); now++; } // グループ4 a.push_back(30); now=2, x=30; while(now*x<=10000){ a.push_back(now*x); now++; } rep(i,n){ cout<<a[i]<<endl; } }
ARC 118 B - Village of M People 解説
問題
Nが のようにK分割されている。Mが与えられるので同様にK分割し、その結果を整数列Bとしたとき、
が最小となるようなものを1つ出力せよ。
制約
解説
数列Bの要素として実数を許すなら、 とすることで、
を全て0にすることができる。
しかし、Bは整数列である必要があるので、上記の差が最小になるように各要素を切り下げ、または切り上げ、合計がMになるようにする必要がある。
一旦、全てのB_iは切り捨てることにすると となる。これらの総和を sum とすると、M-sum 個が切り上げる必要がある要素数になる。
切り上げるべき要素は、
の小数部分がより大きいものを切り上げた方が得なので、各要素を
の値でソートし、大きいものから M-sum 個を切り上げる。
実装
#include <bits/stdc++.h> using namespace std; #define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i) using ll = long long; template<class t> using vc=vector<t>; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int k,n,m; cin>>k>>n>>m; vc<ll> a(k); rep(i,k) cin>>a[i]; vc<pair<int,int>> p(k); rep(i,k){ p[i]={a[i]*m%n,i}; } sort(p.rbegin(),p.rend()); vc<int> b(k); int sum=0; rep(i,k){ int z,h; tie(z,h)=p[i]; int x=a[h]*m/n; b[h]=x; sum+=b[h]; } int rem=m-sum; rep(i,rem){ int z,h; tie(z,h)=p[i]; b[h]++; } rep(i,k) cout<<b[i]<<endl; }
ABC 200 D - Happy Birthday! 2 をDPで解く
問題
N個の要素からなる集合A の2つの異なる部分集合で和の mod 200 の値が等しいものが存在するかを判定し、存在する場合はその2つの集合を出力せよ。
制約
公式解説
公式解説は、鳩ノ巣原理により の場合は必ず解が存在することから、
の範囲において部分集合を全探索すればよい、というもの。
後にDPによる解法を示すが、明らかにこちらの解法の方が実装上楽である。（DP解法では、テーブルからの数列の復元や、Nが大きい場合に扱う数値が64bit整数型すらオーバーフローする可能性がある等、注意すべき点が多い）
の場合は必ず解が存在する理由は、N要素からなる集合から作れる部分集合の個数は空集合を除いて
通りなので、
の時で255通り。この時点で200通りを超えるため、
の範囲において mod 200 の値が重複する部分集合が必ず存在すると示せる。（ここが鳩ノ巣原理）
実装上は、 の範囲においてbit全探索を行い、各部分集合が mod 200 のどの値に対応するかを保持し、mod 200 のある値が2通り以上の部分集合で作れるならば、解が存在するとしてその2つの集合を出力すればよい。
実装
#include <bits/stdc++.h> using namespace std; #define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i) #define irep(i, m, n) for (int i = (int)(m); i < (int)(n); ++i) template<class t> using vc=vector<t>; signed main() { int n; cin >> n; vc<int> a(n); rep(i, n) cin >> a[i], a[i] %= 200; n = min(8, n); vc<vc<vc<int>>> p(200); irep(s, 1, 1 << n) { int sum = 0; vc<int> tmp; rep(i, n) if ((s >> i) & 1) { (sum += a[i]) %= 200; tmp.push_back(i); } p[sum].push_back(tmp); } rep(i, 200) { if (p[i].size() >= 2) { cout << "Yes" << endl; cout << p[i][0].size() << ' '; for (auto e : p[i][0]) cout << e + 1 << ' '; cout << endl; cout << p[i][1].size() << ' '; for (auto e : p[i][1]) cout << e + 1 << ' '; cout << endl; return 0; } } cout << "No" << endl; }
DP解法
以下の2通りについて解の存在判定の仕方を場合分けする。
- Aの全体集合を除く部分集合で和のmod 200が0になるようなものが存在する場合
- 1.以外の場合
この様に場合分けするのは、上記の2通りの場合で、部分集合を復元する手順が異なるため。（ 1) の場合の集合を 2) 場合の手順で復元しようとするとバグる）
それぞれの方法を以下に示す。
1. 部分集合に 0(mod 200) が存在する場合
（簡単のため、以下に例示する集合Aの要素は全てmod 200を取った値とする）
例として以下の様なものがある。
全体集合ではない部分集合の和で 0(mod 200) を作れるならば、その補集合から適当に1要素を取ってきて、Bはその1要素だけ、Cはその1要素と0を作る部分集合とすれば必ず解を作ることができる。
つまり、上の例からB、Cは以下の様になる。
全体集合を除外しているのは、上述の通り 0(mod 200) を作る集合以外に少なくとも1要素は存在しないと、2つの異なる部分集合を作れないため。
2. 1) 以外の場合
1) 以外の場合で解が存在するなら、集合B、Cの総和のmod 200の値は0以外であり、かつその部分集合に 0(mod 200) となるようなものは存在しない。
例として以下の様なもの。
この場合は以下の様にDPテーブルを定義することで、遷移元を辿ってB、Cを復元することができる。
上記から今 dp[i][m] が1以上で、dp[i-1][(m-a[i]+200)%200] も1以上なら総和が m(mod 200) となる部分集合の1つに要素 a[i] を含む者が存在する、ということがわかる。
DPテーブルを更新する過程で、あるテーブルの値が2以上になった瞬間、その m(mod 200) の値を作る部分集合が2つ存在し、かつ、その2つの集合は今見ている a[i] を含むものと、含まないもので2つ存在しているとわかる。
よって実装上は、ある dp[i+1][m] が2以上になったら更新を打ち切り、i+1から遡って集合Bを復元（要素a[i]を含む集合）、iから遡って集合Cを復元（要素a[i]を含まない集合）することになる。
上述した2通り以上になったら更新を打ち切ることはこの実装上の注意点となる。
数え上げDPをしていることになるので、Nが大きく、各要素の mod 200 の値がバラバラである場合などは、最後まで更新しようとするとLong型でもオーバーフローする可能性がある。（N=200なら、最悪2200に近い値を扱うことになる）
それを確認するには、以下のコードを実行してみると良い。
#include <bits/stdc++.h> using namespace std; #define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i) using ll = long long; template<class t> using vc=vector<t>; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); vc<int> a(200); iota(a.begin(), a.end(),0); vc<vc<ll>> dp(201, vc<ll>(200)); dp[0][0]=1; rep(i,200){ rep(j,200){ dp[i+1][j]+=dp[i][j]; dp[i+1][(j+a[i])%200]+=dp[i][j]; } } }
runtime error: signed integer overflow: 5902957071859434052 + 5902959422909828636 cannot be represented in type 'long long int'
実装
#include <bits/stdc++.h> using namespace std; #define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i) #define irep(i, m, n) for (int i = (int)(m); i < (int)(n); ++i) #define all(v) v.begin(), v.end() #define sz(x) (int)(x.size()) template<class t> using vc=vector<t>; signed main() { cin.tie( 0 ); ios::sync_with_stdio( false ); int n; cin>>n; vc<int> a(n); int all=0; rep(i,n){ cin>>a[i]; a[i]%=200; (all+=a[i])%=200; } bool allZero=0; if(all==0) allZero=1; vc<vc<int>> dp(n+1, vc<int>(200)); dp[0][0]=1; rep(i,n){ rep(j,200){ int nj=(j+a[i])%200; dp[i+1][j]+=dp[i][j]; dp[i+1][nj]+=dp[i][j]; } if(dp[i+1][0]>1){ if(i<n-1 or !allZero){ // 部分集合で0を作れる vc<int> pZero; // 0(mod200)を作る部分集合 set<int> s; int pre=0; for(int j=i; j>=0; j--){ int tmp=(pre-a[j]+200)%200; if(dp[j][tmp]){ pZero.push_back(j+1); s.insert(j); pre=tmp; } } int one=-1; rep(j,n) if(s.count(j)==0) { one=j+1; break; } cout<<"Yes"<<endl; cout<<1<<' '<<one<<endl; pZero.push_back(one); sort(all(pZero)); cout<<pZero.size()<<' '; rep(j,pZero.size()) cout<<pZero[j]<<' '; cout<<endl; return 0; } } irep(j,1,200) if(dp[i+1][j]>1){ // 部分集合0を使わずにB,Cを作れる vc<int> b,c; int pre=j; for(int k=i; k>=0; k--){ int tmp=(pre-a[k]+200)%200; if(dp[k][tmp]){ b.push_back(k+1); pre=tmp; } } pre=j; for(int k=i-1; k>=0; k--){ int tmp=(pre-a[k]+200)%200; if(dp[k][tmp]){ c.push_back(k+1); pre=tmp; } } reverse(all(b)); reverse(all(c)); cout<<"Yes"<<endl; cout<<b.size()<<' '; rep(k,b.size()) cout<<b[k]<<' '; cout<<endl; cout<<c.size()<<' '; rep(k,c.size()) cout<<c[k]<<' '; cout<<endl; return 0; } } cout<<"No"<<endl; }
感想
400点問題にしては難しすぎると思ったけど、自身が難しく捉えているだけだった。
燃やす埋める問題についてまとめ
燃やす埋める系の問題をまとめて解いたので、学んだ点について書き残しておく。
燃やす埋める問題とは以下の様な問題である。
いくつかのものが与えられて、それを赤か青に塗り分ける必要がある（燃やすか埋めるか）。それぞれがどちらの色に属するかによって得られるスコアが異なる。
最適に色を塗り分けた時の、スコアを最大値を求めよ。
更に条件として、あるものが赤に属して別のあるものが青に属する場合は罰金、もしくは、あるものが赤に属して別のあるものも赤に属する場合は利得が得られるという依存関係がいくつも与えられる。
この様な問題は解き方が定型化されており、最大流問題に帰着し、適切なグラフを構築してその最小カットを求めることで解くことが可能である。
- いくつかの注意点
- 例題
- ARC 085 E - MUL
- ABC 010 D - 浮気予防
- TopCoder SRM 594 Div I FoxAndGo3
- RUPC 2019 Day 2 H - Ghost
- Educational Codeforces Round 55 G. Petya and Graph
- RUCPC 2019 Day 1 G イルミネーション
- Google Code Jam 2008 World Finals - The Year of Code Jam
- 参考
いくつかの注意点
1. 扱えない条件がある
あるものが赤で別のあるものも赤であるなら罰金、または、あるものが赤で別のあるものが青であるなら利得という依存関係は最小カットで求めることは困難である。一見似た条件に見えるが、これらは最大カットを求める問題になる。（参考）
2. 構築するグラフに負辺は含まれてはいけない
燃やす埋める系問題では、利得や罰金の値が負で与えられる場合がある。負辺があるグラフの最小カット問題はNP困難問題となる（参考に載せた診断人さんのスライドP.22を参照）。そのため条件に負の値が与えられたときは、適宜言い換えを行う必要がある。燃やす埋める系の問題ではよくこの言い換え求められる場面が多く、以下に載せた例題の解法にも使われている。
例題
ARC 085 E - MUL
問題
N個の宝石が与えられてそれぞれの価値が である。そのぞれの宝石について、残すか割るかを選ぶことができ、残した場合は
を、割った場合は0円の得る。さらに、
を割ったならば、
の倍数の宝石も必ず割る必要がある。得られる金額の最大値を求めよ。
制約
解説
この手の利得が負数の場合があり、それを辺のキャパシティとしてそのまま使えないという場合は、一旦全ての利得が得られるものと仮定して、その内の得られないもの最小カットで求めて、罰金として仮定した利得から減算する、という解き方がフレームワークとして使える。
この問題の場合 であるものを利得として全て得られると仮定する。
最小カットを求めるグラフは、必ず残すに属する頂点 S と、必ず割るに属する頂点 T と、N個の宝石の S - (n頂点) - Tというグラフを作り、宝石iを残すなら 、割るなら
の罰金があるとして辺を張る。
宝石iを割るならiの倍数の宝石も必ず割る、という条件は言い換えると、宝石iが割るに属して、iの倍数の宝石のどれかが残すに属する、ということは許されないということになる。
「許されない」という依存関係を表したい場合は、そのようにグラフをカットした場合、それが必ず最小カットにはならないようにすればよく、この場合は頂点iからiの倍数の頂点に無限に相当するcapacityの辺を張ることで表現できる。
上記のグラフのS-T間の最小カットが仮定した利得にたいする罰金の最小値となるので、答えは（仮定した利得）-（グラフの最小カット）となる。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
#define irep(i, m, n) for (int i = (int)(m); i < (int)(n); ++i)
using ll = long long;
template<class t> using vc=vector<t>;
const ll INF = 1e9+10;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
signed main()
{
cin.tie( 0 ); ios::sync_with_stdio( false );
int n; cin>>n;
vc<ll> a(n);
ll ans=0;
rep(i,n) {
cin>>a[i];
ans+=max(a[i], 0LL);
}
auto flow=Dinic<ll>(n+2);
int s=n, t=n+1;
rep(i,n){
flow.add_edge(s,i,max(-a[i],0LL));
flow.add_edge(i,t,max(a[i],0LL));
for(int j=i+(i+1); j<n; j+=(i+1)){
flow.add_edge(i,j,INF);
}
}
ans-=flow.get_max_flow(s,t);
cout<<ans<<endl;
}
ABC 010 D - 浮気予防
これは参考に載せたyosupoさんの記事にも挙げられている問題。
問題
N頂点E辺のグラフがあって、頂点0からアクセスできないようにしたい頂点がG個ある。G個の頂点にアクセスできないようにするためにできる操作が、①任意の辺を1つ取り除く、②任意の頂点を1つ無効化する（接続される辺は維持される）の2つあり、いずれもコスト1かかる。G個全てアクセスできないようにするためのコストの最小値を求めよ。
制約
解説
必ず頂点0側に属する頂点をS(=頂点0)、必ず頂点0側に属さない頂点をT(=頂点N)とする。頂点を無効化する、ということはその頂点から頂点Tへcapacity=1の辺を張ることで表現できるので、元のグラフに加えて、G個のアクセスできないようにしたい頂点と頂点Tに辺を張れば、S-T間の最小カットが求める解になる。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
#define irep(i, m, n) for (int i = (int)(m); i < (int)(n); ++i)
using ll = long long;
template<class t> using vc=vector<t>;
const ll INF = 1e9+10;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
signed main()
{
cin.tie( 0 ); ios::sync_with_stdio( false );
int n,g,e;
cin>>n>>g>>e;
int s=0, t=n;
auto flow=Dinic<int>(n+1);
rep(i,g){
int p; cin>>p;
flow.add_edge(p,t,1);
}
rep(i,e){
int a,b; cin>>a>>b;
flow.add_edge(a,b,1);
flow.add_edge(b,a,1);
}
cout<<flow.get_max_flow(s,t)<<endl;
}
TopCoder SRM 594 Div I FoxAndGo3
https://community.topcoder.com/stat?c=problem_statement&pm=12808&rd=15706
診断人さんのスライドにも載っている問題。（P.36参照）
問題
N*Nの囲碁の盤面の情報が与えられる（各マスは空白、白、黒のいずれか）。任意の空白マスに好きなだけ黒石を置くことができ、黒で白を囲ったときその白石を取ることができる。黒石を最適に配置した時、最終的な空白マスの最大個数を求めよ。
制約
白石は左右上下に連続して配置されることはない
解説
空白と白は最終的に全て空白にできると仮定し、その仮定した個数から空白にできない個数の最小値を引くことで解を求める方針で解く。
空白、白が置かれているマスそれぞれについて、そのままにするか、変更するか（白を取る、空白に黒を置く）かの2つの操作が行えると考える。
必ずそのままにする頂点をS、必ず変更する頂点をTとし、N*N+2頂点のグラフに適切に辺を張ることで、最小カットが求めたい空白にできない個数の最小値になる。
以下はその辺の張り方。
- 空白マスは、そのままにするならコスト0、変更する（黒石を置く）ならコスト1掛かる、と捉えられる
- 白マスは、そのままにするならコスト1、変更する（白石を取る）ならコスト0掛かる、と捉えられる
更に白マスを変更するための条件を整理すると、ある白マスを変更するなら、その周囲の4マスにある空白は全て変更する必要がある、ということになる。 これは言い換えると、ある白マスを変更して、その周辺の空白マスのいずれかが変更されていない、ということは許されないということになる。 許されないという状態はcapacity=無限の辺を張れば表現できるので、白マスの頂点から周辺の空白マスへの頂点へその辺を張る。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
#define irep(i, m, n) for (int i = (int)(m); i < (int)(n); ++i)
using ll = long long;
template<class t> using vc=vector<t>;
const ll INF = 1e9+10;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
const int dy[] = {0, 1, 0, -1};
const int dx[] = {1, 0, -1, 0};
inline bool inside(int y, int x, int H, int W) {
return (y >= 0 && x >= 0 && y < H && x < W);
}
class FoxAndGo3{
public:
int h,w;
int node(int i, int j){
return i*w+j;
}
int maxEmptyCells(vector<string> mat){
h=mat.size(), w=mat[0].size();
int s=h*w, t=h*w+1;
auto flow=Dinic<int>(h*w+2);
int ans=0;
rep(i,h){
rep(j,w){
if(mat[i][j]=='x') continue;
ans++;
int now=node(i,j);
if(mat[i][j]=='.'){
flow.add_edge(s,now,0);
flow.add_edge(now,t,1);
}else{
flow.add_edge(s,now,1);
flow.add_edge(now,t,0);
rep(k,4){
int ni=i+dy[k], nj=j+dx[k];
if(inside(ni,nj,h,w) and mat[ni][nj]=='.'){
flow.add_edge(now,node(ni,nj),INF);
}
}
}
}
}
return ans-flow.get_max_flow(s,t);
}
};
// test
// signed main()
// {
// cin.tie( 0 ); ios::sync_with_stdio( false );
// auto f=FoxAndGo3();
// vc<string> mat(3);
// mat[0]="xox";
// mat[1]="o.o";
// mat[2]="xox";
// cout<<f.maxEmptyCells(mat)<<endl; // 4
// }
RUPC 2019 Day 2 H - Ghost
問題
N人の人がいて、最初それぞれ左右どちらかを向いている。各iについて初期状態と反対の方向を向くときコストA_iがかかる。更にM個の条件が与えられて、(x,y)が向き合っているならコストB_jがかかる。各人について最適に方向を決めた時のコストの最小値を求めよ。
制約
解説
最終的に必ず左を向いている頂点をS、必ず右を向いている頂点をTとして、N+2頂点のグラフに適切に辺を張ることで、S-T間の最小カットが求める解になる。
辺の張り方は以下の様になる。
まず各人iとS,T間にそれぞれ、初期方向と同じ方向ならコスト0、反対方向ならコストA_iの辺を張る。
また(x,y)が向き合っているということは、x<y であるとして、人xが右を向いていて人yが左を向いている状態である。よって、x→yへコストB_jのcapacityの辺を張ればよい。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
template<class t> using vc=vector<t>;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
signed main()
{
cin.tie( 0 ); ios::sync_with_stdio( false );
int n,m; cin>>n>>m;
string u; cin>>u;
auto flow=Dinic<int>(n+2);
int s=n, t=n+1;
rep(i,n) {
int a; cin>>a;
if(u[i]=='L'){
flow.add_edge(s,i,0);
flow.add_edge(i,t,a);
}else{
flow.add_edge(s,i,a);
flow.add_edge(i,t,0);
}
}
rep(i,m){
int x,y,b; cin>>x>>y>>b; x--,y--;
if(x>y) swap(x,y);
flow.add_edge(x,y,b);
}
cout<<flow.get_max_flow(s,t)<<endl;
}
Educational Codeforces Round 55 G. Petya and Graph
問題
N頂点M辺のグラフが与えられる。各頂点、辺にはそれぞれ重みA、Wがあり、部分グラフの重みを（辺の重みの和）-（頂点の重みの和）とするとき、重みの最大値を求めよ。なお、部分グラフは連結である必要はない。
制約
解説
全ての辺の重みのみを得られると仮定し、仮定した合計から、（得られない辺の重み+含まなければならない頂点の重み）を引くことで解を求める。
各頂点、および各辺を部分グラフに含むか、含まないかのいずれかに分けたときの最小カットが求めたい重みの総和になるようにグラフを構築する。
- 辺iが部分グラフに含まれないなら、コストA_i掛かると捉えられる
- 頂点jが部分グラフに含まれるなら、コストW_j掛かると捉えられる
また、ある辺が部分グラフに含まれて、その辺の両端の頂点が含まれない、という事は許されないので、（頂点Sが必ず「部分グラフに含まれる」、頂点Tが必ず「部分グラフに含まれない」とすれば）その2頂点から辺へcapacity=無限の辺を張るようにする。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
template<class t> using vc=vector<t>;
using ll = long long;
const ll LINF = 1e18;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
signed main()
{
cin.tie( 0 ); ios::sync_with_stdio( false );
int n,m; cin>>n>>m;
auto flow=Dinic<ll>(n+m+2);
int s=n+m, t=n+m+1;
rep(i,n){
ll a; cin>>a;
flow.add_edge(s,i,a);
flow.add_edge(i,t,0);
}
ll ans=0;
rep(i,m){
int u,v; ll c; cin>>u>>v>>c;
u--,v--;
ans+=c;
int e=i+n;
flow.add_edge(s,e,0);
flow.add_edge(e,t,c);
flow.add_edge(u,e,LINF);
flow.add_edge(v,e,LINF);
}
ans-=flow.get_max_flow(s,t);
cout<<ans<<endl;
}
RUCPC 2019 Day 1 G イルミネーション
問題
h*wのグリッドが与えられる。各マスには、B[i][j] の得点が配置されており、また (i+j) が偶数のマスにはスイッチがあり、1つのスイッチをONにすると B[i][j]+B[i][j+1]+B[i+1][j]+B[i+1][j+1] の得点が得られるがコストがW掛かる。スイッチのON、OFFを適切に設定した時、得られる得点の最大値を求めよ。
制約
解説
あるマスが2つのスイッチの対象になるなら得点を二重にカウントしないように注意する必要があるが、一旦これを考えないようにすれば、あるスイッチをONにしたとき得られる利得は、max(0, B[i][j]+B[i][j+1]+B[i+1][j]+B[i+1][j+1]-W) となる。全てのスイッチについて、この利得が得られると仮定して、実際には得られない利得の最小値を求めこれを仮定した値から減算することで解を求める。
この問題では、各スイッチについて、ONにするかOFFにするかのS-Tグラフを構築する。あるスイッチをOFFにするなら辺の capacity は上に示した通り max(0, B[i][j]+B[i][j+1]+B[i+1][j]+B[i+1][j+1]-W) となる。また、ある2つのスイッチが共有するマスを持つなら、その2つのスイッチを「両方ONにしたとき共有するマスのスコアだけ罰金」という辺を張りたい。「両方がSに属するとき罰金」はそのままでは扱えないので、偶数行目にあるスイッチと奇数行目にあるスイッチで、頂点S/TのON/OFFを入れ替える。こうすることで、偶数行目にある頂点から奇数行目にある頂点に辺を張ることで、「両方ONのとき共有するマスのスコアだけ罰金」という辺を張ることができる。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
template<class t> using vc=vector<t>;
using ll = long long;
const ll LINF = 1e18;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
int h,w;
int node(int i, int j){
return i*w+j;
}
signed main()
{
cin.tie( 0 ); ios::sync_with_stdio( false );
ll W; cin>>h>>w>>W;
vc<vc<ll>> b(h, vc<ll>(w));
rep(i,h) rep(j,w) cin>>b[i][j];
auto flow=Dinic<ll>(h*w+2);
int s=h*w, t=h*w+1;
ll ans=0;
rep(i,h-1){
rep(j,w-1){
if((i+j)%2==0){
ll sum=0;
sum+=b[i][j];
sum+=b[i][j+1];
sum+=b[i+1][j];
sum+=b[i+1][j+1];
sum-=W;
if(sum>0){
ans+=sum;
}
int now=node(i,j);
if(i%2==0){
flow.add_edge(s,now,max(0LL,sum));
flow.add_edge(now,t,0LL);
if(i+1<h-1 and j+1<w-1){
int nxt=node(i+1,j+1);
flow.add_edge(now,nxt,b[i+1][j+1]);
}
if(i+1<h-1 and j-1>=0){
int nxt=node(i+1,j-1);
flow.add_edge(now,nxt,b[i+1][j]);
}
if(i-1>=0 and j+1<w-1){
int nxt=node(i-1,j+1);
flow.add_edge(now,nxt,b[i][j+1]);
}
if(i-1>=0 and j-1>=0){
int nxt=node(i-1,j-1);
flow.add_edge(now,nxt,b[i][j]);
}
}else{
flow.add_edge(s,now,0LL);
flow.add_edge(now,t,max(0LL,sum));
}
}
}
}
auto res=flow.get_max_flow(s,t);
ans-=res;
cout<<ans<<endl;
}
Google Code Jam 2008 World Finals - The Year of Code Jam
codingcompetitions.withgoogle.com
問題
M*Nのグリッドが与えられ、各マスには、#,.,?のいずれかが記載されている。#マスには初期値として4のスコアが与えられ、上下左右に隣接するマスにある#の数だけスコアが1減算される。各?マスを#、. をいずれかのマスとして適切に割り振ったとき、スコアの合計の最大値を求めよ。
制約
解説
まずスコアが加算される可能性があるマスは?か#のいずれかである。 #の場合は、4-(周囲の#の個数)-(周囲の?のうち#にするものの個数)がそのマスが得るスコアとなり、?マスの場合も自身を#マスとするなら同様となる。 そしていずれの場合でも周囲の?マスは . である場合の方がスコアは大きくなる。 この問題も上述のイルミネーションの解法と同様の方針で、全てのマスについて最大のスコアを得られると仮定して、実際には得られないスコアの最小値を求め減算することで解が得られる。
各?マスについて、#にするか、. にするかで以下の様に変わる。
- #にするなら仮定した値は得られるが、周囲の#の個数*2のスコアが減算される
- .にするなら仮定した値は得られない
また、隣り合う?同士について、「その両方が#になるなら2点減算される」という条件が存在することになる。
「ある2つの要素がどちらともSに属するなら罰金」という条件は扱えないので、これもイルミネーションの時と同様の工夫として、マスの偶奇に応じて、#と.の方向を入れ替えるようにする。これで偶数マスから奇数マスへ辺を張ることで条件を扱えるようになる。
解答コード
#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for (int i = (int)(0); i < (int)(n); ++i)
#define reps(i, n) for (int i = (int)(1); i <= (int)(n); ++i)
template<class t> using vc=vector<t>;
using ll = long long;
const ll LINF = 1e18;
template <class T>
class Dinic {
private:
const int INF = 1e9;
vector<int> level, itr;
struct Edge {
int to, rev; T cap;
Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
};
vector<vector<Edge>> G;
Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
void bfs(int s) {
level.assign(G.size(), -1);
level[s] = 0;
queue<int> q;
q.push(s);
while (!q.empty()) {
auto v = q.front();
q.pop();
for (auto &e : G[v]) {
if (level[e.to] < 0 and e.cap > 0) {
level[e.to] = level[v] + 1;
q.push(e.to);
}
}
}
}
T dfs(int v, int t, T flow) {
if (v == t) return flow;
for (int &i = itr[v]; i < G[v].size(); i++) {
auto &edge = G[v][i];
if (level[v] < level[edge.to] and edge.cap > 0) {
auto f = dfs(edge.to, t, min(flow, edge.cap));
if (f > 0) {
edge.cap -= f;
get_rev(edge).cap += f;
return f;
}
}
}
return 0;
}
public:
Dinic(int n) { G.resize(n); }
void add_edge(int from, int to, T cap) {
G[from].push_back(Edge(to, (int)G[to].size(), cap));
G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
}
T get_max_flow(int s, int t) {
int n = G.size();
T res = 0;
while (true) {
itr.assign(n, 0);
bfs(s);
if (level[t] < 0) break;
while (true) {
T flow = dfs(s, t, INF);
if (flow > 0) res += flow;
else break;
}
}
return res;
}
};
int h,w;
int node(int i, int j){
return i*w+j;
}
const int dy[] = {0, 1, 0, -1};
const int dx[] = {1, 0, -1, 0};
inline bool inside(int y, int x, int H, int W) {
return (y >= 0 && x >= 0 && y < H && x < W);
}
void solve(int idx){
cin>>h>>w;
vc<vc<int>> G(h,vc<int>(w));
rep(y,h) rep(x,w) {
char c; cin>>c;
if(c=='?') G[y][x]=1;
if(c=='#') G[y][x]=2;
}
ll ans=0;
int s=h*w, t=h*w+1;
auto flow=Dinic<ll>(h*w+2);
rep(y,h){
rep(x,w){
int now=node(y,x);
int score=4;
if(G[y][x]==1){
rep(i,4) {
int ny=y+dy[i], nx=x+dx[i];
if(inside(ny,nx,h,w)){
if(G[ny][nx]==1){
if((y+x)%2==0){
int nxt=node(ny,nx);
flow.add_edge(now,nxt,2);
}
}
if(G[ny][nx]==2){
score-=2;
}
}
}
score=max(0,score);
ans+=score;
if((y+x)%2==0){
flow.add_edge(s,now,score);
flow.add_edge(now,t,0);
}else{
flow.add_edge(s,now,0);
flow.add_edge(now,t,score);
}
}else if(G[y][x]==2){
rep(i,4){
int ny=y+dy[i], nx=x+dx[i];
if(inside(ny,nx,h,w) and G[ny][nx]==2) score--;
}
ans+=score;
}
}
}
ll cost=flow.get_max_flow(s,t);
ans-=cost;
cout<<"Case #"<<idx<<": "<<ans<<endl;
}
signed main()
{
cin.tie( 0 ); ios::sync_with_stdio( false );
int t; cin>>t;
reps(i,t){
solve(i);
}
}
参考
www.slideshare.net
ARC 117 C - Tricolor Pyramid の mod3 まわりについての個人的な補足
問題
公式解説
解説にある通り、青白赤を0,1,2に対応させることで各ブロックは二項係数を使って最下段のブロックのから（二項係数の計算量は無視して）O(N)で求めることができるという問題。
この問題には更に 二項係数の計算において mod 3を求める という副題があり、こちらについて学びになった点をメモとして残す。
乗法逆元についておさらい
法Mが与えられたとき、Xの乗法逆元が求められるための条件は、MとXが互いに素であることである。（この条件を、「Mが素数であればよい」と勘違いしないように注意する）
競プロでよくあるのは、法として素数である 109+7 や 998244353 が与えられて、計算に扱う値がそれ未満の自然数なので、法と扱う数が必ず互いに素となる、というパターンでこの場合は逆元求めて計算を行うことができる。
一方で、この問題では法が3で、扱う数が3と互いに素ではない可能性があるため逆元を求めて計算を行おうとすると計算結果が壊れる。
例として を考えると、以下の様に分母に3が出てくるために逆元を求めることができない。
mod 3の計算方法の解説に対する補足
まず前提として二項係数は以下の様に式変形ができる。
この時、右辺も当然整数となる。
が3の倍数になるための条件
解説には という関数を定義して、
の場合、明らかに3の倍数であると言っている。
この条件によって何故3の倍数と言えるかを補足する。
まず、 が整数であることから、
とすると、
であることがわかる。（
となるなら明らかに整数にならない）
ここで、ある数が3の倍数であるなら という形で表すことができるはずで、
をこのように表すことができるための条件が
つまり
である、と解説では言っているのである。（分子側の3の個数が多くないと全体で3の倍数に成り得ないということ）
また上記の説明は、この条件を満たすときに限り は3の倍数になるということも示唆している。
が3の倍数以外の場合について
上記から3の場合とならないなら となることが分かった。
これはつまり、 における分母と分子の3を相殺できること意味し、相殺後の分母は3と互いに素になっているため、mod 3における逆元を求めることができる。
よって の場合は、
から3で割れるだけ割ってできた整数の mod 3 を求め（解説ではこれを
という関数で定義している）、あとは乗法逆元を使って
を求めればよい。（解説ではここを場合分けして求めている）
提出コード
おまけ
問題の副題として、法が 109+7 や 998244353 ではなく、逆元を求めたい数と互いに素ではない可能性があるという問題
問題ネタバレ
行列累乗が本題の問題。 解説: https://pione.hatenablog.com/#ARC-050-C---LCM-111
行列累乗についてのメモと出題例
最近DPの漸化式を行列式に変換するセンスを磨くために関連問題を解いており、問題や行列累乗全般について雑多なメモを残しておこうと思う。
フィボナッチ数列のN項目を高速に求める
行列累乗の有名な活用例としてよくフィボナッチ数列が挙げられる。
オンラインジャッジで全く同じ問題は見つけられなかったが、蟻本にその記載がある。
フィボナッチ数列の漸化式は以下の様になる。
- dp[0]=dp[1]=1
- dp[i+2]=dp[i+1]+dp[i]
この遷移をナイーブに実装すると第N項を求めるためにO(N)の時間が掛かる。
これを高速に求めるため、上記の漸化式を行列積の形に書き換える。（f_iをフィボナッチ数列の第i項目とすると）
{f_1, f_0} = {{1,1}, {1,0}} * {1,0}
更に右辺に {{1,1}, {1,0}} を掛けると {f_2, f_1}、更に掛けると {f_3, f_2} を求めることができる。
つまり第N項目を求めるためには {{1,1}, {1,0}} のN乗を求めればよい。累乗を高速に求める手法に繰り返し二乗法というものがあるが、これが行列にも応用することができる。この高速に行列の累乗を求める手法を競プロ界隈では「行列累乗」と呼称しているらしい。
以下は、行列累乗によりフィボナッチ数列を求めるための実装
フィボナッチ数列
#define mat vector<vector<ll>> /// 行列積 mat mat_mul(mat &a, mat &b) { mat res(a.size(), vector<ll>(b[0].size())); for (int i = 0; i < a.size(); i++) { for (int j = 0; j < b[0].size(); j++) { for (int k = 0; k < b.size(); k++) { (res[i][j] += a[i][k] * b[k][j]) %= MOD; } } } return res; } /// 行列累乗 mat mat_pow(mat a, long long n) { mat res(a.size(), vector<ll>(a.size())); // 単位行列で初期化 for (int i = 0; i < a.size(); i++) res[i][i] = 1; // 繰り返し二乗法 while (n > 0) { if (n & 1) res = mat_mul(a, res); a = mat_mul(a, a); n >>= 1; } return res; } signed main() { ll n; cin >> n; // フィボナッチ数列の行列累乗 mat m(2, vc<ll>(2)); m[0][0] = m[0][1] = m[1][0] = 1; m[1][1] = 0; m = mat_pow(m, n); ll ans = m[1][0]; cout << ans << endl; }
例題をいくつか
以下に行列累乗によって解ける問題を載せておく。
下に行くほど難しくなっている。
ABC 009 D - 漸化式
この問題で漸化式を立てると以下の様になる。
A[i] = (C[0]&A[i-1]) ^ (C[1]&A[i-2]) ^ ...
この遷移を観察すると、自身より前の項によって自身の項が定まるためフィボナッチ数列DPに似ている。
また、求めたい第M項のMの値が大きい（最大109）かつ、自身を求めるために参照する自身より前の項の数Kが小さい（最大100）であるということから計算量から行列累乗が解法の候補に挙がる。（行列累乗で解けるなら、Kが行列の列幅行幅に対応するので、その計算量はO(K3logM)と十分に高速である。）
ここで、この問題では演算が加法、乗法ではなく、排他的論理和（XOR）、論理積（AND）を使うが半環を満たすため行列累乗が使える。
累乗する行列は以下の様になる。
※ 以下の画像はは解説スライドのP.102の引用です。
この問題では乗法が論理積に置き換わっているので、ここでの1とは、論理積を行った結果が変化しない値である。Aの要素は32ビット符号なし整数なので、32ビット1が立った値、(1<<32)-1 などを扱えばよい。
解説スライド
Abc009 from AtCoder Inc.
www.slideshare.net
解答コード
ARC 050 C - LCM 111
X桁のレピュニット数をf(X)と表すことにする。
レピュニット数同士の剰余を観察してみると（A>Bとする）
f(A)%f(B) = f(A%B)
という関係があることがわかる。レピュニット数同士の剰余がその桁数の剰余に対応することから、レピュニット数同士のgcdはその桁数のgcdによって求めることができる。つまり、
gcd(f(A), f(B)) = f(gcd(A, B))
という関係が成り立つ。
gcd(A, B)をGとすると求める解は、
f(A)*f(B)/f(G)
で求めることができるので、後は行列累乗でf(A), f(B), f(G)を高速に求めれば良さそう...に思えるが、この問題で与えられる法Mは適当な正整数であり（f(G)と互いに素であるとは限らないため）f(G)の逆元を求められない可能性がある。
ひとまず f(A) を求めることを考える。
f(i+1) = f(i) * 10 + 1
なので、これを行列式で表すと以下の様になる。
{f(2), f(1)} = {{10,1},{0,1}}*{1,1}
よって、f(N)を求めるには、{{10,1},{0,1}}のN乗を行列累乗で求めればよいとわかる。
残る f(B)/f(G) について、これは逆元を使わずに求めることにする。f(B), f(G)はいずれもレピュニット数で f(B)>=f(G) なので、このような関係性の2数の除法について観察してみる。
f(6) / f(2) = 111111 / 11 = 10101
この 10101 という値はレピュニット数の遷移とほぼ同じで求めることができる。つまり、以下の様に表せる。
g(i+1) = g(i) * 10G + 1
これを行列式で表すと以下の様になる。
{g(2), g(1)} = {{10G,1},{0,1}}*{1,1}
B/G = g とすると、f(B)/f(G) は {{10G,1},{0,1}}*{1,1} をg乗することで求めることができる。
解答コード
EDPC R - Walk
この問題もKという値が大きい（最大1018）かつ、行列の列幅になるであろうNの値が小さい（最大50）ので、行列累乗で解けるならO(N3logK)時間になり間に合いそう。（さらにこの問題は数え上げである）
まず素直に漸化式を考えると、以下の様な遷移が考えられる。
dp[len][pos] := 経路長lenで終点がposであるものの通り数
- dp[0][i] = 1 （初期化）
- dp[len+1][to] += dp[len][from]
簡単のため N=2 という小さい場合について考える。有向辺は以下の様に与えられるとする。
0 0 1 0 // 頂点2->頂点1の有向辺
この行列を(i, j)=(頂点i <- 頂点j の有向辺が存在する)というように変換する。
0 1 // 頂点1<-頂点2の有向辺 0 0
この行列を使えば、以下の様に行列式で漸化式を表すことができる。
{dp[1][0], dp[1][1]} = {{0,1}, {0,0}} * {1,1}
経路長2の長さを求めるに右辺に {{0,1}, {0,0}} を掛ければよく、つまり経路長Kの長さを求めるにはこのK乗を求めればよい。
以上を有向グラフ全般に適用させるためには、与えられる行列Aから、(頂点i <- 頂点j の有向辺が存在する)行列を構築してそのK乗を求めればよい。
解答コード
yukicoder No.840 ほむほむほむら
h(ほ) -> m(む) -> r(ら) の並び順の部分文字列の個数がK倍個であるような文字列の個数を数え上げる問題。
この問題も、文字列の長さとして与えられるNが大きく（最大109）、Kが小さい（最大5）。
まず素直にDPの漸化式を考えてみる。
dp[n][i][j][k] := n文字で、hの個数がi個の倍数、hmの個数がj個の倍数、hmrの個数がkの倍数である場合の通り数（よって配列長はN*K3）
- dp[0][0][0][0] = 1 （初期化）
- dp[n][(i+1)%K][j][k] += dp[n-1][i][j][k] // hを追加
- dp[n][i][(j+i)%K][k] += dp[n-1][i][j][k] // mを追加
- dp[n][i][j][(k+j)%K] += dp[n-1][i][j][k] // rを追加
これを行列式で合わらすために [i][j][k] の3次元配列をK進数によって1次元で表すことを考える。
つまり、hがi個、mがj個、rがk個を iK2+jK+k と表すことにする。この値は高々K3(<=53)と小さく、この行幅列幅の行列を累乗しても十分高速である。
あとは、行列mのm(to, from)を状態toへ状態fromからの遷移があるとして初期化し、このN乗を求めればN文字の各状態の通り数を計算できる。このうち、hmrの個数が0倍の個数の総和を求めればそれが解となる。
解答コード
