Hatena::Grouptopcoder

cafelier@SRM

cafelier のSRM参加記録です。コンテスト中に考えてたことを執拗に全部書き残すとどうなるだろうかという試み本番中にこういうコードが書きたかったなあ、という後で書いた反省コードを書き残す試み

スパムが来たのでしばらくコメント欄をはてなユーザ限定にしています、すみません、

2013-11-24

ACM-ICPC 2013 Aizu Regional 問題I

| 00:36 | はてなブックマーク -  ACM-ICPC 2013 Aizu Regional 問題I - cafelier@SRM

ACM-ICPC 2013 アジア地区予選会津大会 の問題Iを解いてみたよ。

まず、balanced tree はなんであれ、深さ d の葉の重みが w だったら深さ d-1 の重みは 2w、深さ d-2 の重みは 4w、d-3 は 8w、…と倍々に増えていきます。なので例えば、同じ木に重み 1 と 3 の葉が混ざることは絶対にありえない。というわけで、まず最初の列を 1,2,4,8,... だけ含む列と、3,6,12,24,... だけ含む列と、5,10,20,…と、等々に分離してそれぞれ別々に解くことができます。つまり、一般性を失わず、配列の要素が1,2,4,8,...だけの場合のみ考えれば良い。

で、1,2,4,8,... だけの場合、まず過程をすっ飛ばして答えだけ先に書くと、こんな感じに更新式が +1 と max だけの非常にシンプルな DP になります。

#include <iostream>
#include <vector>
#include <map>
#include <numeric>
using namespace std;

// Aの要素が全部 {1,2,4,8,16,...} の場合を解く
int solve_1248(const vector<int>& A)
{
   int sum = accumulate(A.begin(), A.end(), 0);

   vector<int> dp(sum+1, -999999999);
   dp[0] = 0;
   for(int Ai: A)
      for(int x=sum/Ai*Ai; x>=Ai; x-=Ai) // (Aiの倍数すべてについて)
         dp[x] = max(dp[x], dp[x-Ai]+1);

   int best = 0;
   for(int x=1; x<=sum; x*=2)
      best = max(best, dp[x]);
   return best;
}

// A を {1,2,4,8,...} と {3,6,12,24,...} と {5,10,20,40,...} に分けてそれぞれ解く
int solve(const vector<int>& A)
{
   map<int, vector<int>> bucket;
   for(int Ai: A)
   {
      int lsb = Ai & (~Ai+1);  // (Aiの1が立ってる一番下のビットだけ残す)
      bucket[Ai/lsb].push_back(lsb);
   }

   int best = 0;
   for(const auto& kv: bucket)
      best = max(best, solve_1248(kv.second));
   return best;
}

// 入出力
int main()
{
   for(int N; cin>>N, N; )
   {
      vector<int> A(N);
      for(auto& Ai: A)
         cin >> Ai;
      cout << solve(A) << endl;
   }
}

計算量は O(N・ΣA) で、N≦1000、Ai≦500 なので少し大きめに見えますが、ループの中身がものすごく単純なので余裕です。(※あるいは、上に述べたように2の冪乗だけ考えればいいので、Aiの実質の上限は256で、ΣA は 256000 ですが、これも本当は2の冪乗までしか要らないのでつまり 131072 まで見ればよい、等々細かく解析していくと意外と安心できる演算回数が上限として計算できます。)

solve_1248() の導出

上のコードは、実はこういうことをやっています。配列 A を左から見ていって、balanced tree を左からじわじわと作っていく。

int solve_1248(const vector<int>& A)
{
   map<作りかけの木, その木を作るのに最大で使える葉の数> dp;
   dp[空集合] = 0;

   for(int Ai : A)
      for((partial_tree, leaf_count) : dp)
         if(partial_tree の右に葉 Ai をくっつけることができる) {
            ptree2 = partial_tree の右に葉 Ai をくっつけたもの
            dp[ptree2] = max(dp[ptree2], leaf_count+1);
         }

   int best = 0;
   for((ptree, leaf_count) : dp)
      if(ptreeが木として完成している)
         best = max(best, leaf_count);
   return best;
}

作りかけの木、や、右から葉を足す、というのはこんな感じ:

f:id:cafelier:20131125000530p:image:w300

f:id:cafelier:20131125000752p:image:w450

で、こんな木っぽいものを dp のキーにするのは大変そうですしそもそも木の形って指数的にパターンがありすぎて爆発しそうな気が一瞬しますが、そうでもありません。この木はあとで「右から葉を足す」操作しか行わないので、「右から葉を足す」観点で違いがない木たちは同一視して同じものと思って構いません。具体的には、完全に二分木として完成している部分は重みを合計して1ノードにつぶしてしまって差し支えない。

f:id:cafelier:20131125001102p:image:w400

で、こうやって潰した作りかけ木だけを考えると、葉ノードは、常に、左から右に狭義単調減少しているものだけが現れます。上の図なら[8,4]という減少列。なんでかというと減少しないで同じ値が並んでたらそこは木が完成していて潰せるので。同じ並びがないなら、増えるところがあるとbalanced treeにならない。

で、で、[8,4] や [16,2,1] や [8,4,2,1] のような2の冪の狭義減少列は、自然数の二進数表記と1対1に対応します。[8,4]は8+4=12という自然数を、[16,2,1]は19=16+2+1という自然数を表していると思えばよい。逆に自然数から列は一意に復元できます。

つまり、長々と書いてきましたが、ここまで書いたことをまとめると、「作りかけの木」は、「葉の重みを全部足した自然数」で表現できます。ので、それをDPのキーにすれば簡単です。

int solve_1248(const vector<int>& A)
{
   int sum = accumulate(A.begin(), A.end(), 0);

   vector<int> dp(sum+1, -999999999);
   dp[0] = 0;
   for(int Ai: A)
      // 重み総和x-Aiの作りかけの木に、右から葉 Ai をくっつけて重み総和xの木を作る。
      // 右から葉Aiをくっつけるには、単調減少列をキープするには、
      // x-Aiの二進数で一番右の桁がAi以上でないといけない。
      // つまりx-Aiが、あるいはxがAiの倍数でないと。
      for(int x=sum/Ai*Ai; x>=Ai; x-=Ai)
         dp[x] = max(dp[x], dp[x-Ai]+1);

   // 完成している木の重みは2のベキ乗なのでそこのスコアを数える
   int best = 0;
   for(int x=1; x<=sum; x*=2)
      best = max(best, dp[x]);
   return best;
}

以上。

vector<int> dp;という配列を作るのではなくて、map<int,int> や vector<pair<int,int>> を使って実際に作れる可能性のあるsumだけ覚える、というコードを考えた人もいるかもしれませんが、可能性のあるsumは簡単に全通りできるので、この表はちょっとでも意地の悪いデータなら全エントリ埋まります。普通の配列でベタに持つ方が圧倒的に効率がよいです。

solve_1248() の導出:別の見方

この問題は構文解析問題だと思っています。

1. balanced treeを表現する文法を定義する。こういうBNFですね。

2 ::= 1 1
4 ::= 2 2
8 ::= 4 4
... 

2. LR(0)構文解析器のようなボトムアップの構文解析器を作る

3. 使わない数値は捨てるという処理が必要なので非決定的な構文解析になりますが、その非決定性を全探索…は非効率的なので

4. ボトムアップ構文解析のスタックをキーにしてメモ化する

5. スタックの要素は必ず [8,4,1] 的な単調減少列になっているので、(上でも書いた二進数と対応する理論により)高々、作れる最大ツリーの重みくらいの可能性しかないので十分なメモ化が利く

トラックバック - http://topcoder.g.hatena.ne.jp/cafelier/20131124

2013-07-13

ACM-ICPC 2013 Tokyo Domestic 問題F

| 20:49 | はてなブックマーク -  ACM-ICPC 2013 Tokyo Domestic 問題F - cafelier@SRM

ACM-ICPC 2013 国内予選 の問題Fを解いてみたよ。

O(n^3m^3 + (n^3+m^3)r|x|).

#include <iostream>
#include <vector>
#include <set>
#include <tuple>
using namespace std;

// x∈dp[i][j]  ⇔  A[i..j) を書き換えまくったら x 1つにできる
vector<vector<set<int>>> make_rewritability_table(
  const vector<int>& A, const vector<tuple<int,int,int>>& R)
{
  const int N = A.size();
  vector<vector<set<int>>> dp(N, vector<set<int>>(N));
  for(int i=0; i<N; ++i)
    dp[i][(i+1)%N].insert(A[i]);

  for(int k=2; k<=N; ++k)
  for(int i=0; i<N; ++i)
  for(int m=1; m<k; ++m)
  for(auto& rule : R) {
    int x1,x2,y; tie(x1,x2,y) = rule;
    if(dp[i][(i+m)%N].count(x1) && dp[(i+m)%N][(i+k)%N].count(x2))
      dp[i][(i+k)%N].insert(y);
  }
  return dp;
}

// 高速化のため set<int> を int のビット表現に変換
vector<vector<int>> as_bitset(const vector<vector<set<int>>>& table)
{
  vector<vector<int>> bit_table(table.size(), vector<int>(table[0].size()));
  for(int i=0; i<table.size(); ++i)
  for(int k=0; k<table[i].size(); ++k)
  for(int v : table[i][k])
    if(v > 0)
      bit_table[i][k] |= 1<<v;
  return bit_table;
}

// A[ao]が先頭に来るようにAを回転したもの vs B[bo]以下略
// を書き換えまくって作れる最長共通列を考える。
//
// A[as..ae] と B[bs..be) が同じ数 1 つに書き換えられるときに
// (as,bs)==>(ae,be) という長さ 1 の辺を引いた DAG を考えて、
// DAG の最長経路を求めればよい。
int solve_rotated(const vector<vector<int>>& rA, int ao,
                  const vector<vector<int>>& rB, int bo)
{
  const int N = rA.size(), M = rB.size();
  vector< vector<int> > dp(N+1, vector<int>(M+1, -99999));
  dp[N][M] = 0;
  for(int ae=N; ae>=1; --ae)
  for(int be=M; be>=1; --be)
  for(int as=0; as<ae; ++as)
  for(int bs=0; bs<be; ++bs)
    if(rA[(as+ao)%N][(ae+ao)%N] & rB[(bs+bo)%M][(be+bo)%M])
      dp[as][bs] = max(dp[as][bs], 1 + dp[ae][be]);
  return dp[0][0];
}

// 回転は最初に全部やれば十分なので、
// 全通りの回し方を考えて solve_rotated に投げる
int solve(const vector<int>& A, const vector<int>& B,
          const vector<tuple<int,int,int>>& R)
{
  vector<vector<int>> rA = as_bitset(make_rewritability_table(A, R));
  vector<vector<int>> rB = as_bitset(make_rewritability_table(B, R));

  int best = -1;
  for(int ao=0; ao<A.size(); ++ao)
  for(int bo=0; bo<B.size(); ++bo)
    best = max(best, solve_rotated(rA, ao, rB, bo));
  return best;
}

// 入出力
int main()
{
  for(int N,M,R; cin>>N>>M>>R, N|M|R; )
  {
    vector<int> A(N); for(int& Ai: A) cin>>Ai;
    vector<int> B(M); for(int& Bi: B) cin>>Bi;

    vector<tuple<int,int,int>> P;
    for(int i=0,ID=-1; i<R; ++i)
    {
      int K; cin>>K;
      vector<int> X(K); for(int& Xi: X) cin>>Xi;
      int Y; cin>>Y;

      // 左辺が長いルールは、左辺が長さ2になるように分解する
      // 例: a b c d -> e ==> {a b -> -1,  -1 c -> -2,  -2 d -> e}
      for(int k=0; k+1<K; ++k)
        P.emplace_back(k?ID-k:X[0], X[k+1], k==K-2?Y:ID-k-1);
      ID -= K-2;
    }

    cout << solve(A, B, P) << endl;
  }
}

各所で指摘されているけれど make_rewritability_table() は、CYK法 と呼ばれる構文解析アルゴリズムですね。全部分区間の処理が済んでいれば全体の処理も書き換え1回分だけ考えればいい、という区間DP

入力列の長さに関する影響だけを考えて O(n^3) であると書かれることが多いのですけど、よほど実装を工夫でもしない限り、書き換え規則のサイズ|G|まで考慮にいれるとこれが掛かって O(n^3・|G|) なので、回転の扱いをさぼるために solve_rotated の中で毎回計算し直したりすると意外と時間がかかったりする。

トラックバック - http://topcoder.g.hatena.ne.jp/cafelier/20130713

2012-12-02

ACM-ICPC 2012 Tokyo Regional 問題F, I

| 23:42 | はてなブックマーク -  ACM-ICPC 2012 Tokyo Regional 問題F, I - cafelier@SRM

ACM-ICPC 2012 アジア地区予選 の問題FとIを解いてみたよ。

F

  • 「物体xと物体yの重さの差はwグラムです」という形式の報告と「物体xと物体yの重さの差はいくつ?」という形式の質問がたくさん来ます。質問にそのたびに答えを返しなさい。
    • x+5=y, y+10=z みたいなのが報告済みなら「xとzの差は?」に15と答えないといけません。
    • わからないときは "UNKNOWN" と答えれてよい。
#include <iostream>
#include <vector>
#include <utility>
using namespace std;

typedef int Node, Weight;

struct UnionFind
{
   vector< pair<Node, Weight> > uf; // <parent, offset_from_parent>

   UnionFind(int N)
   {
      for(int i=0; i<N; ++i)
         uf.push_back(make_pair(i, 0));
   }

   // Returns <root, offset_from_root> while compressing the path.
   pair<Node, Weight> Find(int a)
   {
      if(uf[a].first != a) {
         pair<Node, Weight> p = Find(uf[a].first);
         uf[a] = make_pair(p.first, uf[a].second+p.second);
      }
      return uf[a];
   }

   void Union(int a, int b, Weight b_minus_a) // [a] + b_minus_a = [b]
   {
      pair<Node, Weight> aa = Find(a); // [aa.first] + aa.second = [a]
      pair<Node, Weight> bb = Find(b); // [bb.first] + bb.second = [b]
      // Then, [bb.first] + (bb.second - b_minus_a - aa.second) = [aa.first]
      if( aa.first != bb.first )
         uf[aa.first] = make_pair(bb.first, bb.second - b_minus_a - aa.second);
   }
};

int main()
{
   for(int N,Q; cin>>N>>Q,N|Q; )
   {
      UnionFind uf(N);
      while( Q --> 0 ) {
         char cmd; int a,b; cin>>cmd>>a>>b; --a,--b;
         if(cmd == '!') {
            int w; cin>>w;
            uf.Union(a,b,w);
         } else {
            pair<Node, Weight> aa = uf.Find(a), bb = uf.Find(b);
            if(aa.first != bb.first)
               cout << "UNKNOWN" << endl;
            else
               cout << (bb.second - aa.second) << endl;
         }
      }
   }
}
  • UNKNOWNかそうでないかは、報告された関係で繋がっているかどうかでわかるので、同値類分類、いわゆる Union-Find データ構造で解けます。
  • で、KNOWN の時に具体的な差の値を返さないといけないんですが、こういう類いの"差"の情報はUnion-Findのedgeに乗せてやることができます。
    • 大抵の人は Union-Find は同値類関係を判定するだけの分類処理としてライブラリ化していると思うので、それを果たしてちゃんと拡張できるかどうか、的な問題

union-findの枝に情報を載せるんじゃなくて、各物体の「仮の」重さを決めておいて、factの報告が来て矛盾していたら、依存関係の小さい方の「仮の」重さを調整する、みたいな方法もありです(Union-Findのバランス木を使ってO(n log n)にする方法と実質的に同じです。)

I

  • 単語の列を決まった幅の紙面にレイアウトしたい。単語間空白の最大幅ができるだけ小さくなるように改行位置を決定せよ。

これは「最大値の最小化にぶんたんさく」というルビを振って二分探索で解いてもいいんですがとりあえず O(n) 解。

#include <iostream>
#include <vector>
#include <stack>
using namespace std;

int solve(const int W, const vector<int>& x)
{
   const int N = x.size();

   // xsum[k] := Σ x[0..k)
   vector<int> xsum(1, 0);
   for(int k=0; k<N; ++k)
      xsum.push_back(xsum.back()+x[k]);

   // DP高速化のための準備
   int line_start = 0;
   vector<int> next_smaller(N, 0);
   stack<int> no_next_smaller; no_next_smaller.push(0); // (番兵)

   // dp[k] := 単語 x[k] が行頭に来るように x[0..k) を最適配置した場合のそこまでの最適解
   vector<int> dp(N); dp[1] = W;
   for(int k=2; k<N; ++k)
   {
      // line_start := 直前の行の行頭になり得る一番最初の単語
      while(xsum[k]-xsum[line_start]+(k-line_start-1) > W)
         ++line_start;

      dp[k] = W;
      for(int ls=line_start; ls+1<k;) {
         // s := x[ls..k) を一行に配置したときの、その行の最大連続空白
         int s = (W-(xsum[k]-xsum[ls])-1)/(k-ls-1) + 1;

         // DP表の更新
         dp[k] = min(dp[k], max(dp[ls],s));
         if(dp[ls] <= s) {
            // lsを増やすとsは広義単調増加するので、
            // これ以上回してもdp[k]が更新されることはない。おしまい
            break;
         } else {
            // ここに来るとき、dp[line_start..ls) >= dp[ls] > s である。
            // さて、x[ls..k) より x[ls..k+1) の方が一行に配置したときの空白"s"は小さい。
            // なので、次周で dp[k+1] を埋めるときも max(dp[ls],s) の値を支配するのは
            // dp[ls] の方なので、dp[line_start..ls) のゾーンは見る必要がない
            line_start = ls;

            // lsを増やすとsは広義単調増加するので、
            // dp[ls]が小さくなるとこまで行かないとdp[k]は更新されない。じゃんぷ
            if( !(ls=next_smaller[ls]) )
               break; // なかった
         }
      }

      // next_smaller の更新
      for(; dp[no_next_smaller.top()]>dp[k]; no_next_smaller.pop())
         next_smaller[no_next_smaller.top()] = k;
      no_next_smaller.push(k);
   }

   // 最終行が単語 x[k] から始まる場合
   int best = W;
   for(int k=0; k<N; ++k)
      if( xsum[N]-xsum[k]+(N-k-1) <= W ) {
         int s = (k+1==N ? 0 : 1);
         best = min(best, max(dp[k], s));
      }
   return best;
}

int main()
{
   for(int W, N; cin>>W>>N, W|N; ) {
      vector<int> x(N);
      for(int i=0; i<N; ++i)
         cin >> x[i];
      cout << solve(W, x) << endl;
   }
}

next_smallerなしでもなんとかなる気もするんですが、はてさて。

まあでも、そんなに頑張らなくても二分探索+しゃくとり法で O(|x| log W) で解けます。

#include <iostream>
#include <vector>
using namespace std;

bool is_possible(const int W, const vector<int>& x, const int MaxSpace)
{
   // CanStartLine[k] := 単語 x[k] から後ろを空白幅MaxSpace以下でレイアウトできますか
   vector<bool> CanStartLine(x.size());

   int CanStartCnt = 0;
   int TightWidth = -1,        TightNextStart = x.size();
   int LooseWidth = -MaxSpace, LooseNextStart = x.size();

   for(int k=x.size()-1; k>=0; --k) {

      // LooseNextStart := (x[k]+MaxSpace+...+MaxSpace+x[e] >= W になる最小の e+1)
      // LooseWidth     := その時の x[k]+MaxSpace+...+MaxSpace+x[e]

      LooseWidth += x[k] + MaxSpace;
      while( W <= LooseWidth - (x[LooseNextStart-1]+MaxSpace) ) {
         LooseWidth  -= x[--LooseNextStart] + MaxSpace;
         CanStartCnt += (CanStartLine[LooseNextStart] ? 1 : 0);
      }

      // TightNextStart := (x[k]+1+...+1+x[e] > W となる最小の e+1)
      // TightWidth     := その時の x[k]+1+...+1+x[e]

      TightWidth += x[k] + 1;
      while( W < TightWidth - (x[TightNextStart-1]+1) ) {
         TightWidth  -= x[--TightNextStart] + 1;
         CanStartCnt -= (CanStartLine[TightNextStart] ? 1 : 0);
      }

      // x[k] が最初の1単語だとすると、
      // 次の行の先頭に来れる単語は [LooseNextStart, TightNextStart) のどれか&どれでも。
      // CanStartCnt := その範囲にある CanStartLine の true の個数。

      CanStartLine[k] = TightWidth<=W || CanStartCnt>0;
   }

   return CanStartLine[0];
}

int solve(const int W, const vector<int>& x)
{
   // 許す最大空白幅を二分探索
   int L=0, R=W; // (L,R]
   while( R-L>1 )
      (is_possible(W, x, (L+R)/2) ? R : L) = (L+R)/2;
   return R;
}

int main()
{
   for(int W, N; cin>>W>>N, W|N; ) {
      vector<int> x(N);
      for(int i=0; i<N; ++i)
         cin >> x[i];
      cout << solve(W, x) << endl;
   }
}

しゃくとりの替わりに累積和上の二分探索やBITなどなどでもう一個 O(log W) を掛けても間に合うようです。

別解。計算量的に謎だけどジャッジデータは余裕で通る。

#include <iostream>
#include <vector>
using namespace std;

bool is_possible(const int W, const vector<int>& x, const int MaxSpace)
{
   // x[i] が置かれることがあり得る column の全リスト
   vector< pair<int,int> > range; // 重なりのない閉区間のリストとして保持
   bool CanStartLine = true;      // ただし [0,0] だけ特別扱いする

   for(int i=0; i+1<x.size(); ++i) {
      // x[i+1] の可能性をさぐる
      vector< pair<int,int> > range2;
      bool CanStartLine2 = false;
      for(int k=(CanStartLine ? -1 : 0); k<(int)range.size(); ++k)
      {
         int L = (k==-1 ? 0 : range[k].first);
         int R = (k==-1 ? 0 : range[k].second);

         // x[i] で行を終えることが可能なら x[i+1] は行頭に来れる
         if( L<=W-x[i] && W-x[i]<=R )
            CanStartLine2 = true;

         // x[i]が[L,R]に置かれるなら次の単語は[L+x[i]+1, R+x[i]+MaxSpace]である
         L += x[i] + 1;
         R += x[i] + MaxSpace;
         R = min(R, W-x[i+1]); // でもあんまり右には置けないよ
         if( L <= R )
            range2.push_back(make_pair(L,R));
      }

      // rangeがLでソート済みならrange2もそうであることに注意。
      // range2に突っ込まれた区間をできるだけマージする
      CanStartLine = CanStartLine2;
      range.clear();
      for(int k=0,e; k<range2.size(); k=e)
      {
         int L = range2[k].first;
         int R = range2[k].second;
         for(e=k+1; e<range2.size(); ++e)
            if(range2[e].first <= R+1)
               R = max(R, range2[e].second);
            else
               break;
         range.push_back(make_pair(L,R));
      }
   }

   // どこでもいいけど最後の単語がどっかに置ける可能性があるならOK
   return (CanStartLine || !range.empty());
}

int solve(const int W, const vector<int>& x)
{
   // 許す最大空白幅を二分探索
   int L=0, R=W; // (L,R]
   while( R-L>1 )
      (is_possible(W, x, (L+R)/2) ? R : L) = (L+R)/2;
   return R;
}

int main()
{
   for(int W, N; cin>>W>>N, W|N; ) {
      vector<int> x(N);
      for(int i=0; i<N; ++i)
         cin >> x[i];
      cout << solve(W, x) << endl;
   }
}

先頭行から割と素直に単語を配置していく感じですが、MaxSpace=1 の時は配置が一意に決まりますし、MaxSpace>1 の時は行の後ろの方にいくにつれ可能な区間の幅が伸びていくことを考えると、区間のdisjointなリスト range の長さは高々 √W にしかならないので、O(|x|・√W・log W) には少なくともなっていて、実際にはもっと低く抑えられる気がしています。

トラックバック - http://topcoder.g.hatena.ne.jp/cafelier/20121202

2010-12-16

ACM-ICPC 2010 Tokyo Regional 問題F, G

| 21:12 | はてなブックマーク -  ACM-ICPC 2010 Tokyo Regional 問題F, G - cafelier@SRM

ACM-ICPC 2010 アジア地区予選 の問題FとGを解いてみたよ。

F

  • 『素数 Q と、0-9 の数字が N 個並んだ列がある。その中で十進数として読んだら Q の倍数になる区間は何個?』
    • Q ≦ 1億
    • N ≦ 10万
  • Q=3 で 12345 なら 12, 123, 12345, 234, 3, 345 の6個
  • ゼロで始まる区間はカウントしません

  • まず入力の最大サイズを見ると
    • O(N3) や O(NQ) や O(N2) の時間計算量では時間かかりすぎです
    • O(N) か O(N log N) か。
easy case
  • とりあえず Q が 2 か 5 なら簡単で
  • 「下一桁が0,2,4,6,8な区間を全部数える」や「下一桁が0,5な区間を全部数える」
  • それだけですね。
int solve(const vector<int>& a, int Q)
{
   if( Q==2 || Q==5 )
      return easyCase(a, Q);
   else
      return hardCase(a, Q);
}

int easyCase(const vector<int>& a, int Q)
{
   int answer = 0;
   {
      int nonZero = 0;
      for(int i=0; i<a.size(); ++i)
      {
         nonZero += (a[i]!=0 ? 1 : 0);
         answer  += (a[i]%Q==0 ? nonZero : 0);
      }
   }
   return answer;
}
  • Qの倍数を見かけたら、それより左にあったゼロじゃない桁の個数を足し合わせ。

hard case
  • 問題は、Q が 2 でも 5 でもないとき。
  • どうしましょう。問題を数式で整理すると
    • 10k*a[i] + 10k-1*a[i+1] + ... + 101*a[i+k-1] + a[i+k]
    • が Q の倍数になってるような i, k は全部で何個?
  • でした。

  • 10k という余計なものが掛け算されてますが無視すると
  • これは、a[i] + a[i+1] + ... + a[i+k-1] + a[i+k] という 配列の区間の和の話をしています。

  • 区間和を扱うときにはよくでてくるアイデアがあって
    • b[0] = a[0]
    • b[1] = b[0] + a[1] (= a[0]+a[1])
    • b[2] = b[1] + a[2] (= a[0]+a[1]+a[2])
    • ...
  • まず、こういう風に「端からの区間和」だけを全部求めておくと
    • a[i]+a[i+1]+...+a[i+k] = b[i+k] - b[i-1]
  • あとは、どこでも好きな場所の区間和を引き算一発で求められます。

  • これに似てるので、何か応用できないでしょうか。
  • 問題には、単なる区間和ではなくて、
  • 左にいくほど10が沢山かかった重み付け区間和が出てきているので
    • a[0], a[1], a[2], ..., a[N-2], a[N-1]
  • これにとりあえず10を重みづけて掛けてみる。
    • 10N-1*a[0], 10N-2*a[1], 10N-3*a[2], ..., 10*a[N-2], a[N-1]
  • これにさっきの区間和計算を適用すると…?
  • 欲しかった
    • 10k*a[i] + 10k-1*a[i+1] + ... + 101*a[i+k-1] + a[i+k]
    • (10k*a[i] + 10k-1*a[i+1] + ... + 101*a[i+k-1] + a[i+k]) * 10N-k-1
  • 余計な 10N-k-1 が掛け算された和なら、引き算一発で求まっちゃいます。

  • 実はこれさえわかれば十分なのです。
  • 「Q の倍数かどうか」しか気にしないので。
  • Q が2でも5でもない素数なら、10を掛けても掛けなくても Qの倍数かどうかはかわりません!

  • というわけで
    • a[] に 10 を重み付け掛け算
    • b[] = 区間和計算用の「端からの和」
    • b[i+k] - b[i-1] が Q の倍数になる i, k の組の個数を数える
  • とやれば答えが求まります。
  • 最後の組の個数を数える部分は2重ループすると時間が足りませんが、
  • b[i] の値を Q で割った余りで分類してヒストグラム作っておけば1重ループで計算できます
int hardCase(const vector<int>& a, int Q)
{
   int answer = 0;
   {
      map<int,int> freq;
      freq[0] = 1;

      int base = 1, sum = 0;
      for(int i=a.size()-1; i>=0; --i)
      {
         sum  = (a[i]*base + sum) % Q;
         base = (base*10) % Q;

         if( a[i] != 0 )
            answer += freq[sum];
         freq[sum] ++;
      }
   }
   return answer;
}
  • というのを1回のループに全部まとめると、こうなりました。
advanced case
  • 実は
    • 『Q が素数とは限らない』
  • という問題設定でも、解けます。
  • easyCase と hardCase の合わせ技。お暇な方はチャレンジしてみては。


G

  • 『「有向グラフの最短経路長を求めよ」という問題をコンテストで出題しようと思ったんだけど、テストデータをランダム生成したら、面白くない答えばっかりになってしまいました。"2010" や "42" や "6502" のような意味ありげな答えになるネタを仕込みたいです。頑張ろう。』
    • 入力は
      • 点数 100 以下、辺数 1000 以下の、辺に非負重みが乗った有向グラフ
      • その頂点 s と t
      • 目標の最短経路長 c
    • 出力は
      • sからtへの最短経路長をピッタリ c にするために書き換えなきゃ行けない辺の本数の最小値
    • 条件として
      • 「元々の最短経路長より c の方が小さい」場合しか来ないそうです。

  • 「ピッタリ c」と考えると難しいですが
  • 「最短経路長を c 以下にする」という目標だとどうでしょうか。
  • そうすると、辺の長さを書き換える時は徹底的に短く変えた方がお得になるので
  • 0 に変えるケースだけ考えればよくなります。

  • というわけで
    • (頂点, 辺を0に変える操作を使った回数)
  • を頂点にした拡大グラフで最短路を求めればよし。
typedef int vert, cost, modify;
typedef vector< vector< pair<vert,cost> > > graph;

modify dijkstra(const graph& G, cost Ct, vert Src, vert Dst)
{
   modify mMin = G.size();

   typedef pair<vert, modify>  state;
   typedef pair<cost, state> c_state;
   priority_queue< c_state, vector<c_state>, greater<c_state> > Q;
   set<state> V;

   Q.push( c_state(0, state(Src,0)) );
   while( !Q.empty() )
   {
      cost   c = Q.top().first;
      vert   v = Q.top().second.first;
      modify m = Q.top().second.second;
      Q.pop();

      if( c>Ct || m>mMin || V.count(state(v,m)) )
         continue;
      V.insert(state(v,m));

      if( v==Dst && c<=Ct )
         mMin = min(mMin, m);

      for(int i=0; i<G[v].size(); ++i)
      {
         vert u  = G[v][i].first;
         cost cc = G[v][i].second;
         c_state next[2] = {c_state(c+cc,state(u,m)), c_state(c,state(u,m+1))};
         for(int j=0; j<2; ++j)
            if( !V.count(next[j].second) )
               Q.push(next[j]);
      }
   }
   return mMin;
}

int main()
{
   fstream fin("G.in");
   for(int Nv,Ne,Ct; fin>>Nv>>Ne>>Ct,Nv|Ne|Ct; )
   {
      graph G(Nv);
      for(int i=0; i<Ne; ++i)
      {
         vert f,t;
         cost c;
         cin >> f >> t >> c;
         G[f-1].push_back( make_pair(t-1,c) );
      }
      cout << dijkstra(G, Ct, 0, Nv-1) << endl;
   }
}
  • この実装よりも、
  • 「書き換え回数が少ないノードを先に探索して、距離c以下でtに辿り着けたらそこで探索打ち切り」
  • という方が多分効率的ですかね。

  • 実はこれで問題は解けていて、
  • というのは「最短経路長をc以下にできる」なら
  • 減らしすぎた部分をちょびちょび増やしていけば必ず「ピッタリcにもできる」ので。
  • 難しく言うと、最短経路長は辺の長さに関する連続関数なので、中間値の定理によりOKです。

advanced case
  • 『c は元々の最短経路長より長い』とすると、どうなるでしょう。
    • 同じように「c 以上にする」と考えて、
    • 「辺の長さを∞に変えるパターンだけ」扱うのは同じなのですが
    • その先は…?
  • 私も答えがわからないので、どなたかすごい人考えてみて下さい
  • 追記: 投機的撃墜例。最適解は真ん中の2本変えればいいので2です。

f:id:cafelier:20101218094913p:image

トラックバック - http://topcoder.g.hatena.ne.jp/cafelier/20101216

2010-07-03

ACM-ICPC 2010 Japan Domestic 問題E

| 02:08 | はてなブックマーク -  ACM-ICPC 2010 Japan Domestic 問題E - cafelier@SRM

ACM-ICPC 2010 国内予選 の問題Eを解いてみたよ。



解法1:最良優先

問題に書いてあることに忠実に、スター(ト)からゴール(ド)に向かって、辞書順で早い方早い方を優先して進んで行く方法です。注意しないといけないのは「何も考えずに、まっしぐらに辞書順で一番早いラベルの矢印だけを選んで進む」のは間違いということです。二つ理由があって…


  • 辞書順最小の方に進むとゴールにたどりつけなくなる場合
3 2 0 2
0 1 a
0 2 b
    • 答えは "b" ですが、スタート 0 から出てる最速の矢印は "a" だからといってそっちに行っちゃうと未来がありません。行き止まり。
    • 解決策としては、本体の探索を始める前に、各節点が「行き止まり」かどうか判定しといて、そっちには決して行かないようにするなど

  • 接頭辞に要注意
3 3 0 2
0 1 a
0 2 ab
1 2 c
    • 上の例の答えは "ab" ですが、スタート 0 から生えてる矢印のラベルだけで比べると、"a" < "ab" なので、うっかり "a" 方向に突き進んで "ab" のことを忘れてしまうと、間違えて "ac" を答えてしまいます。
    • 解決策としては、矢印一本だけ見ると遅くみえるルートでも「探索候補」としてとりあえず覚えておいて、必要になったらとりだす感じ。

言葉で説明するよりもコード見る方がわかりやすいですね。

// Java
import java.util.*;

public class PowerfulSpell
{
   public static void main(String[] arg) throws java.io.IOException
   {
      // データセットを1個読み込んでは solve() を呼び出すループ
      for(Scanner sc = new Scanner(System.in) ;;)
      {
         int N = sc.nextInt();
         int A = sc.nextInt();
         int S = sc.nextInt();
         int G = sc.nextInt();
         if( N==0 && A==0 && S==0 && G==0 )
            break;

         Arrow[] arr = new Arrow[A];
         for(int i=0; i<A; ++i)
            arr[i] = new Arrow(sc);

         System.out.println( solve(N,S,G,arr) );
      }
   }

   // 矢印データ
   static class Arrow
   {
      final int    x;
      final int    y;
      final String lab;

      Arrow(Scanner sc) throws java.io.IOException
      {
         x   = sc.nextInt();
         y   = sc.nextInt();
         lab = sc.next();
      }
   }

   // スタートからの経路を表すデータ
   // どういう文字列で来たか(spell)と、今どの節点にいるか(lastNode)
   static class Path implements Comparable<Path>
   {
      final String spell;
      final int lastNode;
      Path(String s, int n) { spell=s; lastNode=n; }

      // 辞書順で早いほうが「良い」経路
      public int compareTo( Path rhs ) {
         int c = spell.compareTo(rhs.spell);
         return c!=0 ? c : lastNode-rhs.lastNode;
      }
   }

   // 本体
   static String solve(int N, int S, int G, Arrow[] arr)
   {
      // とりあえず、ゴールに着けない節点に迷い込んではいけないので、
      // まず全節点について、ゴールまで行けるかどうか調べておく
      //
      // 方法は、幅優先探索(BFS)でも深さ優先探索(DFS)でも、
      // なんでもお好みの方法でいいと思います。

      boolean[][] reachable = new boolean[N][N];

      for(Arrow a : arr)
         reachable[a.x][a.y] = true;
      for(int k=0; k<N; ++k)
      {
         reachable[k][k] = true;
         for(int i=0; i<N; ++i)
            for(int j=0; j<N; ++j)
               reachable[i][j] |= reachable[i][k] & reachable[k][j];
      }

      // そもそもスタートからゴールに行けない場合は "NO"

      if( !reachable[S][G] )
         return "NO";

      // 行ける場合は、
      // スタートから始めて、「辞書順で早い順に」経路文字列を全探索していきます。

      TreeSet<Path> q = new TreeSet<Path>(); // 探索候補
      q.add( new Path("",S) ); // 「スタート S にいます」状態から探索開始

      for(;;)
      {
         Path p = q.pollFirst(); // 候補のなかから一番早いのを取り出す

         if( p.lastNode == G ) // ゴールについた
            return p.spell;

         if( p.spell.length() >= N*6 ) // ループしてる
            return "NO";

         for(Arrow a : arr) // それ以外:今の節点から一つパス延ばしたものを候補に入れる
            if( p.lastNode==a.x && reachable[a.y][G] )
               q.add( new Path(p.spell+a.lab, a.y) );
          // 本当は、"辞書順で最小の矢印と、それを接頭辞とするもの" だけ
          // 入れれば十分なんですが、コード書くの面倒なので全部入り…
      }
   }
}

複数の候補をためておいて、良い方から順に探索の舌を伸ばしていく、いわゆる

  最良優先探索

の一種になります。最良の候補を管理するデータ構造として、PriorityQueue ではなく TreeSet を使っているのは、こうしないと「経路は違うけど文字列は同じでlastNodeも同じ」なパスが大量にありえるため。TreeSet なら重複を自動的に除いてくれます。


「文字列の長さが 6N を超えたらループしてる」というのは、節点数 N で、ラベルの長さの最大が 6 なので、文字列の長さが 6N ということは、少なくともこの経路 p は同じ接点を2回以上通っているからです。

辞書順最小のループからは、すぐ抜け出すよりも1周回ってから抜け出す方が、それよりも2周回ってから抜け出す方が(以下略)お得なので、これが出てきた時点で "NO" が確定。


解法2:文字列より文字で

ちょっと考えてみると、解法1の面倒な部分は、ほとんどすべて「ラベルが文字列」だったことに起因していることがわかります。もし仮にラベルが全部一文字だったら、

3 3 0 2
0 1 a
0 2 ab
1 2 c

こんなデータで「うっかり接頭辞だから辞書順で早そうに見える方に突き進む」こともありません。全部1文字なので接頭辞だったりそうじゃなかったりしませんので。


つまり1文字なら TreeSet を用意して最良優先…のようなことを考えなくてよくなります。常に辞書順最小の矢印を進めばよい。ただし、最小ならどれも良いわけではなく

4 4 0 3
0 1 a
0 2 a
1 3 c
2 3 b

こんなデータで、0 から 1 と 2 どっちに行こうかなーと思って 1 だけに進んで 2 を忘れてしまうと "ab" が見つからずアウトです。ここは幅優先探索的に「最小矢印で行ける先は全部見る」とすればOKです。


というわけで、入力を無理矢理1文字ラベルに変換してから解く方法です。

/* C */
#include <stdio.h>
#include <memory.h>
typedef enum { false, true } bool;

enum { MAX_LAB  = 6 };
enum { MAX_EDGE = 400*MAX_LAB };
enum { MAX_NODE = 40+MAX_EDGE };
enum { MAX_ANS  = MAX_NODE };

/* 矢印データ */
typedef struct { int x, y; char lab; } arrow;

/* とりあえず、ゴールに着けない節点に迷い込んではいけないので、
 * まず全節点について、ゴールまで行けるかどうか調べておく */
void compute_reachabilty(arrow arr[], int N, int A, int G, bool r[])
{
   bool changed;

   r[G] = true;
   do {
      int i;
      changed = false;
      for(i=0; i<A; ++i)
         if( !r[arr[i].x] && r[arr[i].y] )
            r[arr[i].x] = changed = true;
   } while( changed );
}

void solve(arrow arr[], int N, int A, int S, int G, int SpellLimit)
{
   int  i,j,k;
   char answer[MAX_ANS+1] = {};
   bool cur[MAX_NODE] = {};
   bool r[MAX_NODE] = {};

   /* そもそもスタートからゴールに行けない場合は "NO" */
   compute_reachabilty(arr, N, A, G, r);
   if( !r[S] ) {
      puts("NO");
      return;
   }

   /* 「スタート S にいます」状態から探索開始 */
   cur[S] = true;

   for(i=0; i<SpellLimit; ++i) {
      if( cur[G] ) {
         /* ゴールについた */
         puts(answer);
         return;
      }

      /* 今の節点から辞書順最小の矢印で行けるところ全部に並列に進む */
      bool next[MAX_NODE] = {};
      char nextLab = 127;
      for(k=0; k<A; ++k)
         if( cur[arr[k].x] && r[arr[k].y] && arr[k].lab<=nextLab ) {
            if( arr[k].lab < nextLab ) {
               memset(next, 0, sizeof(next));
               nextLab = arr[k].lab;
            }
            next[arr[k].y] = true;
         }
      answer[i] = nextLab;
      memcpy(cur, next, sizeof(cur));
   }

   /* ループしてる */
   puts("NO");
}

int main()
{
   /* データセットを1個読み込んでは solve() を呼び出すループ */
   int N,A,S,G;
   while( scanf("%d%d%d%d",&N,&A,&S,&G), N|A|S|G ) {
      int i,j,ai=0,SpellLimit = MAX_LAB*N;
      arrow arr[MAX_EDGE];

      for(i=0; i<A; ++i) {
         int  fr, to;
         char lab[99];
         scanf("%d %d %s", &fr, &to, lab);

         /* 文字列ラベルを文字ラベルに分解。
          * つまり
          *    ・---abc--->・
          * のかわりに、新しく節点を二つ導入して
          *    ・---a--->・---b--->・---c--->・
          * という魔法の紋様と思い込むことにします。 */

         for(j=0; lab[j]; ++j) {
            int   x = j==0 ? fr : N-1;
            int   y = lab[j+1] ? N++ : to;
            arrow a = {x, y, lab[j]};
            arr[ai++] = a;
         }
      }
      solve(arr, N, ai, S, G, SpellLimit);
   }
   return 0;
}

ところで、ラベルを全部1文字に分解したものは、正規表現の実装に使われる 非決定性オートマトン によく似ています。実はそのものです。


上記のコードの cur や next のように「ある文字をたどって行ける先の節点を全部覚えておいて、全節点並列に1文字1文字進んでいく」方法で非決定性オートマトンを動かすやり方は、最近公開された Google 発の正規表現ライブラリ RE2 の中の人が、ちょっと昔に "Regular Expression Matching Can Be Simple And Fast" という記事で大プッシュしていました。読んでみると面白いかも。


解法3:最短経路の話

この問題はグラフの最短経路問題に見える…という見方もあります。


ふつう、最短経路というとグラフで

  「辺の重み⇔自然数や実数、最短⇔経路の重みの和が最小」

という場合ですが、それ以外の構造でも同じ問題を考えることができます。

  「辺の重み⇔文字列、最短⇔経路の重みの結合が辞書順最小」

とか。


大雑把に言うと、最短経路の有名アルゴリズムである ダイクストラ法ベルマン・フォード法

  min(x,y) + z = min(x+z, y+z)

を満たすような min 演算と + 演算が辺の重みに対して定義できれば、使えます。数値であることが本質ではなくて、最小値を選ぶ演算と、重みを結合する演算が存在して分配則を満たすことが、このあたりのアルゴリズムの本質です。日本語でいうと「v1→v2→...→vkの最短経路は、v1→...→v[k-1]の最短経路に辺を継ぎ足した形に必ずなってる」。ダイクストラは追加で

  min(x,x+z) = x

が必要だったかも。いわゆる「辺の重みが負でない」という。


と、えらそーに書いておきながら、文字列の辞書順の場合

  min(x,y) + z = min(x+z, y+z)

これ、成り立ちません。接頭辞の呪い。

  min("a","ab")+"c"="ac" ≠ "abc"=min("a"+"c","ab"+"c")

でも、逆から足すと

  z + min(x,y) = min(z+x, z+y)

これは成り立ちます。


というわけで、逆から、つまりスタートからゴールに向かってではなく、ゴールからスタートに向かって最短路を探すアルゴリズムが使えます。

// C++
#include <iostream>
#include <vector>
#include <string>
using namespace std;

struct arrow { int x, y; string lab; };
static const string INF = "NO";

string bwd_bellman_ford(const vector<arrow>& arr, int N, int A, int S, int G, int LOOP)
{
   vector<string> d(N, INF);
   d[G] = "";
   for(int t=0; t<N*LOOP; ++t)
   {
      vector<string> d2 = d;
      for(int i=0; i<A; ++i)
         if( d[arr[i].y]!=INF && (d2[arr[i].x]==INF || d2[arr[i].x]>arr[i].lab+d[arr[i].y]) )
         {
            d2[arr[i].x] = arr[i].lab + d[arr[i].y];
            if( arr[i].x==S && t>=N )
               return "NO"; // スタートからの最短経路にサイクルがある
         }
      d.swap(d2);
   }
   return d[S];
}

int main()
{
   for(int N,A,S,G,maxLab=0; cin>>N>>A>>S>>G, N|A|S|G;)
   {
      vector<arrow> arr(A);
      for(int i=0; i<A; ++i)
      {
         cin >> arr[i].x >> arr[i].y >> arr[i].lab;
         maxLab = max<int>(maxLab, arr[i].lab.size());
      }
      cout << bwd_bellman_ford(arr, N, A, S, G, 1+maxLab) << endl;
   }
}

結構これで挑んだチームは多いんじゃないかと思うんですが、ハマりどころがあるとすれば、

  • 最外周ループ for(int t=0; t<N*LOOP; ++t) のループ回数が足りない
40 41 0 39
0 1 x
1 1 x
1 39 y
0 2 xxxxxx
2 3 xxxxxx
(略)
37 38 xxxxxx
38 39 xxxxxz
    • 回れば回るほど辞書順でよくなるサイクル(最短経路のたとえで言えば"負の"サイクル)の検出は N 回まわせばいい、という知識を元に N 回で止めて "NO" を返してしまうと、「スタートからゴールへの辞書順最小経路とは関係ないところにあるサイクル」を検出して間違えて "NO" にしてしまう可能性があります。
    • 「負のサイクルの影響がスタートからの最短路まで波及した」ことを検知してはじめて "NO" を返す方が安全
      • N回目までにどっかでサイクルができて、N本矢印をたどればスタートにそのサイクルの影響が来るから 2N 回で安全!?
      • というのは間違いで、上の例のように、できた負のサイクルを何周かしないとスタートに影響できない場合があります。
      • これもループしない呪文の長さの上限 6N まで回してみればわかるので、そのくらいで打ち切ります。

  • 非負制約の「逆バージョン」 min(x,z+x)=x は文字列では成り立たない(min("z","a"+"z")="az"≠"z")ので、後ろからダイクストラは使えません。
    • ダイクストラするなら前からですが、接頭辞やサイクルに気をつける必要があります。
    • 気をつけると、要するに解法1になります。

  • 「前からベルマン・フォード」を殺す入力例としては
3 3 0 2
0 1 a
0 1 ab
1 2 c
    • こんなものがあります。
    • 0 から 1 への最短路は常に "a" なので、何周更新しても 1 への最短路は "a" のまま変わりません。ですが、0 から 1 を経由して 2 に行く最短路の時には "ab" を経由する必要があって、困ります。
    • 「最短路と、それを接頭辞として持つもの全部」を記憶しておく形に変えた「前からベルマン・フォード」なら上手く行くかも。

cafeliercafelier2010/07/06 14:53解法4:長さ L で行ける辞書順最小路、を各長さ L (≦6N)各ノードについてもとめる動的計画法(あるいは拡張ダイクストラ)
というのがあるらしい。
http://d.hatena.ne.jp/Tayama/20100703/1278141167
http://d.hatena.ne.jp/kusano_prog/20100705/1278361318
http://d.hatena.ne.jp/pes_magic/20100706/1278367137
なるほどなー。

トラックバック - http://topcoder.g.hatena.ne.jp/cafelier/20100703

presented by cafelier/k.inaba under CC0