phyllo’s algorithm note

レッドコーダーへの道のりは遠い。休んでる場合じゃない!

yukicoder No.492 IOI数列

問題

IOI数列は、1,101,10101,1010101,...とa_1=1, a_i=100*a_{i-1}+1であるような数列である。
第N項について、1000000007で割ったあまりと101010101010101010101で割ったあまりを求めよ。

制約

1 <= N <= 10^18

解説

第N項を計算したい。
後者の10101...で割ったあまりは、64bit整数に収まらないが、パターンの繰り返しなので、それで答えればよい。
前者の1000000007で割ったあまりを考える。

特性方程式

「a_n = b*a_n + c」の形をしているので、特性方程式x=bx+cを解いてあげれば等比数列として考えられる。

x=100x+1 <=> x = -1/99

よって、「(a_{n+1}-x)=100*(a_{n}-x)」となるので、一般項は、「a_{n+1} = 100^n * (a_1 - x) + x」となる。
展開すると、a_{n+1}=(100^{n+1} - 1)/99なので、「a_n=(100^n - 1)/99」になる。

この計算を、MODでの累乗(繰り返し二乗法)や逆元(拡張ユークリッド互除法を使う方法)を使ってあげれば高速に求められる。

行列累乗

より一般に、m項間漸化式について、一つ次の項を求めるような行列計算を立てると、第N項を行列累乗で求められる。

v_{n+m} = (a_{n+m}, a_{n+m-1}, ..., a_{n+1})^Tとして、「v_{n+m} = X * v_{n+m-1}」となるような行列Xを作れれば、第n項は「v_{n+m-1} = X^n * v_{m}」のように書ける。

今回の漸化式は定数項があるので、
\begin{pmatrix} a_{n+1} \\ 1 \end{pmatrix} = \begin{pmatrix} 100&1 \\ 0&1 \end{pmatrix} * \begin{pmatrix} a_{n} \\ 1 \end{pmatrix}
から
\begin{pmatrix} a_{n+1} \\ 1 \end{pmatrix} = \begin{pmatrix} 100&1 \\ 0&1 \end{pmatrix}^n * \begin{pmatrix} a_{1} \\ 1 \end{pmatrix}
という式になる。行列の累乗は繰り返し二乗法を使ってO(m^3 log n)または高速にO(m^2 log n)の計算量で求められる。

反省

  • 漸化式は行列累乗を検討

KUPC2016 E. 柵 / Fences

問題

H*Wのグリッドのいくつかのマスにヤギがいる。
ヤギは上下左右の隣接マスに移動できるが、移動先となるマスに柵がある場合は、その方向へは移動できない。
ヤギがグリッドの外に出ないようにするために必要な柵の最小個数を求めよ。

制約

1 <= H <= 100
1 <= W <= 100

解説

http://www.kupc.jp/editorial/2016/E.pdf


最小カットで求めることができる。

マスに柵を置くかどうかを扱うために、各マスをin頂点とout頂点に分解し、ヤギのいないマスの場合はinからoutを容量1でつなぎ、ヤギがいる場合は容量infでつなぐ。
ヤギのいないマスの場合は、さらに、out頂点から上下左右のマスのin頂点へ容量infでつなぐ。

最後に、s-tフローで扱うために、頂点sからグリッドの縁のマスのin頂点へ容量infでつなぎ、ヤギのいるマスのout頂点から頂点tへつないで、s-t最小カットを求めればよい。

AGC010 C. Cleaning

問題

Nノードからなる木が与えらえる。
各ノードにはA[i]個の石が置かれている。
「2つの(異なる)葉ノードを選び、そのパス上のノードすべてから1つずつ取り除く」という操作を繰り返して、すべての石を取り除けるか答えよ。
ただし、パス上に石のないノードがある場合は操作は行えない。

制約

2 <= N <= 10^5
0 <= A[i] <= 10^9

解説

https://atcoder.jp/img/agc010/editorial.pdf

適当にノードを選んで根にして、根付き木として考える。
あるノードvについて、親ノードがw、子ノードがc_1,c_2,...,c_kであるような場合を考える。

vとc_iの間の辺が何回パスとして通過するかをp[c_i]とすると、vとwの間の辺が何回パスとして通過するかp[w]は、ノードvの通過回数A[v]が決まっているので、
p[w] = 2*A[v] - Σp[c_i]
でなければならない。
このとき、条件として、

  • p[w]は0以上A[v]以下の数でなければいけない
  • もし、子ノードが1つの場合は、p[w]=A[v]=p[c_1]でなければならない
  • vの周りの各辺の通過回数p[x]は、辺の通過回数合計/2回 = (p[w]+p[c_1]+...+p[c_k])/2以下でなけれなならない

を満たす必要がある。

これを再帰的に繰り返し(dfs)、根まで求め、

  • 根が葉ノードの場合、根の疑似的な親ノードへの辺の通過回数=A[root]になっている
  • 根が葉ノードでない場合、根の疑似的な親ノードへの辺の通過回数=0になっている

のどちらかになっていれば、すべての石を取り除くことができる。

反省

  • 根付き木で考える
  • 葉ノードから決定的に決まるものを見つける

ARC068 E. Snuke Line

問題

直線上の鉄道で、0~Mまでの番号が付いたM+1個の駅がある。
N種類の名産品が与えられ、名産品iは駅lから駅rまでの区間で売られている。
最初0番の駅にいるとして、d駅ごとに停車する電車の場合、購入可能な名産品の種類数を知りたい。
dが1~Mの場合それぞれについて、種類数を答えよ。

制約

1 <= N <= 3*10^5
1 <= M <= 10^5
1 <= l_i <= r_i <= M

解説

AtCoder Regular Contest 068/ Beginner Contest 053 解説放送 - YouTube
解説放送の方の方法で解いてみた。


よく知られた事実として、以下のループはO(M log M)程度であることが使える。

for(int d=1; d<=M; d++){
  for(int p=1; d*(long long)p<=M; p++){
    //d*pについて、何かする
  }
}


単純な方法だと、「l_i <= d*p <= r_i」であるかをN個の名産品でチェックする方法がある。
しかし、これだとTLEしてしまうし、気を付けないと、d*pもd*(p+1)も範囲内であるような「l_i <= d*(p+1) <= r_i」のケースでダブルカウントしてしまうかもしれない。


範囲l_i~r_iに対して、d*pが1回しかカウントしないようにするうまい方法として、「d*(p-1)は範囲外で、d*pが範囲内の場合カウント」することで、範囲内の停車駅のうち、一番左の駅を見つけることができる。
X=d*(p-1)、Y=d*pとおいて考えると、「X < l_i、かつ、l_i <= Y <= r_i」を満たすような特産品を高速に求められれば良いことになる。


範囲は扱いにくいので、「l~r」という情報を2次元平面の点p(l,r)として考える。
上のXとYの条件はこの2次元平面上でどのような条件になるか考えると、分解して、

  • X < l
  • l <= Y
  • Y <= r

を満たす範囲の点の数を見つけることに相当する。
この3つの条件は、以下の*の部分のような平面上の長方形になっている。

M |-----+-------+
  |     |*******|
r |-----+***p***|
  |     |*******|
Y |-----+---+---|
  |     |   |   |
  +--------------
        X   l   Y

これを高速に処理することを考える。


今回、問い合わせクエリ(X,Y)と特産品の点(l,r)が先に求めておける。
そこで、問い合わせと特産品の点を合わせてl軸方向でソートしておき、r軸方向に対する1次元BITを使って、l軸の小さい方から点を置いたり、問い合わせに答えたりをすれば必要な(X,Y)の答えを求めておける。
あとは、X,Yの個数を使って該当する特産品数を求めて表示すればよい。

反省

(X,Y)の条件を満たす点の数を保存するのに、

map<pair<int,int>,int> memo;

としてしまったが、O(M log M)のループとはいえ、10^6個以上要素がある。
このままだと、memo[make_pair(X,Y)]などとアクセスした場合、木が深いため、時間が無視できない程度かかってしまう。
代わりに、

map<int,int> memo[100005];

などと木の深さが深くならないようにすることで、時間をかけずにアクセスできるようになる。

コード

#include <bits/stdc++.h>
using namespace std;
#define REP(i,a,n) for(int i=(a); i<(int)(n); i++)
#define rep(i,n) REP(i,0,n)
#define FOR(it,c) for(__typeof((c).begin()) it=(c).begin(); it!=(c).end(); ++it)
#define ALLOF(c) (c).begin(), (c).end()
typedef long long ll;
typedef unsigned long long ull;

#define MAX 100005

class BIT{
  static const int MAX_N = 1 << 18;
  int n, bit[MAX_N+1];
public:
  BIT(int n_){
    n = n_;
    for(int i=0; i<n+1; i++) bit[i] = 0;
  }
  int sum(int i){
    int s = 0;
    while(i>0){
      s+=bit[i];
      i-=i&(-i);
    }
    return s;
  }
  void add(int i, int x){
    while(i<=n){
      bit[i]+=x;
      i+=i&(-i);
    }
  }
};
 
struct ST {
  int type;
  int l, r;
};
bool operator<(const ST& a, const ST& b){
  if(a.l == b.l) return a.type < b.type;
  return a.l < b.l;
}
 
map<int,int> query[MAX];
 
int main(){
  ios::sync_with_stdio(false);
 
  int N, M;
  cin >> N >> M;
 
  vector<ST> pts;
  rep(i,N){
    int l, r;
    cin >> l >> r;
    pts.push_back((ST){0,l,r});
  }

  for(int d=2; d<=M; d++){
    for(int p=1; d*(ll)p<=M; p++){
      int X = d*(p-1);
      int Y = d*p;
 
      if(query[X].count(Y-1)==0){
        query[X][Y-1] = 1;
        pts.push_back((ST){1,X,Y-1});
      }
      if(query[X].count(MAX)==0){
        query[X][MAX] = 1;
        pts.push_back((ST){1,X,MAX});
      }
      if(query[Y].count(Y-1)==0){
        query[Y][Y-1] = 1;
        pts.push_back((ST){1,Y,Y-1});
      }
      if(query[Y].count(MAX)==0){
        query[Y][MAX] = 1;
        pts.push_back((ST){1,Y,MAX});
      }
    }
  }

  sort(ALLOF(pts));
 
  BIT bit(1<<18);
 
  rep(i,pts.size()){
    if(pts[i].type == 0){
      bit.add(pts[i].r+1, 1);
    }else{
      query[pts[i].l][pts[i].r] = bit.sum(pts[i].r+1);
    }
  }

  for(int d=1; d<=M; d++){
    if(d==1){
      cout << N << endl;
    }else{
      int ret = 0;
      for(int p=1; d*(ll)p<=M; p++){
        int X = d*(p-1);
        int Y = d*p;
 
        ret += query[Y][MAX] - query[Y][Y-1] - query[X][MAX] + query[X][Y-1];
      }
      cout << ret << endl;
    }
  }
 
  return 0;
}

CodeFestival2016 Final D.Pair Cards

問題

N枚のカードがあり、カードiには整数X_iが書かれている。
今、「書かれた数字が同じ」または「和がMの倍数」になるような2枚のカードをペアにする。
最大で何ペア作れるか?

制約

2 <= N <= 10^5
1 <= M <= 10^5
1 <= X_i <= 10^5

解説

https://cf16-final.contest.atcoder.jp/data/other/cf16-final/editorial.pdf


「和がMの倍数」になるようなペアの処理は、X_i mod Mでグルーピングすると、余りがjのグループとM-jのグループのペアを取ればMの倍数になる。
(ただし、mod M == 0か、mod M == M/2(ただし、M mod 2==0)の場合は、同じグループ内でペアを取る)


こうすれば、グループAとグループBでの最大マッチングを求めて足し合わせればよくなる。


|A|==|B|の場合は、|A|個ペアが作れる。
|A|<|B|の場合、|A|個のペアは作れることがわかるが、|B|-|A|個のうち、「書かれた数字が同じ」場合はペアを取ることができる。
したがって、グループBの方で、|B|-|A|枚以下のカードを使って、書かれた数字が同じになるペアの最大ペア数をmとすると、「m+|A|」ペアが最大マッチングになる。

yukicoder No.374 コイン

問題

2人で以下のゲームを行う。
・半径Aの円の中に、半径Bの円を交互に置いていく
・置いていく円は、半径Aの円からはみ出したり、すでに置いた円とかさなったりしてはいけない
・先に置けなくなった方が負け
・置かれた円は動かせない
それぞれ最適に行動した場合、どちらが勝つか答える。

制約

1 <= A,B <= 2^32

解説

A < Bの場合、先手がすでにおけないので、後手の勝ち。
A >= Bの場合を考える。


2*A/3 < Bの場合は、真ん中に置くと、その後、どこにも置けなくなるので、先手の勝ちということがわかる。
2*A/3 >= Bの場合は、先手がどこに置いても後手に手番が回る。


最適な行動は、まず、真ん中に置いて、後手に手番をわたす。
後手がどこに置いたとしても、「円の中心の点対称の位置に円を置く」ということを繰り返すと、必ず先手はおくことができる。
したがって、先手がこの行動をとった場合、必ず勝てるので、先手の勝ちとなる。

AtCoder Regular Contest 055 C. ABCAC

問題

文字列Sが与えられる。
この文字列は、それぞれ長さ1以上の文字列A,B,Cを「ABCAC」のように連結してできている。
このような分割は何通りあるか?

制約

5<=Sの長さ<=200,000
(部分点:Sの長さ<=2,000)

解法

説明のために、「A B C A' C'」と置く。N=(Sの長さ)とする。


もし、A'が出現する位置がxだとわかったとして、A'の長さが決まれば、C'の長さが決まる。
なので、A=A'となる最大の長さa、C'=Cとなる最大の長さcが高速に求まれば、以下から判定できる。
・a+c<(N-x)ならば、作ることができないので、無し
・x<=(N-x)ならば、Bの部分を作ることができないので、無し
候補となるA'とC'の最大の長さは、1文字以上である必要があるので、それぞれ、min(a,N-x-1)、min(c,N-x-1)
・このとき、組み合わせの数は、「min(a,N-x-1) + min(c,N-x-1) - (N-x) + 1」で求められる
 ・min(a,N-x-1)+min(c,N-x-1)==N-xの場合は、1通り
 ・min(a,N-x-1)+min(c,N-x-1)==N-x+1の場合は、2通り
 ・、、、


ということで、位置xにおける「A=A'となる最大の長さa」と「C'=Cとなる最大の長さc」を高速に求める必要がある。
エディトリアルではいくつかの方法が紹介されているが、SuffixArrayと高さ配列を使っても解けるので、それを試してみる。


高さ配列は、SuffixArrayが構築されている場合、次の接尾辞とのprefixの一致数をO(N)で求められる。
SuffixArrayの中で、インデクスが0の部分が「Sの先頭からの文字列」を表すので、そこから前後に高さ配列の値が降順になるようにすれば、すべての位置に対する「A=A'となる最大の長さ」を求められる。

インデクス 接尾辞 高さ配列 高さ配列の値をインデクス0から降順になるようにしたもの
11 0 0
10 a 0 0
7 abra 1 1
0 abracadabra 4 4
3 acadabra 1 1
5 adabra 1 1
8 bra 0 0
1 bracadabra 3 0
4 cadabra 0 0
6 dabra 0 0
9 ra 0 0
2 racadabra 2 0

インデクスが0の部分から、高さ配列を前後にそれぞれ降順になるようにしたものが表の一番右。
これは、min()を取りながらO(N)で求められ、これは位置xに対する「A=A'となる最大の長さ」の情報になっている。(ので、O(1)で求められる)


「C'=Cとなる最大の長さc」については、文字列Sをreverseしたものについて同様に行えば、位置x-1で終わるようなCの最大の長さがO(1)で求められる。
したがって、
・SuffixArray構築にO(N)~O(N log^2 N)
・位置xを決めるためにO(N)
・aとcを求めるのにO(1)
・位置xでの組み合わせの数を求めるのにO(1)
なので、全体でO(N)~O(N log^2 N)程度で求められる。

解法2

Z algorithmは上のをO(N)で作ることができる。

//z[i]:=文字列sと、位置iから始まる部分文字列との共通接頭辞の長さ(ただしz[0]=0)
std::vector<int> Zalgorithm(const std::string& s){
  int n = s.length();
  std::vector<int> z(n, 0);

  int L = 0, R = 0;
  for(int i=1; i<n; i++){
    if(i>R){
      L = R = i;
      while(R<n && s[R-L] == s[R]) R++;
      z[i] = R-L;
      R--;
    }else{
      int k = i-L;
      if(z[k] < R-i+1){
        z[i] = z[k];
      }else{
        L = i;
        while(R<n && s[R-L] == s[R]) R++;
        z[i] = R-L;
        R--;
      }
    }
  }
  return z;
}

テンプレ

コードも残したいと思います。


使っているテンプレは以下です。

#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <iostream>
#include <queue>
#include <list>
#include <stack>
#include <map>
#include <numeric>
#include <set>
#include <sstream>
#include <string>
#include <vector>
using namespace std;
#define REP(i,a,n) for(int i=(a); i<(int)(n); i++)
#define rep(i,n) REP(i,0,n)
#define FOR(it,c) for(__typeof((c).begin()) it=(c).begin(); it!=(c).end(); ++it)
#define ALLOF(c) (c).begin(), (c).end()
typedef long long ll;