読者です 読者をやめる 読者になる 読者になる

撹乱順列(完全順列)の数について理解する

撹乱順列とは

{1,2,3,...n}の順列で,どの要素も元の位置にないような順列を撹乱順列(またの名を完全順列)といいます.例えばn=3の撹乱順列として{2,3,1},{3,2,1},{3,1,2}とかが挙げられます.{2,1,3}とかは撹乱順列ではありません.3番目の要素が元の位置にあります.

つまり,{1,2,...n} の順列 P = {p1,p2,...,pn}について考えるとき,全ての i についてpi≠iが成り立つような順列を撹乱順列と呼びます.

そしてその撹乱順列の個数はモンモール数という名前がついており,以下の漸化式で求められます.

個数の漸化式

  • a1 = 1
  • a2 = 1
  • an = (n-1)(an-1+an-2) (n ≧ 3)

なぜこのような漸化式になるのかについて考えてみましょう.説明のために遷移とかいう言葉を使っていますが再帰をイメージしてください.

まず,長さnの撹乱順列 {1,2,..,n}のうち,1の位置をiで確定させることを考えます(n-1通りあります).別に確定させる数字のはnとかでも何でもいいですが,とりあえず1に着目してその位置を考えます. f:id:kyuridenamida:20160427011957p:plain 確定としては,こんな感じですね.いま,位置 i が埋まってしまったのですが,数字 i はフリーです.これはちょっと厄介な状態です. 撹乱順列の漸化式は,ある使われてない数字の元の位置が既に埋まっているような状態に対する式ではありませんでした.

iをどこに置くか考えてみましょう.

数字 i の位置を 1 に行かせるケース

数字 i の位置を 1 に行かせるケースを考えます. f:id:kyuridenamida:20160427012252p:plain

この場合,なんか明らかに残ったやつを再ラベリングすれば{1,2,...n-2}の撹乱順列を求めることと同等になるので,an-2に遷移したら良さそうな感じしますね.このケースはとても簡単でした.

数字 i の位置を 1以外 に行かせるケース

さて,もうひとつが厄介です.以下のような1じゃないところに遷移するケースです. f:id:kyuridenamida:20160427021733p:plain

これ頭の良い人がグッと睨むと,長さn-1の撹乱順列と状態が等価らしいんですが,「は?」状態です.いやいや,下の図の2つの状態が等価てお前. f:id:kyuridenamida:20160427014752p:plain

何言ってるのって気分になりました(図で本来 1 2 3 … i-1 i+1 … nとなるべきところが1 2 3 n とはしょられてるのはご愛嬌)

いやだってさ,上の状態,数字1使われてるやん,しかも数字iの元の位置埋まってるのに数字iフリーじゃん.下の綺麗な状態と等価だなんてとても思えない.

ただ,よく見てみると「数字iの位置は1以外」という条件があります.これがキモっぽいです.

よく考えてみると,上の状態における数字iを数字1に再ラベルしても問題なさそうな感じがしませんか? だって,数字iを1にしたところで,「iの位置は1以外」という性質が「1の位置は1以外」と言い変わるだけで,その性質というのは常に撹乱順列だから元々成り立っているはずです.

この変換によって,上の状態が下の状態と等価な状態になり,長さn-1の撹乱順列の状態であることがわかりました.

まとめ

まとめると,iを選ぶ候補として(n-1)通りあって,あとはiを選ぶとそれぞれのiについてan-1+an-2通りだけ組み合わせの数があることが分かるので,冒頭で述べた漸化式になります.

参考にしたサイト

ちなみに,ほとんどこちら様の受け売りです.こちら様の素晴らしい親切な説明ですら少し悩んだので,悩んだ部分について皆さんと共有しようと思いました.

for-e-study.com

JAG Contest 2016 Domestic C - みさわさんの根付き木

jag2016-domestic.contest.atcoder.jp

誤って書いた記事を消してしまったので再投稿です.

問題文

省略

解法(反転して表示)

LL(1)による再帰構文解析を知っていると特に思考をせず木構造に直すことができる. 木に落としたらマージをするが,これは再帰でやるのが実装が楽でおすすめ.平衡二分木をmerge-splitベースで実装したことがあるとイメージがつきやすいかも. 余談ですが,構文解析,「次に来る文字が△◎×だったら~」みたいな処理を結構アドホックだと思っていたら,first集合やfollow集合とかいう概念があったりするので,そういうのを学んでみて,強く意識すると迷ったり漏れが発生することが少なくなるかも.

ソース

#include <bits/stdc++.h>
using namespace std;
 
 
string S;
int p;
struct T{
    int v;
    T *l,*r; 
    T(){
        l = r = NULL;
    }
}; 
 
T *merge(T *a,T *b){
    if( a == NULL || b == NULL ) return NULL;
    T *r = new T();
    r->v = a->v + b->v;
    r->l = merge(a->l,b->l);
    r->r = merge(a->r,b->r);
    return r;
}
void view(T *x){
    
    if( !x ) return;
    cout << "(";
    view(x->l);
    cout << ")";
    cout << "[";
    cout << x->v;
    cout << "]";
    cout << "(";
    view(x->r);
    cout << ")";
    
}
 
T* expr(){
    if( S[p] != '(' ) return NULL;
    T *res = new T();
    //←の子
    p++;
    res->l = expr();
    p++;
    
    p++;
    int x = 0;
    while( S[p] >= '0' && S[p] <= '9' ) x = x * 10 + S[p++] - '0';
    res->v = x;
    p++;
    //→の子
    p++;
    res->r = expr();
    p++;
    return res;
}
int main(){
    string s,t;
    cin >> s >> t;
    S = s;
    p = 0;
    auto one = expr();
    S = t;
    p = 0;
    auto two = expr();
    view(merge(one,two));
    cout << endl;
}

JAG Contest 2016 Domestic F - 土地相続

jag2016-domestic.contest.atcoder.jp

本番バグらせて解けなかった悲しいね.

問題文

省略

解法(反転して表示)

基本的な方針はDPと場合分け.累積和とかである区間の値の総和がO(1)で求めれるように前処理しとく. マスは全て非負なので,空き空間は特に作って得することがないことにまあ気づく.

なので,N=1のときは,与えられた長方形全ての総和が答え.

N=2のときは,垂直もしくは水平な直線でぶった切ってみて,その2つの総和の最小値が答え.

N=3のとき,構成される長方形群はどう配置しても,必ず垂直もしくは水平に分割できる場所があるので,dp[着目長方形の左上][左下][あと何個の長方形を取り出さなければいけないかk]という大雑把な見積もりではO(k(WH)2)な状態数を扱うDPで解ける.遷移は「垂直・水平の分割位置」と「分かれた上下(or左右)の長方形にkをどう割り振るか」を全て試す.

実際の状態数は4乗にはならない.垂直もしくは水平にぶった切るごとに端としてとりうる値の自由度が増えていくのだけれど,一回分割するごとに自由度が1つ増えるだけなので,N<=3なら自由度が高々2なので,大体二乗のdpになっているはず.

問題は,N=4のケースだが,これは,水平or垂直にぶった切れる場所が1つでもあれば,N=3の時のDPと同様にすれば,状態数O(k(WWH+WHH)),遷移O(W+H)のdpで解ける. しかし,ぶった切れないときがある.どういうケースかというと,以下のようなケースである.これさえ対策すれば良い. 僕は「下限値の二分探索」+「左上の長方形の形を決め打ちすると連鎖的に他の3つの長方形の形が貪欲で定まる」みたいな感じでこのケースの最適値を求めた.人によっては二分探索使ってない人もいるらしい. しかも,時計回りのケースを回転・反転したものに対して適用することで反時計回りの法の実装の手を抜いた. f:id:kyuridenamida:20160424215047p:plain

まとめると,N<=3ならDPするだけ,N>=4なら上記の特殊ケースを試してみて大きいほうを返す. 実際の実装ではN!=4でなくともN=4と仮定して特殊ケースを試したり,W=H=200と拡張している.

ソース

#include <bits/stdc++.h>
using namespace std;
 
 
map< array<int,4>,int> memo[5];
 
int A[210][210];
int B[210][210];
 
int sum[210][210];
int H,W,N;
 
 
inline int S(int b,int a,int d,int c){
    return sum[c][d] + sum[a-1][b-1] - sum[a-1][d] - sum[c][b-1];
}
int f(const int x1,const int x2,const int y1,const int y2,const int n){
    if( n == 1 ){
        return S(x1,y1,x2,y2);
    }
    
    if( memo[n].count({x1,x2,y1,y2}) ) return memo[n][{x1,x2,y1,y2}];
    int ans = 0;
    for(int i = x1+1 ; i <= x2 ; i++){
        int sub = 0;
        
        for(int j = 1 ; j < n ; j++){
            sub = max(sub,min(f(x1,i-1,y1,y2,n-j),f(i,x2,y1,y2,j)));
        }
        ans = max(sub,ans);
    }
    for(int i = y1+1 ; i <= y2 ; i++){
        int sub = 0;
        for(int j = 1 ; j < n ; j++){
            sub = max(sub,min(f(x1,x2,y1,i-1,n-j),f(x1,x2,i,y2,j)));
        }
        ans = max(sub,ans);
    }
    return memo[n][{x1,x2,y1,y2}] = ans;
}
 
bool check(int lower){
    if( lower == 0 ) return true;
    
    for(int h1 = 1 ; h1 < H ; h1++){
        for(int w1 = 1 ; w1 < W ; w1++){
            if( S(1,1,w1,h1) < lower ) continue;
            int h2 = h1;
            while( h2 < H and S(w1+1,1,W,h2) < lower ) h2++;
            if( S(w1+1,1,W,h2) < lower ) continue;
            int w3 = w1;
            while( w3 > 1 and S(w3,h2+1,W,H) < lower ) w3--;
            if( S(w3,h2+1,W,H) < lower ) continue;
            if( S(1,h1+1,w3-1,H) < lower ) continue;
            return true;
        }
    }
    return false;
    
}

 
int special(){
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            sum[i][j] = A[i][j];
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            sum[i][j] += sum[i-1][j];
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            sum[i][j] += sum[i][j-1];

    int l = 0 , r = 1e9;
    while( l != r ){
        int m = (l+r+1) / 2;
        if( check(m) ){
            l = m;
        }else{
            r = m-1;
        }
    }
    return l;
}
int main(){
    cin >> H >> W >> N;
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            cin >> A[i][j];
    H = W = 200;
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            sum[i][j] = A[i][j];
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            sum[i][j] += sum[i-1][j];
    for(int i = 1 ; i <= H ; i++)
        for(int j = 1 ; j <= W ; j++)
            sum[i][j] += sum[i][j-1];
            
    int ans = f(1,W,1,H,N);

        // N = 4 の入り組んだケース対策.盤面を回転・反転して同じ手法を試してる.無駄があるかも.
    for(int __ = 0 ; __ < 2 ; __++){
        for(int _ = 0 ; _ < 4 ; _++){
            for(int i = 1 ; i <= H ; i++){
                for(int j = 1 ; j <= W ; j++){
                    B[i][j] = A[j][i];
                }
            }
            for(int i = 1 ; i <= H ; i++){
                for(int j = 1 ; j <= W ; j++){
                    A[i][j] = B[i][j];
                }
            }
            ans = max(special(),ans);
            // cout << hoge() << endl;
        }
        for(int i = 1 ; i <= H ; i++){
            for(int j = 1 ; j <= W/2 ; j++){
                swap(A[i][j],A[i][W-j+1]);
            }
        }
    }
    cout << ans << endl;
 
}

JAG Contest 2016 Domestic E - 選挙活動

ICPC系

jag2016-domestic.contest.atcoder.jp

問題文

省略

解法(反転して表示)

オーソドックスな幾何. 列挙して意味のある候補点を列挙する.この候補点に意味のない点が含まれていてもいいので雑に. 候補点を列挙するために,有権者から各多角形の頂点への直線を伸ばして,それを列挙してみる(見える範囲が変化する角度に対する線分を考えるとそれを列挙するのが合理的そうなことが分かる). で,得られた直線全てについて,2つの直線の交点を求めて,それを候補点とする.

候補点が求まれば,各点について,それを見ることのできる有権者の数を数えれば良い. 条件の「多角形の頂点または有権者のいる場所の座標のうち3点が同一直線状に存在することはない.」がありがたい制約で,これのおかげで幾何パートがかなり楽になっていると言っても良い. 多分実装の軽い方針で多角形と線分の交差判定をする必要があるんですが,その制約のおかげで以下のような図のケースを考えなくて良いため,楽になる. f:id:kyuridenamida:20160424212830p:plain

計算量はよく見積もってないけど多項式時間で終わるしそんなに大きくない.

ソース

#include <bits/stdc++.h>
using namespace std;
 
const double EPS = 1e-8;
const double INF = 1e12;
typedef complex<double> P;
typedef vector<P> G;
namespace std {
  bool operator < (const P& a, const P& b) {
    return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b);
  }
}
 
double cross(const P& a, const P& b) {
  return imag(conj(a)*b);
}
 
double dot(const P& a, const P& b) {
  return real(conj(a)*b);
}
struct L : public vector<P> {
  L(const P &a, const P &b) {
    push_back(a); push_back(b);
  }
};
 
int ccw(P a, P b, P c) {
  b -= a; c -= a;
  if (cross(b, c) > 0)   return +1;     
  if (cross(b, c) < 0)   return -1;      
  if (dot(b, c) < 0)     return +2;       
  if (norm(b) < norm(c)) return -2;       
  return 0;
}
 
P crosspoint(const L &l, const L &m) {
  double A = cross(l[1] - l[0], m[1] - m[0]);
  double B = cross(l[1] - l[0], l[1] - m[0]);
  if (abs(A) < EPS && abs(B) < EPS) return m[0]; 
  if (abs(A) < EPS) assert(false); 
  return m[0] + B / A * (m[1] - m[0]);
}
 
// 以下の条件式で,0未満としているのは(ccwの返り値は-EPSなのは全く意味ない),線分が端点で接しているときは当たっていないことにしたいため.
bool intersectSS(const L &s, const L &t) {
  return ccw(s[0],s[1],t[0])*ccw(s[0],s[1],t[1]) < -EPS &&
         ccw(t[0],t[1],s[0])*ccw(t[0],t[1],s[1]) < -EPS;
}
bool intersectLL(const L &l, const L &m) {
  return abs(cross(l[1]-l[0], m[1]-m[0])) > EPS || 
         abs(cross(l[1]-l[0], m[0]-l[0])) < EPS;   
}
 
 
P in(){
    double x,y;
    cin >> x >> y;
    return P(x,y);
}
int main(){
    int N,M;
    cin >> N >> M;
    vector<L> ls;
    vector<G> gs; 
    vector<P> check;
    vector<P> look;
    for(int i = 0 ; i < N ; i++){
        int L;
        cin >> L;
        G g;
        for(int j = 0 ; j < L ; j++){
            g.push_back(in());
            check.push_back(g.back());
            look.push_back(g.back());
        }
        gs.push_back(g);
    }
    vector<P> people;
    for(int i = 0 ; i < M ; i++){
        people.push_back(in());
    }
    for(int i = 0 ; i < M ; i++){
        for(int j = 0 ; j < check.size() ; j++){
            ls.push_back(L(check[j],people[i]));
        }
    }
    
    for(int i = 0 ; i < ls.size() ; i++){
        for(int j = i + 1 ; j < ls.size() ; j++){
            if( intersectLL(ls[i],ls[j]) ){
                look.push_back(crosspoint(ls[i],ls[j]));
            }
        }
    }
    int ans = 0;
    for( auto lp : look ){
        int sub = 0;
        for(int i = 0 ; i < M ; i++){
            int f = 1;
            for(int j = 0 ; j < gs.size() ; j++){
                for(int k = 0 ; k < gs[j].size() ; k++){
                    L l(gs[j][k],gs[j][(k+1)%gs[j].size()]);
                    if( intersectSS(l,L(people[i],lp)) ) f = 0;
                }       
            }
            if( f ) sub++;
        }
        ans = max(ans,sub);
    }
    cout << ans << endl;
}

JAG Contest 2016 Domestic D - インビジブル

ICPC系

jag2016-domestic.contest.atcoder.jp

問題文

省略

解法(反転して表示)

よくあるゲーム木のmin-max探索,後手は差を最小化,先手は差を最大化するように実装する. 状態数がどれくらいになるのかを見積もるのが本質. 実際山札の残り枚数と,どの区間(※)がスタックに残ってるかみたいなのを考えると,状態として O(NM*(N+M)2)程度しか無い.実際はパスの回数とかどっちのターンかみたいな状態も持つんだけど,それはとても小さい定数倍だからまあ問題ない. 実装を複雑化させないために構造体とmapを用いて実装したので重いけど通る.

(※) 手札の引き方はパスを無視すれば一意で,a_1b_1a_2b_2 ,...というふうな列になるので,この列の区間を考える.

「スタックが空になってから」パスを2回が終了条件な注意しないとWAする.

ソース

#include <bits/stdc++.h>
using namespace std;
 
 
 
struct State{
    vector<int> p1,p2;
    vector< pair<int,int> > stack;
    int turn;
    int passed;
};
bool operator < (const State &a,const State &b){
    if( a.passed != b.passed ) return a.passed < b.passed;
    if( a.turn != b.turn ) return a.turn < b.turn;
    
    if( a.p1 != b.p1 ) return a.p1 < b.p1;
    if( a.p2 != b.p2 ) return a.p2 < b.p2;
    return a.stack < b.stack;
}
 
map<State,int> memo;
int dfs(const State s){
    if( s.passed == 2 ) return 0;
    if( memo.count(s) ) return memo[s];
    int ans;
    if( s.turn == 0 ){
        ans = -1e9;
        if( s.p1.size() ){
            State cpy = s;
            cpy.stack.push_back({cpy.p1.back(),s.turn});
            cpy.p1.pop_back();
            cpy.passed = 0;
            cpy.turn ^= 1;
            ans = max(dfs(cpy),ans);
        }
        {
            State cpy = s;
            int profit = 0;
            int coef1 = 1;
            int coef2 = 1;
            for(int i = cpy.stack.size() - 1 ; i >= 0 ; i--){
                if( cpy.stack[i].second != s.turn &&  cpy.stack[i].first == -1 ) coef1 = 0;
                if( cpy.stack[i].second == s.turn &&  cpy.stack[i].first == -1 ) coef2 = 0;
                
                if( cpy.stack[i].second == s.turn &&  cpy.stack[i].first != -1 )
                    profit +=  coef1 * cpy.stack[i].first;
                if( cpy.stack[i].second != s.turn &&  cpy.stack[i].first != -1 )
                    profit -=  coef2 * cpy.stack[i].first;
            }
            if( cpy.stack.size() == 0 ) cpy.passed++;
            cpy.stack.clear();
            cpy.turn ^= 1;
            ans = max(dfs(cpy)+profit,ans);
        }
    }else{
        ans = +1e9;
        if( s.p2.size() ){
            State cpy = s;
            cpy.stack.push_back({cpy.p2.back(),s.turn});
            cpy.p2.pop_back();
            cpy.passed = 0;
            cpy.turn ^= 1;
            ans = min(dfs(cpy),ans);
        }
        {
            State cpy = s;
            int profit = 0;
            int coef1 = 1;
            int coef2 = 1;
            for(int i = cpy.stack.size() - 1 ; i >= 0 ; i--){
                if( cpy.stack[i].second != s.turn &&  cpy.stack[i].first == -1 ) coef1 = 0;
                if( cpy.stack[i].second == s.turn &&  cpy.stack[i].first == -1 ) coef2 = 0;
                
                if( cpy.stack[i].second == s.turn &&  cpy.stack[i].first != -1 )
                    profit +=  coef1 * cpy.stack[i].first;
                if( cpy.stack[i].second != s.turn &&  cpy.stack[i].first != -1 )
                    profit -=  coef2 * cpy.stack[i].first;
            }
            if( cpy.stack.size() == 0 ) cpy.passed++;
            cpy.stack.clear();
            cpy.turn ^= 1;
            ans = min(dfs(cpy)-profit,ans);
            
        }
    }
    return memo[s] = ans;
    
}
 
int main(){
    int n,m;
    cin >> n >> m;
    vector<int> a(n), b(m);
    for(int i = 0 ; i < n ; i++) cin >> a[i];
    for(int i = 0 ; i < m ; i++) cin >> b[i];
    reverse(a.begin(),a.end());
    reverse(b.begin(),b.end());
    
    State s = {a,b,{},0,0};
    cout << dfs(s) << endl;
    
}

JAG Contest 2016 Domestic B - 豪邸と宅配便

ICPC系

jag2016-domestic.contest.atcoder.jp

問題文

省略

解法(反転して表示)

時刻 a に宅配便が届くと,[a-M,a+M)の範囲に書斎には居られない. なので,[0,T)の範囲で区間に含まれていない時刻を列挙すれば良い. 境界に気をつけつつ,愚直にシミュレーションすると時間計算量 O(NM+T) で解ける.いもす法を使うと O(N+T) になるけど制約的にどうでもいい.

ソース

#include <bits/stdc++.h>
using namespace std;
 
 
int imos[20100];
int main(){
    int N,M,T;
    cin >> N >> M >> T;
    for(int i = 0 ; i < N ; i++){
        int a;
        cin >> a;
        for(int j = 0 ; j <= M ; j++)
            if( a - j >= 0 ) imos[a-j] = true;
        for(int j = 0 ; j < M ; j++)
            if( a + j < T ) imos[a+j] = true;
    }
    
    int ok = 0;
    for(int i = 0 ; i < T ; i++){
        // cout << imos[i] << " " << i << endl;;
        if( !imos[i] ){
            ok++;
        }
    }
    cout << ok << endl;
}

JAG Contest 2016 Domestic A - 阿吽の呼吸

ICPC系

jag2016-domestic.contest.atcoder.jp

問題文

省略

解法(反転して表示)

本質的には括弧列が与えられるのでvalidな括弧列ですかという問題に言い換えれば見通し良く解ける. カウンタを定義する.「A」が来たらカウンタをインクリメント,「Un」が来たらデクリメント,途中でカウンタが負の数になったらNO,最後まで処理してカウント!=0ならNO.

ソース

#include <bits/stdc++.h>
using namespace std;
 
int main(){
    int a = 0;
    int n;
    cin >> n;
    for(int i = 0 ; i < n ; i++){
        string t;
        cin >> t;
        if( t == "A" ){
            a++;
        }else{
            a--;
            if( a < 0 ){
                cout << "NO" << endl;
                return 0;
            }
        }
    }
    if( a == 0 ) cout << "YES" << endl;
    else cout << "NO" << endl;
}