土下座しながら探索中

主に競技プログラミング

AOJ 2437 : DNA

問題リンク : http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2437&lang=jp

問題概要 :
日本語なので略

解法:
DP + 半分全列挙で解いた

まず、i文字目(1<=i<=Na+Nt+Ng+Nc)にくる文字が何かを計算する
与えられる文法に含まれる非終端記号間の関係をグラフにするとDAGなので、これは再帰などでその通り計算すると良い
が計算途中で文字数がNa+Nt+Ng+Ncを超えるようなら、答えは0なので枝刈すること ( 出ないとTLEする )
文字数がNa+Nt+Ng+Ncを超えるようなら、枝刈するので、計算にはそのくらいの時間しかかからない

次に、区間[0,(Na+Nt+Ng+Nc)/2)と[(Na+Nt+Ng+Nc)/2,(Na+Nt+Ng+Nc))のそれぞれについて、
その区間で作れる文字列とそれが何通りあるかを計算する
これは動的計画法で求めると良い dp[i文字目][Aの数][Tの数][Gの数][Cの数] := A, T, G, Cがそれぞれ(Aの数), (Tの数), (Gの数), (Cの数)だけ出現する文字列が何通りあるか
(i文字目)はその最大値分だけ配列を用意するとMLEするので、2にして使い回すこと
O((Na+Nt+Ng+Nc) * Na * Nt * Ng * Nc)
半分に分けたので2.5 * 10^9くらい
AOJなら3 * 10^9は一瞬なのでOK

最後に、2つの区間について半分全列挙で答えを求める
O(Na * Nt * Ng * Nc)

コード (デバッグコード付き debug = 1にするとデバッグモード) :

#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;++i)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-8)
#define equals(a,b) (fabs((a)-(b))<EPS)
#define pb push_back
#define ALL(x) x.begin(),x.end()

using namespace std;

typedef long long ll;

bool LT(double a,double b) { return !equals(a,b) && a < b; }
bool LTE(double a,double b) { return equals(a,b) || a < b; }

#define A_ 0
#define T_ 1
#define G_ 2
#define C_ 3
#define NONTERMINAL 0
#define TERMINAL 1

struct Data {
  int type; // NONTERMINAL or TERMINAL
  int dst; // for NONTERMINAL
  int bit; // for TERMINAL
};

#define MAX 60
vector<int> N;
int m;
vector<Data> G[MAX];

const ll mod = 1000000007LL;

bool debug = 0;

inline int toIDX(char c) { if( c == 'A' ) return A_; if( c == 'T' ) return T_; if( c == 'G' ) return G_; if( c == 'C' ) return C_; assert(false); }
inline char toCHR(int  i) { if( i == A_ ) return 'A'; if( i == T_ ) return 'T'; if( i == G_ ) return 'G'; if( i == C_ ) return 'C'; assert(false); }


bool fail;
vector<int> chars;
int enum_char_dfs(int cur) {
  if( fail ) return -1;
  int ret = 0;
  rep(i,(int)G[cur].size()) {
    Data &d = G[cur][i];
    if( d.type == NONTERMINAL ) {
      ret += enum_char_dfs(d.dst);
      if( ret >= N[0] + N[1] + N[2] + N[3] + 1 ) { fail = true; return -1; }
    } else {
      chars.pb(d.bit);
      //++len;
      ++ret;
    }
  }
  return ret;
}

int dp[2][52][52][52][52];
int solve_DP(int L,int R) { //[L,R)
  memset(dp,0,sizeof dp);
  dp[0][0][0][0][0] = 1;
  int idx = 0;
  REP(i,L,R) {
    rep(w,N[0]+1) rep(x,N[1]+1) rep(y,N[2]+1) rep(z,N[3]+1) dp[(idx+1)&1][w][x][y][z] = 0;
    rep(a,N[0]+1) {
      rep(t,N[1]+1) {
	rep(g,N[2]+1) {
	  rep(c,N[3]+1) {
	    int &score = dp[idx&1][a][t][g][c];
	    if( score == 0 ) continue;
	    if( (chars[i]>>A_) & 1 ) ( dp[(idx+1)&1][a+1][t][g][c] += score ) %= mod;
	    if( (chars[i]>>T_) & 1 ) ( dp[(idx+1)&1][a][t+1][g][c] += score ) %= mod;
	    if( (chars[i]>>G_) & 1 ) ( dp[(idx+1)&1][a][t][g+1][c] += score ) %= mod;
	    if( (chars[i]>>C_) & 1 ) ( dp[(idx+1)&1][a][t][g][c+1] += score ) %= mod;
	  }
	}
      }
    }
    ++idx;
  }
  return idx;
}

const ll ccoef = 1000000;
const ll gcoef = 10000;
const ll tcoef = 100;
inline ll toHash(int a,int t,int g,int c) { return a + t * tcoef + g * gcoef + c * ccoef; }

void toMap(int idx, map<ll,ll> &mp) {
  rep(a,N[0]+1) {
    rep(t,N[1]+1) {
      rep(g,N[2]+1) {
	rep(c,N[3]+1) {
	  ll hs = toHash(a,t,g,c);
	  ll score = dp[idx][a][t][g][c];
	  if( score == 0 ) continue;
	  mp[hs] = score;
	}
      }
    }
  }
}

void solve() {
  if( debug ) {
    cout << "* Graph ---" << endl;
    rep(i,m) {
      cout << "  - i = " << i << "-th:" << endl;
      rep(j,(int)G[i].size()) {
	Data &d = G[i][j];
	if( d.type == NONTERMINAL ) {
	  cout << "  [NONTERMINAL] : to " << d.dst << endl;
	} else {
	  assert( d.bit != -1 );
	  bitset<4> BIT(d.bit);
	  cout << "  [   TERMINAL] : bit = " << BIT << endl;
	}
      }
      puts("");
    }
  }
  
  // enumerate chars
  chars.clear();
  fail = false;
  int len = enum_char_dfs(0);
  if( fail ) { puts("0"); return ; }
  if( debug ) {
    assert( len == (int)chars.size() );
    cout << "* enumerate chars ---" << endl;
    cout << "  - len = " << len << endl;
    cout << "    ";
    rep(i,len) {
      cout << "[";
      rep(j,4) {
	if( (chars[i]>>j) & 1 ) {
	  cout << toCHR(j);
	}
      }
      cout << "] ";
    }
    puts("");
  }

  // Dynamic Programming!!
  if( len > N[0] + N[1] + N[2] + N[3] ) { puts("0"); return ; }
  int sz = len / 2;
  if( debug ) {
    cout << "[SPLIT] -- [" << 0 << "," << sz << ") and [" << sz << "," << len << ")" << endl;
  }
  map<ll,ll> mp[2];
  int lst;
  lst = solve_DP(0,sz);
  toMap(lst&1,mp[0]);



  if( debug ) {
    cout << "* mp[0] Info ---" << endl;
    for(auto v : mp[0]) {
      ll hs = v.first;
      int c = hs / ccoef;
      hs -= ( c * ccoef );
      int g = hs / gcoef;
      hs -= ( g * gcoef );
      int t = hs / tcoef;
      hs -= ( t * tcoef );
      int a = hs;
      ll score = v.second;
      cout << "  (A,T,G,C) = (" << a << "," << t << "," << g << "," << c << ")" << endl;
      cout << "      score = " << score << endl;
    }
  }

  lst = solve_DP(sz,len);
  toMap(lst&1,mp[1]);
  if( debug ) {
    cout << "* mp[1] Info ---" << endl;
    for(auto v : mp[1]) {
      ll hs = v.first;
      int c = hs / ccoef;
      hs -= ( c * ccoef );
      int g = hs / gcoef;
      hs -= ( g * gcoef );
      int t = hs / tcoef;
      hs -= ( t * tcoef );
      int a = hs;
      ll score = v.second;
      cout << "  (A,T,G,C) = (" << a << "," << t << "," << g << "," << c << ")" << endl;
      cout << "      score = " << score << endl;
    }
  }

  // hanbun-zennrekkyo
  ll answer = 0;
  for(auto v : mp[0]) {
    ll hs = v.first;
    int c = hs / ccoef;
    hs -= ( c * ccoef );
    int g = hs / gcoef;
    hs -= ( g * gcoef );
    int t = hs / tcoef;
    hs -= ( t * tcoef );
    int a = hs;
    ll score = v.second;
    if( a > N[0] || t > N[1] || g > N[2] || c > N[3] ) continue;
    int rem_a = N[0] - a;
    int rem_t = N[1] - t;
    int rem_g = N[2] - g;
    int rem_c = N[3] - c;
    ll next_hs = toHash(rem_a,rem_t,rem_g,rem_c);
    if( mp[1].count(next_hs) ) {
      ( answer += ( ( score * mp[1][next_hs] ) % mod ) ) %= mod;
    }
  }
  cout << answer << endl;
}

void parse(string s, map<string,int> &mp, int &id) {
  rep(i,(int)s.size()) if( s[i] == ':' ) { s[i] = ' '; break; }
  stringstream ss;
  ss << s;
  vector<string> vec;
  while( ss >> s ) vec.pb(s);
  if( !mp.count(vec[0]) ) { mp[vec[0]] = id++; }
  int cur = mp[vec[0]];
  REP(i,1,(int)vec.size()) {
    string &t = vec[i];
    if( t[0] == '[' ) {
      int ptr = 1, bit = 0;
      while( ptr < (int)t.size() && t[ptr] != ']' ) bit = bit | (1<<toIDX(t[ptr++]));
      G[cur].pb((Data){TERMINAL,-1,bit});
    } else {
      if( !mp.count(vec[i]) ) { mp[vec[i]] = id++; }
      int nex = mp[vec[i]];
      G[cur].pb((Data){NONTERMINAL,nex,-1});
    }
  }
}

int main() {
  N.resize(4);
  rep(i,4) cin >> N[i];
  cin >> m;
  cin.ignore();
  string line;
  map<string,int> mp;
  int id=0;
  rep(i,m) {
    getline(cin,line);
    parse(line,mp,id);
  }
  solve();
  return 0;
}

LiveArchive 6396 : Factors

問題リンク : https://icpcarchive.ecs.baylor.edu/external/63/6396.pdf

問題概要 :
関数 f は2^63未満の正の整数 k を引数にとり、
k の素因数の異なるarrangementの数 n を返す関数である
例えば、 f(10) = 2 である
なぜなら、 10 = 2 * 5 = 5 * 2 であるから
また、 f(20) = 3 である
なぜなら、 20 = 2 * 2 * 5 = 2 * 5 * 2 = 5 * 2 * 2 であるから

2^63未満の正の整数 n が与えられるので、 f(k) = n となるような最小の k を求めよ
入力で与えられる n は k が 2^63未満となるものに限られる

解法 :
解説読みました..
http://www.csc.kth.se/~austrin/icpc/finals2013solutions.pdf

k = p1^e1 * p2^e2 * ... * pt^et とすると、 n = ( e1 + e2 + ... + et )! / ( e1! * e2! * ... * et! )
見ての通り、 n はp1, p2, ..., pt の値に影響しない
ならば、kを最小化したいなら p1, p2, ..., pt は小さい素数から順に割り当ててしまうのが良い
( p1 = 2, p2 = 3, ... )
また、 e1 >= e2 >= ... >= et としてもよい
ei < ej ( i >= j ) であるようなものがあっても良いのだが、
ei ejの順序はnに影響しないので、そのような場合はeiとejをスワップしてしまったほうが k を小さくできる
まとめると、
k = p1^e1 * p2^e2 * ... * pt^et ( p1 = 2, p2 = 3, ... かつ e1 >= e2 >= ... かつ k < 2^63 )
さて、このk、ここまで制約を厳しくすると、意外とそこまで多くないような気がする
使う素数は高々15個くらいまで ( それ以上だと、eiを全て1にしても2^63以上になる
また、e1から順番に全探索で決め打ちしていくとしても、eiの値は前よりも減少していくわけだし、kが2^63以上ならそこで終わりなので、
ちょっと全探索を書いて終わりそうか見てみると爆速なので、メモ化して投げるとAC
( e1は高々63で、e2は高々その半分もなくて、e3はe2のその半分もなくて、...と考えていくと、かなり大雑把に見積もっても2 * 10^6以下になるので、全探索で大丈夫そう )
(ちなみに解説によると、 上記の制約を満たすkは43606しかないらしい)


ちなみに、long longや__int128を使っても、素直に ( e1 + e2 + ... + et )! / ( e1! * e2! * ... * et! )を計算するとオーバーフローする
多倍長を書くか、JavaのBigIntegerを使うか、割り算をいい感じにするかしないとダメ

コード :

#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;++i)
#define rep(i,n) REP(i,0,n)
#define ALL(x) x.begin(),x.end()

using namespace std;

std::ostream &operator<<(std::ostream &dest, __int128_t value) {
  std::ostream::sentry s(dest);
  if (s) {
    __uint128_t tmp = value < 0 ? -value : value;
    char buffer[128];
    char *d = std::end(buffer);
    do {
      --d;
      *d = "0123456789"[tmp % 10];
      tmp /= 10;
    } while (tmp != 0);
    if (value < 0) {
      --d;
      *d = '-';
    }
    int len = std::end(buffer) - d;
    if (dest.rdbuf()->sputn(d, len) != len) {
      dest.setstate(std::ios_base::badbit);
    }
  }
  return dest;
}

const int m = 16;
__int128 maxi;
__int128 prim[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53};

map<__int128,__int128> rf;
__int128 calc_k(deque<int> &inp) {
  deque<int> deq;
  rep(i,(int)inp.size()) REP(j,2,inp[i]+1) deq.push_back(j);
  int limit = 0;
  rep(i,(int)inp.size()) limit += inp[i];
  sort(ALL(deq));
  __int128 k = 1;
  int i = 2;
  while( i <= limit ) {
    k *= (__int128)i;
    while( !deq.empty() && ( k % deq.front() == 0 ) ) {
      k /= deq.front();
      deq.pop_front();
    }
    while( !deq.empty() && ( k % deq.back() == 0 ) ) {
      k /= deq.back();
      deq.pop_back();
    }
    ++i;
  }
  while( !deq.empty() && ( k % deq.front() == 0 ) ) {
    k /= deq.front();
    deq.pop_front();
  }
  assert(deq.empty());
  return k;
}

void dfs(int cur, __int128 v,int limit,deque<int> &deq) {
  if( cur >= m ) return;
  __int128 w = v;
  rep(i,(limit+1)) {
    w = w * prim[cur];
    if( w >= maxi ) break;
    deq.push_back(i+1);
    __int128 k = calc_k(deq);
    assert( k );
    if( !rf.count(k) ) rf[k] = w;
    else rf[k] = min(rf[k],w);
    dfs(cur+1,w,i+1,deq);
    deq.pop_back();
  }
}

__int128 parse(string s) {
  __int128 n = 0;
  rep(i,(int)s.size()) {
    n *= (__int128)10;
    n += (__int128)( s[i] - '0' );
  }
  return n;
}

int main() {
  maxi = 1;
  REP(i,1,64) maxi = maxi * (__int128)2;
  deque<int> deq;
  dfs(0,1,63,deq);
  string s;
  while( cin >> s ) {
    __int128 n = parse(s);
    cout << n << " " << rf[n] << endl;
  }
  return 0;
}

SRM 730 easy : StonesOnATree

問題概要:
各頂点に重みのついた木が与えられる

.. 木の情報は以下の通り
.... 木の頂点数はnで表される
.... 木の各頂点には0からn-1までの番号が割り振られており、
.... i番目(0<=i<=n-1)の頂点の重みはw[i]で表される
.... また、各頂点の子の数は高々2である
この木を使ってゲームをする
ゲームで行える操作は次の2つ
1. どれか1つ石の置かれていない頂点を選び、石を置く
..ただし、石を置くには以下の制約を満たさなければならない
....石を置きたい頂点の全ての子供に石が置かれている
..石を置くと、重みにコストに足す
2. どれか1つ石の置かれている頂点を選び、その頂点から石を取り除く
..石を取り除くと、その頂点に石を置いたときに加算した重みをコストから引く
コストの初期値は0である
木の根に石を置くまでにかかったコストの最大値の最小値を求めよ

解説:
葉の部分から数えていく これは冗談でございます

解法の基本的な方針は次の通り
ある頂点に石を置いたら、その頂点の子の石は即座に取り除く、といった感じで葉から順に計算していく


もうすこししっかり書くと、
答えを求める関数をdfsとする
dfsは0以上n-1以下の整数iを受け取り、
そのiを根とする部分木について概要のゲームをといたときの答えを返す
このとき、dfsは次のように定義できる

BASIS :
.. Case : 頂点が葉のとき
頂点の番号をiとしたとき、答えは明らかにw[i]なのでw[i]を返す

INDUCTION:
.. Case : 頂点が1つしか子を持たないとき
今見ている頂点の番号をi,子の頂点の番号をjとすると
答えとなりうる値はw[i]+w[j]か頂点iに石を置くために必要なコストの最大値のどちらか
選択の余地はないので大きい方をかえす

.. Case : 頂点が2つ子を持つとき
答えとなりうる最大値は以下の4つのいずれか
左の子→右の子→今見ている頂点の順番で石を置いたときの最大値か
右の子→左の子→今見ている頂点の順番で石を置いたときの最大値か
左の子に石を置くために必要なコストの最大値か
右の子に石を置くために必要なコストの最大値か
これは選択できるので、もっとも小さい値を返す

コード:

const int IINF = INT_MAX;
#define MAX 20000
vector<int> G[MAX];
vector<int> w;
int V;

int dfs(int cur) {
  if( G[cur].size() == 0 ) {
    return w[cur];
  } if( G[cur].size() == 1 ) {
    return max(w[cur]+w[G[cur][0]],dfs(G[cur][0]));
  }
  assert( G[cur].size() == 2 );
  int ret = 0;
  int L = G[cur][0], R = G[cur][1];
  int Lmaxi = dfs(L), Rmaxi = dfs(R);
  ret = max(Lmaxi,Rmaxi);
  ret = max(ret,min(Lmaxi+w[R],Rmaxi+w[L]));
  ret = max(ret,w[cur]+w[L]+w[R]);
  return ret;
}

class StonesOnATree {
public:
  int minStones(vector <int> p, vector <int> _w) {
    V = _w.size();
    rep(i,V) w.push_back(_w[i]);
    rep(i,(int)p.size()) G[p[i]].push_back(i+1);
    return dfs(0);
  }
};

AOJ 2099 : Walk under a Scorching Sun

問題リンク : Walk under a Scorching Sun | Aizu Online Judge

問題概要 :
いくつかのの道と建物に加え太陽がある。
道、建物、太陽はそれぞれ線分、凸多角形+高さ、点として表される。

道の上にあるある点xと太陽を結ぶ直線が建物と交差するなら、
その点xは日陰である。
それ以外の点は日向である。

道の上にある点sから別の道の上にある点tへ移動するときに通らなければならない日向の距離の総和の最小値を求めよ。

解法:
1. 建物からできる影を計算
2. 道を表す線分と影の交点を列挙
3. 線分アレンジメント
4. アレンジメントでできたグラフの辺のうち、影にあたる辺の重みを0にする
5. dijkstraでsからtまでの最短距離を求める

以下備忘録
→ 影を計算する
f:id:tei3344:20170704004837p:plain

ちなみに、問題文中に書いてあったかどうかわからないが、
道と建物が接点をもつことはないらしい

コード :

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cassert>
#include<cmath>
#include<complex>
#include<cstring>
#include<climits>
#include<vector>
#include<algorithm>
#include<queue>
#include<deque>
#include<bitset>
#include<sstream>
#include<map>
#include<set>

#define REP(i,s,n) for(int i=s;i<n;++i)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-7)
#define equals(a,b) (fabs((a)-(b))<EPS)
#define COUNTER_CLOCKWISE 1
#define CLOCKWISE -1 
#define ONLINE_BACK 2
#define ONLINE_FRONT -2
#define ON_SEGMENT 0


using namespace std;

const bool debug = false; // DEBUG - DEBUG - DEBUG - DEBUG - DEBUG - DEBUG - DEBUG - DEBUG - DEBUG
 
inline double toRad(double theta){
  return theta * M_PI / 180.0;
}

// x >= y
bool GTE(double x, double y){ return equals(x,y) || x > y; }
 
// x > y
bool GT(double x, double y){ return !equals(x,y) && x > y; }
 
// x < y
bool LT(double x, double y){ return !equals(x,y) && x < y; }
 
// x <= y
bool LTE(double x, double y){ return equals(x,y) || x < y; }

// -- Library 2d -- BEGIN ----------------------------------------
class Point{
public:
  double x,y;
 
  Point(double x = 0,double y = 0): x(x),y(y){}
 
  Point operator + (Point p){return Point(x+p.x,y+p.y);}
  Point operator - (Point p){return Point(x-p.x,y-p.y);}
  Point operator * (double a){return Point(a*x,a*y);}
  Point operator / (double a){return Point(x/a,y/a);}
  Point operator * (Point a){ return Point(x * a.x - y * a.y, x * a.y + y * a.x); }
 
  bool operator < (const Point& p) const{ return !equals(x,p.x)?x<p.x:(!equals(y,p.y)&&y<p.y); }
 
  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }
 
};
 
struct Segment{
  Point p1,p2;
  Segment(Point p1 = Point(),Point p2 = Point()):p1(p1),p2(p2){}
  bool operator <  (const Segment& s) const { return ( p2 == s.p2 ) ? p1 < s.p1 : p2 < s.p2; }
  bool operator == (const Segment& s) const { return ( s.p1 == p1 && s.p2 == p2 ) || ( s.p1 == p2 && s.p2 == p1 ); }
 
};
 
typedef Point Vector;
typedef Segment Line;
typedef vector<Point> Polygon;
 
ostream& operator << (ostream& os,const Point& a){ return os << "(" << a.x << "," << a.y << ")"; }
 
ostream& operator << (ostream& os,const Segment& a){ return os << "( " << a.p1 << " , " << a.p2 << " )"; }
 
double dot(Point a,Point b){ return a.x*b.x + a.y*b.y; }
 
double cross(Point a,Point b){ return a.x*b.y - a.y*b.x; }
 
double norm(Point a){ return a.x*a.x+a.y*a.y; }
 
double abs(Point a){ return sqrt(norm(a)); }
 
//rad は角度をラジアンで持たせること
Point rotate(Point a,double rad){ return Point(cos(rad)*a.x - sin(rad)*a.y,sin(rad)*a.x + cos(rad)*a.y); }
 
// a => prev, b => cur, c=> next
// prev から cur へ行って next へ行く際の角度を求める
double getArg(Point a,Point b,Point c){
  double arg1 = atan2(b.y-a.y,b.x-a.x);
  double arg2 = atan2(c.y-b.y,c.x-b.x);
  double arg = fabs( arg1 - arg2 );
  while( arg > M_PI ) arg -= 2.0 * M_PI;
  return fabs(arg);
}

int ccw(Point p0,Point p1,Point p2){
  Point a = p1-p0;
  Point b = p2-p0;
  if(cross(a,b) > EPS)return COUNTER_CLOCKWISE;
  if(cross(a,b) < -EPS)return CLOCKWISE;
  if(dot(a,b) < -EPS)return ONLINE_BACK;
  if(norm(a) < norm(b))return ONLINE_FRONT;
  return ON_SEGMENT;
}

bool intersectLL(Line l, Line m) {
  return abs(cross(l.p2-l.p1, m.p2-m.p1)) > EPS || // non-parallel
         abs(cross(l.p2-l.p1, m.p1-l.p1)) < EPS;   // same line
}
bool intersectLS(Line l, Line s) {
  return cross(l.p2-l.p1, s.p1-l.p1)*       // s[0] is left of l
         cross(l.p2-l.p1, s.p2-l.p1) < EPS; // s[1] is right of l
}
bool intersectLP(Line l,Point p) {
  return abs(cross(l.p2-p, l.p1-p)) < EPS;
}
bool intersectSS(Line s, Line t) {
  return ccw(s.p1,s.p2,t.p1)*ccw(s.p1,s.p2,t.p2) <= 0 &&
         ccw(t.p1,t.p2,s.p1)*ccw(t.p1,t.p2,s.p2) <= 0;
}
bool intersectSP(Line s, Point p) {
  return abs(s.p1-p)+abs(s.p2-p)-abs(s.p2-s.p1) < EPS; // triangle inequality
}
 
Point projection(Line l,Point p) {
  double t = dot(p-l.p1, l.p1-l.p2) / norm(l.p1-l.p2);
  return l.p1 + (l.p1-l.p2)*t;
}
Point reflection(Line l,Point p) {
  return p + (projection(l, p) - p) * 2;
}
double distanceLP(Line l, Point p) {
  return abs(p - projection(l, p));
}
double distanceLL(Line l, Line m) {
  return intersectLL(l, m) ? 0 : distanceLP(l, m.p1);
}
 
double distanceLS(Line l, Line s) {
  if (intersectLS(l, s)) return 0;
  return min(distanceLP(l, s.p1), distanceLP(l, s.p2));
}
double distanceSP(Line s, Point p) {
  Point r = projection(s, p);
  if (intersectSP(s, r)) return abs(r - p);
  return min(abs(s.p1 - p), abs(s.p2 - p));
}
 
double distanceSS(Line s, Line t) {
  if (intersectSS(s, t)) return 0;
  return min(min(distanceSP(s, t.p1), distanceSP(s, t.p2)),
             min(distanceSP(t, s.p1), distanceSP(t, s.p2)));
}
 
Point crosspoint(Line l,Line m){
  double A = cross(l.p2-l.p1,m.p2-m.p1);
  double B = cross(l.p2-l.p1,l.p2-m.p1);
  if(abs(A) < EPS && abs(B) < EPS){
    vector<Point> vec;
    vec.push_back(l.p1),vec.push_back(l.p2),vec.push_back(m.p1),vec.push_back(m.p2);
    sort(vec.begin(),vec.end());
    assert(vec[1] == vec[2]); //同じセグメントかもよ
    return vec[1];
    //return m.p1;
  }
  if(abs(A) < EPS)assert(false);
  return m.p1 + (m.p2-m.p1)*(B/A);
}

struct Edge {
  int from,to;
  double cost;
  Edge(int from=0,int to=0,double cost=0):from(from),to(to),cost(cost){}
  bool operator < (const Edge& a)const { return !equals(cost,a.cost) && cost < a.cost; }
};
 
vector<vector<Edge> > segmentArrangement(vector<Segment> vs,vector<Point> &ps) {
/*
端点もいれたいなら
ここであらかじめ端点だけpsにいれる
*/
  rep(i,(int)vs.size())
    REP(j,i+1,(int)vs.size())
    if(intersectSS(vs[i],vs[j]))
      ps.push_back(Point(crosspoint(vs[i],vs[j])));
 
  sort(ps.begin(),ps.end());  
  ps.erase(unique(ps.begin(),ps.end()),ps.end());
 
  vector<vector<Edge> > ret(ps.size());
 
  for(int i=0;i<(int)vs.size();i++){
    vector<pair<double,int> > list;
    rep(j,(int)ps.size()) if(intersectSP(vs[i],ps[j]))list.push_back(pair<double,int>(norm(vs[i].p1-ps[j]),j));
    sort(list.begin(),list.end());
    for(int j=0;j+1<(int)list.size();++j) {
      int from = list[j].second, to = list[j+1].second;
      double cost = abs(ps[from]-ps[to]);
      ret[from].push_back(Edge(from,to,cost));
      ret[to].push_back(Edge(to,from,cost));
    }
  }  
  return ret;
}
 
//cross product of pq and pr
double cross3p(Point p,Point q,Point r) { return (r.x-q.x) * (p.y -q.y) - (r.y - q.y) * (p.x - q.x); }
  
//returns true if point r is on the same line as the line pq
bool collinear(Point p,Point q,Point r) { return fabs(cross3p(p,q,r)) < EPS; }
  
//returns true if point t is on the left side of line pq
bool ccwtest(Point p,Point q,Point r){
  return cross3p(p,q,r) > 0; //can be modified to accept collinear points
}
 
bool onSegment(Point p,Point q,Point r){
  return collinear(p,q,r) && equals(sqrt(pow(p.x-r.x,2)+pow(p.y-r.y,2)) + sqrt(pow(r.x-q.x,2) + pow(r.y-q.y,2) ),sqrt(pow(p.x-q.x,2)+pow(p.y-q.y,2)) ) ;
}
  
 
bool isConvex(vector<Point> p) {
  int sz = (int)p.size();
  if(sz < 3)return false;//boundary case, we treat a point or a line as not convex
  bool isLeft = ccwtest(p[0],p[1],p[2]);
  for(int i=1; i<(int)p.size();i++)
    if(ccwtest(p[i],p[(i+1)%sz],p[(i+2)%sz]) != isLeft) return false;
  return true;
}
 
 
double angle(Point a,Point b,Point c) {
  double ux = b.x - a.x, uy = b.y - a.y;
  double vx = c.x - a.x, vy = c.y - a.y;
  return acos((ux*vx + uy*vy)/sqrt((ux*ux + uy*uy) * (vx*vx + vy*vy)));
}  
  
//多角形poly内(線分上も含む)に点pが存在するかどうは判定する  
bool inPolygon(Polygon poly,Point p){
  if((int)poly.size() == 0)return false;
  rep(i,(int)poly.size()) if(onSegment(poly[i],poly[(i+1)%poly.size()],p))return true;
  double sum = 0;
  for(int i=0; i < (int)poly.size() ;i++) {
    if( equals(cross(poly[i]-p,poly[(i+1)%poly.size()]-p),0.0) ) continue; // 3点が平行だとangleがnanを返しsumがnanになり死ぬ
    if( cross3p(p,poly[i],poly[(i+1)%poly.size()]) < 0 ) sum -= angle(p,poly[i],poly[(i+1)%poly.size()]);
    else                                                 sum += angle(p,poly[i],poly[(i+1)%poly.size()]);
  }
  // あまり誤差を厳しくしすぎると良くないので以下のほうが良い 
  const double eps = 1e-5;
  return (fabs(sum - 2*M_PI) < eps || fabs(sum + 2*M_PI) < eps);
}  


bool inPolygon(const Polygon &poly,Segment seg){
  vector<Point> vp;
  vp.push_back(seg.p1);
  vp.push_back(seg.p2);
  rep(i,(int)poly.size()) {
    Segment edge = Segment(poly[i],poly[(i+1)%(int)poly.size()]);
    if( equals(cross(seg.p1-seg.p2,edge.p1-edge.p2),0.0) ) {
      if( onSegment(seg.p1,seg.p2,edge.p1) ) vp.push_back(edge.p1);
      if( onSegment(seg.p1,seg.p2,edge.p2) ) vp.push_back(edge.p2);
    } else {
      if( intersectSS(seg,edge) ) vp.push_back(crosspoint(seg,edge));
    }
  }
  sort(vp.begin(),vp.end());
  vp.erase(unique(vp.begin(),vp.end()),vp.end());
  rep(i,(int)vp.size()-1) {
    Point middle_point = ( vp[i] + vp[i+1] ) / 2.0;
    if( !inPolygon(poly,middle_point) ) return false;
  }
  return true;  
}
 
Polygon andrewScan(Polygon s)
{
  Polygon u,l;
  if((int)s.size() < 3)return s;
 
  sort(s.begin(),s.end());
 
  u.push_back(s[0]);
  u.push_back(s[1]);
  l.push_back(s[s.size()-1]);
  l.push_back(s[s.size()-2]);
 
  for(int i=2;i<(int)s.size();i++)
    {
      for(int n=(int)u.size();n >= 2 && ccw(u[n-2],u[n-1],s[i]) != CLOCKWISE; n--)
	u.pop_back();
      u.push_back(s[i]);
    }
 
  for(int i=(int)s.size()-3; i>=0 ; i--)
    {
      for(int n=(int)l.size(); n >= 2 && ccw(l[n-2],l[n-1],s[i]) != CLOCKWISE; n--)
	l.pop_back();
      l.push_back(s[i]);
    }
 
  reverse(l.begin(),l.end());
 
  for(int i = (int)u.size()-2; i >= 1; i--) l.push_back(u[i]);
 
  return l;
}


// -- Library 2d -- END   ----------------------------------------

struct Data {
  int cur;
  double v;
  bool operator < (const Data &data) const {
    //return v > data.v;
    return LT(data.v,v); // data.v < v
  }
};

double dijkstra(vector<vector<Edge>> &G,int s,int t) {

  int V = G.size();
  vector<double> mini(V,1e8);

  priority_queue<Data> Q;
  Q.push((Data){s,0});
  mini[s] = 0;

  while( !Q.empty() ) {
    Data data = Q.top(); Q.pop();
    rep(i,(int)G[data.cur].size()) {
      Edge &e = G[data.cur][i];
      if( LT(data.v+e.cost,mini[e.to]) ) {
	mini[e.to] = data.v + e.cost;
	Q.push((Data){e.to,mini[e.to]});
      }
    }
  }
  assert( !equals(mini[t],1e8) );
  return mini[t];
}

inline bool is_parallel(Segment s,Segment t) { return equals(cross(s.p1-s.p2,t.p1-t.p2),0.0); }

const double rho = 1e5;
Polygon calc_shadow(Polygon &poly,double H,double theta,double phi) {
  assert( poly.size() >= 3 );
  assert( !( is_parallel(Segment(poly[0],poly[1]),Segment(poly[1],poly[2])) ) );
  
  theta = toRad(theta);
  phi   = toRad(phi);
  if( debug ) {
    cout << "theta = " << theta << " and phi = " << phi << endl;
  }
  Polygon shadow_poly;
  rep(i,(int)poly.size()) {
    double r = H / tan(phi);
    Vector v = Vector(r*cos(theta),r*sin(theta));
    shadow_poly.push_back(poly[i]-v);
  }
  rep(i,(int)poly.size()) shadow_poly.push_back(poly[i]);
  return andrewScan(shadow_poly);
}

void compute(vector<double> &H,vector<Polygon> &polies,vector<Segment> &segs,double theta,double phi,Point &s,Point &t) {  
  // 影のグラフを構築
  vector<Polygon> shadow_polies;
  rep(i,(int)polies.size()) {
    shadow_polies.push_back(calc_shadow(polies[i],H[i],theta,phi));
    if( debug ) {
      puts("");
      cout << i << "-th polygon (MINE)" << endl;
      rep(j,(int)shadow_polies[i].size()) {
	cout << shadow_polies[i][j] << endl;
      }
      puts("");
    }
  }
  
  if( debug ) {
    puts("");
    cout << "--- graph construction was finished! ---" << endl;
    puts("");
  }

  // 交点列挙
  vector<Point> ps;
  rep(i,(int)segs.size()) {
    ps.push_back(segs[i].p1);
    ps.push_back(segs[i].p2);
  }
  ps.push_back(s);
  ps.push_back(t);
  rep(i,(int)shadow_polies.size()) {
    int V = shadow_polies[i].size();
    rep(j,V) {
      Segment seg = Segment(shadow_polies[i][j],shadow_polies[i][(j+1)%V]);
      rep(k,(int)segs.size()) {
	if( intersectSS(seg,segs[k]) ) {
	  ps.push_back(crosspoint(seg,segs[k]));
	}
      }
    }
  }

  if( debug ) {
    puts("");
    cout << "--- enumeration was finished! ---" << endl;
    puts("");
  }

  // グラフ構築
  vector<vector<Edge>> G = segmentArrangement(segs,ps);
  assert( G.size() == ps.size() );

  if( debug ) {
    puts("");
    cout << "--- arrangement was finished! ---" << endl;
    puts("");
  }

  // 重み付け (中間地点が影に含まれてたら重み0)
  rep(i,(int)G.size()) {
    rep(j,(int)G[i].size()) {
      Edge &e = G[i][j];
      Point mp = (ps[e.from]+ps[e.to])/2.0;
      bool in_shadow = false;
      rep(k,(int)shadow_polies.size()) {
	if( inPolygon(shadow_polies[k],mp) ) {
	  in_shadow = true;
	  break;
	}
      }
      if( in_shadow ) {
	e.cost = 0;
      }
    }
  }

  if( debug ) {
    puts("");
    cout << "--- weighting was finished! ---" << endl;
    puts("");
  }
  
  // dijkstra
  int sp=-1,ep=-1;
  rep(i,(int)ps.size()) {
    if( ps[i] == s ) sp = i;
    if( ps[i] == t ) ep = i;
  }
  assert( sp != -1 );
  assert( ep != -1 );
  printf("%.10f\n",dijkstra(G,sp,ep));
}


int main() {
  int N,M;
  while( cin >> N >> M, N|M ) {
    vector<double> H(N);
    vector<Polygon> polies(N);
    rep(i,N) {
      int V;
      cin >> V >> H[i];
      polies[i].resize(V);
      rep(j,V) cin >> polies[i][j].x >> polies[i][j].y;
    }
    vector<Segment> segs(M);
    rep(i,M) cin >> segs[i].p1.x >> segs[i].p1.y >> segs[i].p2.x >> segs[i].p2.y;
    double theta,phi;
    cin >> theta >> phi;
    Point s,t;
    cin >> s.x >> s.y;
    cin >> t.x >> t.y;

    /*
    rep(i,N) {
      int V = polies[i].size();
      rep(j,V) {
	Segment seg = Segment(polies[i][j],polies[i][(j+1)%V]);
	rep(k,M) {
	  assert( !intersectSS(seg,segs[k]) );
	}
      }
    }
    */

    compute(H,polies,segs,theta,phi,s,t);
  }
  return 0;
}

TCO17 Algorithm Round 2A DistanceZeroAndOne

以下備忘録
解法:
d0,d1からグラフを構築し、それらのグラフを1つにまとめる
ある辺を答えとなるグラフに追加するかどうかは次のように決定する
1. d0,d1から構築したグラフの両方がその辺を持つなら答えのグラフにもその辺を追加
2. 片方しか持たないなら、その辺を追加した場合もう片方の最短距離に影響しないかどうかをチェックし、影響しないなら追加
最後にこうしてできたグラフの0,1からの最短距離がd0,d1と一致するか調べる

コード:

#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)

using namespace std;

typedef long long ll;
const int IINF = INT_MAX;

class DistanceZeroAndOne {
public:
  
  vector<int> bfs(vector<string> &G,int sp) {
    int n = G.size();
    vector<int> mini(n,IINF);
    mini[sp] = 0;
    deque<int> deq;
    deq.push_back(sp);
    while( !deq.empty() ) {
      int c = deq.front(); deq.pop_front();
      rep(i,n) {
        if( G[c][i] == 'Y' )  {
          if( mini[c] + 1 < mini[i] ) {
            mini[i] = mini[c] + 1;
            deq.push_back(i);
          }
        }
      }
    }
    return mini;
  }

  bool check(vector<string> &G,vector<int> &d0,vector<int> &d1) {
    vector<int> m0 = bfs(G,0);
    vector<int> m1 = bfs(G,1);
    if( m0 == d0 && m1 == d1 ) return true;
    return false;
  }

  void makeG(vector<string> &G,vector<int> &d,int sp) {
    int n = G.size();
    vector<int> pre;
    pre.push_back(sp);
    REP(i,1,n) {
      vector<int> nex;
      rep(j,n) if( d[j] == i ) {
        nex.push_back(j);
        rep(k,(int)pre.size()) {
          G[j][pre[k]] = 'Y';
          G[pre[k]][j] = 'Y';
        }
      }
      pre = nex;
    }
  }

  vector <string> construct( vector <int> d0, vector <int> d1 ) {
    int n = d0.size();
    vector<string> ans(n);
    vector<string> G0(n),G1(n);
    rep(i,n) ans[i] = G0[i] = G1[i] = string(n,'N');
    makeG(G0,d0,0);
    makeG(G1,d1,1);
    rep(i,n) {
      rep(j,n) {
        if( G0[i][j] != G1[i][j] ) {
          char tmp = 'N';
          if( G0[i][j] == 'Y' ) {
            int vi = d1[i], vj = d1[j];
            if( vi + 1 < vj || vj + 1 < vi )  tmp = 'N';
            else tmp = 'Y';
          } else {
            int vi = d0[i], vj = d0[j];
            if( vi + 1 < vj || vj + 1 < vi )  tmp = 'N';
            else tmp = 'Y';
          }
          ans[i][j] = tmp;
        } else {
          ans[i][j] = G0[i][j];
        }
      }
    }
    if( check(ans,d0,d1) ) return ans;
    return vector<string>();
  }
};

TCO17 Algorithm Round 2A Easy FoxAndCake2

以下備忘録

解法 :
与えられる整数c,sの偶奇で場合分けする。
1. cとsがともに偶数ならば Possible ( 分割しなくても gcd(c,s) != 1 なので
2. cとsがともに奇数ならば
- c と s がともに 3 以上なら Possible ( (3,3) (c-3,s-3)とすればよい
- そうでないなら、 c か s の一方は1なので Impossible
3. cとsのうち片方だけが偶数ならば
- これを (偶、偶)(偶、奇)に分割する
  奇はできるだけ小さくしたほうがよいので 3
偶も(3とのgcdが1でないようなもののなかで)できるだけ小さくしたほうがよいので 6
  cとsから(3,6)のペアを作ったあとにc,sがそれぞれ2以上なら Possible
それ以外は Impossible

コード :

const string YES = "Possible";
const string NO  = "Impossible";

class FoxAndCake2 {
public:
  string isPossible( int _c, int _s ) {
    ll c = _c, s = _s;
    if( __gcd(c,s) != 1LL ) return YES;
    if( c % 2 == 1 && s % 2 == 1 && min(c,s) >= 3 )  return YES;
    ll even = c ,odd = s;
    if( !( even % 2LL == 0 ) ) swap(even,odd);
    odd -= 3LL;
    even -= 6LL;
    if( odd >= 2LL && even >= 2LL ) return YES;
    return NO;
  }
};

SRM 714 div 1 easy : ParenthesisRemoval

問題概要 :
正しい対応関係を持つ括弧の文字列が与えられる
その文字列の先頭の開き括弧 ( とそれより後ろにあるいずれかの閉じ括弧 ) を削除する
削除後の文字列の括弧の対応関係も正しくなければならない
このとき、括弧の消し方は何通りあるか?

解法:
後ろから数え上げる
文字列を後ろから見ていき、 開き括弧がきたら
それより後ろの閉じ括弧の数 - それより後ろの開き括弧の数 # コメント : これはいま見ている開き括弧と一緒に消すことができる括弧の数
を答えにどんどんかけていく

コード :

#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS 1e-10
#define equals(a,b) (fabs((a)-(b)) < EPS)

using namespace std;

typedef long long ll;

const int IINF = INT_MAX;

const ll mod = 1000000007;

class ParenthesisRemoval {
public:
  int countWays(string s) {
    ll ans = 1;
    int L = 0, R = 0;
    for(int i=(int)s.size()-1;i>=0;--i) {
      if( s[i] == ')' ) ++R;
      else {
	ans = ( ans * ( R - L ) ) % mod;
	++L;
      }
    }
    return ans;
  }
};