cafelier@SRM

SRM425 Div1 1000 (2008/12/7)

あとで | 17:43 |

425の1000点問題を考える。

一昨日くらいにようやくd:id:tsukunoの言ってたDPというのを理解したので、あんど、O(3^(N-1) * N) にはなるんじゃないだろうかと思ったので組んでみてました。最悪ケース（16頂点、確率0%や100%の辺が一個もない）で3秒ちょい。少し足りない…

#include <vector>
#include <string>
#include <map>
using namespace std;
typedef unsigned int uint;

static const uint OUT=0, IN=1, ROOT=1, LEAF=2, HAVELEAF=0xAAAAAAAA;
#define is(mm,j,t)    (((mm>>(j*2))&3)==(t))
#define isnot(mm,j,t) (((mm>>(j*2))&3)!=(t))

class RoadsOfKingdom
{
public:
  map<uint,double> memo;
  double P[16][16];

  double getProbability(vector <string> roads) 
  {
    // 入力
    int N = roads.size();
    for(int i=0; i<N; ++i)
      for(int j=0; j<N; ++j)
        P[i][j] = (roads[i][j]-'0')/8.0;

    // メモ化再帰
    double ans = 0.0;
    for(uint m=0; m<(1<<N); m+=2) { // すべての全域木mmに関して...
      uint mm = ROOT;
      for(int i=1; i<N; ++i)
        mm |= ((m>>i)&1 ? IN : LEAF) << i*2;
      ans += rec(mm, N); // その全域木ができる確率
    }
    return ans;
  }

  // 部分木mmを、元グラフの各ノードを
  //  (0: mmに入ってない, 1: mmの内部ノード, 2: mmのリーフ)
  // というビットパターンで表す。16頂点なので32bitで足りる。
  double rec(uint mm, int N)
  {
    // trivial cases
    if( mm == ROOT )
      return 1.0;
    if( isnot(mm,0,ROOT) || !(mm&HAVELEAF) )
      return 0.0;

    // memoized?
    map<uint,double>::iterator it = memo.find(mm);
    if( it != memo.end() )
      return it->second;

    // 一番頂点番号の小さいリーフ f をとりのぞく
    int f = 0;
    for(; f<N; ++f)
      if( is(mm,f,LEAF) )
        break;

    // P[f][*] が全部ぶっ壊れる確率
    double alldead = 1.0;
    int numAlive = 0;
    for(int j=0; j<N; ++j)
      if( isnot(mm,j,OUT) )
        if( P[f][j] == 1.0 )
          ++numAlive;
        else
          alldead *= 1.0 - P[f][j];
    if( numAlive >= 2 )
      return 0.0;

    // 「fの親がjだった場合」をmmの全ての内部ノードjについて調べて足し算
    double px = rec( mm^(LEAF<<f*2), N ); // jがf以外にも子供を持ってる場合の確率

    double ans = 0.0;
    for(int j=0; j<N; ++j)
      if( is(mm,j,IN) && P[f][j] != 0.0 )
      {
        double pdead = alldead;
        if( P[f][j] != 1.0 )
          if( numAlive ) continue;
          else pdead /= 1.0 - P[f][j];
        double prest = rec( mm^(LEAF<<f*2)^(3u<<j*2), N ) + px; // fがjの唯一の子の場合
        ans += P[f][j] * pdead * prest;
      }
    return memo[mm] = ans;
  }
};

1.5倍くらいだと小手先の技でもどうにかなりそうな気がしなくもないのですが、まあそれよりはアルゴリズムの改良でなんとかしたい…。そうこうしているうちにEditorialが出てしまったけれど見ないぞー。

実際は 3^(N-1) もなくて、memoに入る状態数は高々120万弱なんですよね。常に最初のLeafを除く方法で出てくる木構造はそのくらいしかない。なのでそういう状態だけもっと直接的に列挙する方法とかがあれば良さそうなんだけど。うむむ