2017-12-26

POJ 2420 A Star not a Tree?

POJ 数理最適化

問題

二次元座標上に点がN個与えられる。この座標上から一点を指定し、N個の点との距離の総和を最小化せよ。

やりかた

解析的に答えが求まりそうにないので、数理最適化。自分は最急降下法でやった。
選んだ座標を $(X, Y)$ とし、与えられたN個の点の座標を $(a_i, b_i)$ すると距離の総和は $L = \sum \sqrt{(X - a_i)^2 + (Y - b_i)^2}$ となる。

なのでこの関数の勾配方向は
$\left( \frac{\partial L}{\partial X}, \frac{\partial L}{\partial Y} \right) = \left( \sum \frac{X-a_i}{\sqrt{(X - a_i)^2 + (Y - b_i)^2}}, \sum \frac{Y-b_i}{\sqrt{(X - a_i)^2 + (Y - b_i)^2}} \right)$ となるので、ここから更新式 $X \leftarrow X - \alpha * {\partial L}/{\partial X}, \quad Y \leftarrow Y - \alpha * {\partial L}/{\partial Y}$ に従って予め決めた回数だけX, Yを更新する。

関数が凸関数でないような気がしたため、局所最小値に入ってしまうかもしれないと思ったが、初期値を重心にすれば大丈夫だろとたかをくくって出した。

以下ソース。

int x[100], y[100];

int main(){
  int N; cin >> N;

  double sumX = 0, sumY = 0;
  for(int i = 0; i < N; i++){
    cin >> x[i] >> y[i];
    sumX += x[i], sumY += y[i];
  }
  
  double X = sumX / N;
  double Y = sumY / N;

  for(int i = 0; i < 100000; i++){
    double alpha = 1e-2;

    double dx = 0, dy = 0;
    for(int j = 0; j < N; j++){
      dx += (X - x[j]) / sqrt(pow(X - x[j], 2.0) + pow(Y - y[j], 2.0));
      dy += (Y - y[j]) / sqrt(pow(X - x[j], 2.0) + pow(Y - y[j], 2.0));
    }
    
    X = X - alpha * dx;
    Y = Y - alpha * dy;
  }

  double L = 0;
  for(int j = 0; j < N; j++)
    L += sqrt(pow(X - x[j], 2.0) + pow(Y - y[j], 2.0));
  printf("%d\n", (int)(L + 0.5));
  return 0;
}

最初、学習率が小さすぎて最小値までたどりつかずWA。提出する前は収束するかきちんと確かめた上で、制限時間内で収束するような学習率にするべきだとわかった。（あたりまえだが）

2017-12-18

POJ 3422 Kaka's Matrix Travels

POJ 最小費用流

問題

NxNのグリッドがあり、各マスに0以上の整数が書かれている。グリッドの左上のマスからスタートして右か下に一つずつ移動して右下のマスまでたどり着きたい。通ったマスの数字は得点となり、一度通ったらそのマスの数字は0になる。K回スタートからゴールに移動したときに得られる総得点を最大化せよ。

やりかた

最小（最大）費用流。単純に最大経路をK回調べるのでは駄目。
費用流の問題として考えるとソースから左上のマスにつながり、右下のマスからシンクにつながるネットワークにKだけ流すときの最大費用流として考えることができる。ここで問題になるのは各マスが一回通ったあとは得点が0になることであり、ネットワークフローで考えればこれは（得点があるときの）そのセルに流量1の制限があることと同じになる。そこで下の図のように各マスのコピーを作り、コピーに対して1の流量のある得点がつく辺と、流量無限大のコスト0の辺を加えてやることで、得点のある辺は一回しか使われないようにすることができる。
f:id:Area1:20171218005854p:plain

以下ソース。

int N, K;
struct min_cost_flow{
  struct edge{
    int to, cap, cost, rev;
    edge(int to, int cap, int cost, int rev) : 
      to(to), cap(cap), cost(cost), rev(rev) {};
  };
  vector<vector<edge> > G;

  min_cost_flow(int V) : G(V) {}

  void add_edge(int from, int to, int cap, int cost){
    G[from].push_back(edge(to, cap, cost, G[to].size()));
    G[to].push_back(edge(from, 0, -cost, G[from].size() - 1));
  }

  int flow(int s, int t, int f){
    typedef pair<int, int> P;
    int V = G.size();
    vector<int> po(V);
    int res = 0;
    while(f > 0){
      vector<int> dist(V, INF), prevv(V, -1), preve(V, -1);
      priority_queue<P, vector<P>, greater<P> > Q;
      Q.push(make_pair(0, s));
      dist[s] = 0;
      while(!Q.empty()){
	P q = Q.top(); Q.pop();
	int u = q.second;
	if(dist[u] < q.first) continue;
	for(int i = 0; i < G[u].size(); i++){
	  edge &e = G[u][i];
	  if(e.cap > 0 && dist[e.to] > dist[u] + e.cost + po[u] - po[e.to]){
	    dist[e.to] = dist[u] + e.cost + po[u] - po[e.to];
	    prevv[e.to] = u;   preve[e.to] = i;
	    Q.push(P(dist[e.to], e.to));
	  }
	}
      }
      if(dist[t] == INF) return -1;
      for(int v = 0; v < V; v++) po[v] += dist[v];
      
      int scap = f;
      for(int v = t; v != s; v = prevv[v])
	scap = min(scap, G[prevv[v]][preve[v]].cap);
      
      f -= scap;
      res += scap * po[t];
      for(int v = t; v != s; v = prevv[v]){
	edge &e = G[prevv[v]][preve[v]];
	e.cap -= scap;
	G[v][e.rev].cap += scap;
      }
    }
    return res;
  }
};

int main(){
  scanf("%d %d", &N, &K);

  int V = N * N;
  int s = 2 * V, t = s + 1;
  
  min_cost_flow mcf(2 * V + 2);
  for(int i = 0; i < N; i++){
    for(int j = 0; j < N; j++){
      int x;
      scanf("%d", &x);

      mcf.add_edge(i * N + j, i * N + j + V, 1, -x);//最小費用流でとくため符号を反転
      mcf.add_edge(i * N + j, i * N + j + V, INF, 0);
      
      if(j != N - 1)
	mcf.add_edge(i * N + j + V, i * N + j + 1, INF, 0);
      if(i != N - 1)
      	mcf.add_edge(i * N + j + V, (i + 1) * N + j, INF, 0);
    }
  }
  mcf.add_edge(s, 0, INF, 0);
  mcf.add_edge(2 * V - 1, t, INF, 0);
  
  printf("%d\n", -mcf.flow(s, t, K));////最小費用流でとくため符号を反転
  return 0;
}

図が薄くてすみません。

2017-12-13

POJ 3790 Recursively Palindromic Partitions

POJ 分割統治深さ優先探索

問題

与えられた自然数を、それ以下の自然数の数列の和で表現する。これを自然数の分割と呼ぶ。回文となるような分割をPalindromic Partitionという。そのうち以下の条件を満たす分割をRecursively Palindromic Partitionsという。
条件：分割の左右それぞれがRecursively Palindromic Partitionsである。ただし分割数1の分割は常にRecursively Palindromic Partitionである。

与えられる自然数のRecursively Palindromic Partitionsの数をもとめよ。

やりかた

分割統治で真ん中から左右に分けていき、分けきれなくなるまでわけられたらひとつと数える。メモ化して計算量を落とす。

以下ソース。

int N;
ll dp[1001][1001];

ll rec(int idx, int rem){
  if(rem == 0) return 1;
  if(dp[idx][rem] != -1) return dp[idx][rem];
  
  ll ret = 0;
  for(int i = 0; i <= rem; i++){
    if((rem - i) % 2) continue;
    ret += rec(idx + 1, (rem - i) / 2);
  }
  return dp[idx][rem] = ret;
}

int main(){
  int T;
  cin >> T;
  for(int i = 1; i <= T; i++){
    cin >> N;
    FILL(dp, -1);
    cout << i << " " << rec(0, N) << endl;
  }
  return 0;
}

2017-11-17

PKU 3182 The Grove

POJ 幾何

問題

グリッド上に連結した池がある。これをぐるっと回って戻ってくるための最短経路を求めよ。ただし8連結方向に移動できる。

やりかた

通ったけど多分嘘解法。

凸包。まず池の周りを4連結で囲むように印をつける。その印にたいして凸包を計算して池を一周するための最短経路数を求める。
続いて凸包で囲まれる円周の一辺を取り除いて始点を迂回するルートにした場合の経路長を計算する（下図の点線ルート、もしくは波線ルート）。この取り除く辺を円周の全ての辺に対して試してみて、最小値を計算する。

正しくないような気がしたが、これで通った。

以下ソース。

const int NUM = 100;
int dy[] = {1, 0, -1, 0, 1, -1, -1, 1};
int dx[] = {0, 1, 0, -1, 1, 1, -1, -1};
int H, W;
int d[NUM][NUM];
char f[NUM][NUM];

int bfs(int sy, int sx, int gy, int gx){
  FILL(d, -1);
  d[sy][sx] = 0;
  queue<int> Q;
  Q.push(sy * W + sx);
  while(!Q.empty()){
    int y = Q.front() / W, x = Q.front() % W;
    Q.pop();
    if(gy == y && gx == x) return d[y][x];
    for(int k = 0; k < 8; k++){
      int ny = y + dy[k], nx = x + dx[k];
      if(0 <= ny && ny < H && 0 <= nx && nx < W &&
	 d[ny][nx] == -1 && f[ny][nx] == '.'){
	d[ny][nx] = d[y][x] + 1;
	Q.push(ny * W + nx);
      }
    }
  }
}


int add(int a, int b){
  if(abs(a + b) < EPS * (abs(a) + abs(b))) return 0;
  return a + b;
}

struct P{
  int x, y;
  P() {}
  P(int x, int y) : x(x), y(y) {}
  P operator + (P p){
    return P(add(x, p.x), add(y, p.y));
  }
  P operator - (P p){
    return P(add(x, -p.x), add(y, -p.y));
  }
  P operator * (int d){
    return P(x * d, y * d);
  }
  int dot(P p){
    return add(x * p.x, y * p.y);
  }
  int det(P p){
    return add(x * p.y, -y * p.x);
  }
};

bool cmp_x(P a, P b){
  if(a.x != b.x) return a.x < b.x;
  return a.y < b.y;
}

vector<P> convex_hull(vector<P> &p){
  int N = p.size();
  sort(p.begin(), p.end(), cmp_x);
  int k = 0;
  vector<P> ret(N * 2);
  for(int i = 0; i < N; i++){
    while(k > 1 && (ret[k - 1] - ret[k - 2]).det(p[i] - ret[k - 1]) <= 0) k--;
    ret[k++] = p[i];
  }
  for(int i = N - 2, t = k; i >= 0; i--){
    while(k > t && (ret[k - 1] - ret[k - 2]).det(p[i] - ret[k - 1]) <= 0) k--;
    ret[k++] = p[i];
  }
  ret.resize(k - 1);
  return ret;
}

int main(int argc, char **argv){
  vector<P> p;
  cin >> H >> W;
  int sy, sx;
  for(int i = 0; i < H; i++){
    for(int j = 0; j < W; j++){
      cin >> f[i][j];
      if(f[i][j] == '*')
	sy = i, sx = j;
    }
  }
  for(int i = 0; i < H; i++){
    for(int j = 0; j < W; j++){
      for(int k = 0; k < 4; k++){
	int ny = i + dy[k], nx = j + dx[k];
	if(0 <= ny && ny < H && 0 <= nx && nx < W &&
	   f[ny][nx] == 'X'){
	  P tmp;
	  tmp.y = i, tmp.x = j;
	  p.push_back(tmp);
	  break;
	}
      }
    }
  }

  bfs(sy, sx, -1, -1);
  
  vector<P> ch = convex_hull(p);
  int M = ch.size();

  int round = 0;
  for(int i = 0; i < M; i++){
    f[ch[i].y][ch[i].x] = '@';
    round += max(abs(ch[i].y - ch[(i + 1) % M].y),
		 abs(ch[i].x - ch[(i + 1) % M].x));
  }

  int ans = INF;
  for(int i = 0; i < M; i++){
    int l = max(abs(ch[i].y - ch[(i + 1) % M].y),
		abs(ch[i].x - ch[(i + 1) % M].x));
    int hige = d[ch[i].y][ch[i].x]
      + d[ch[(i + 1) % M].y][ch[(i + 1) % M].x];
    ans = min(ans, round - l + hige);
  }
  cout << ans << endl;
  return 0;
}

アップした図はiPad pro+Applepenで書いてみたもの。これは便利！

2017-11-17

POJ 1387 Labyrinth

POJ グラフ幅優先探索

問題

やりかた

木の直径を求める問題。グラフ上の適当な点からはじめて最も遠い点を見つける。続いてその点から最も遠い点を探すとその間が直径になる。
メモリ量、制限時間ともにシビア。グラフは隣接リストで持たずに、ある頂点に対する４方向の連結情報をビット状態で持つようにして、一つの頂点の情報をcharで表現してメモリを減らした。

以下ソース。

int R, C;
char f[1001][1001];
int dy[] = {1, 0, -1, 0};
int dx[] = {0, 1, 0, -1};

char G[1000001];
int d[1000001];

void bfs(int idx){
  memset(d, -1, sizeof(d));

  queue<int> Q;
  Q.push(idx);
  d[idx] = 0;
  while(!Q.empty()){
    int x = Q.front(); Q.pop();
    for(int k = 0; k < 4; k++){
      if(G[x] >> k & 1){
	int v = x + dy[k] * C + dx[k];
	if(d[v] == -1){
	  Q.push(v);
	  d[v] = d[x] + 1;
	}
      }
    }
  }

  return;
}

int main(){
  int T;
  scanf("%d", &T);
  while(T--){
    scanf("%d %d", &C, &R);
    for(int i = 0; i < R; i++){
      for(int j = 0; j < C; j++){
	f[i][j] = getchar();
      }
      getchar();
    }
    memset(G, 0, sizeof(G));
    
    int start;
    for(int i = 0; i < R; i++){
      for(int j = 0; j < C; j++){
	if(f[i][j] == '.'){
	  start = i * C + j;
	  for(int k = 0; k < 4; k++){
	    int ny = i + dy[k], nx = j + dx[k];
	    if(0 <= ny && ny < R
	       && 0 <= nx && nx < C && f[ny][nx] == '.'){
	      G[i * C + j] |= 1 << k;
	    }
	  }
	}
      }
    }

    //木の直径
    bfs(start);
    int fur = max_element(d, d + R * C) - d;
    bfs(fur);
    int ans = *max_element(d, d + R * C);
    printf("Maximum rope length is %d.\n", ans);
  }
  return 0;
}

2017-11-17

POJ 1376 Robot

POJ 幅優先探索

問題

グリッド上の格子点を直径1.6mの円形のロボットがスタート地点からゴール地点まで移動する。グリッドの一部に障害物があり、ロボットを障害物にぶつけてはいけない。ロボットは一回の移動で左右に向きを変えるか、1~3つ分の格子点を直進することができる。ゴール地点に行くための最短の移動回数を求めよ。

やりかた

グリッド上の位置、向きを状態にもつ幅優先。
以下ソース。

int H, W;
int sy, sx, gy, gx;
string dir;
bool f[60][60];
int d[60][60][4];
int dy[] = {-1, 0, 1, 0};
int dx[] = {0, 1, 0, -1};

bool in(int y, int x){
  if(0 < y && y < H && 0 < x && x < W){
    if(!f[y][x] && !f[y - 1][x - 1] &&
       !f[y - 1][x] && !f[y][x - 1]){
      return true;
    }
  }
  return false;
}

int main(){
  while(cin >> H >> W, H){
    for(int i = 0; i < H * W; i++){
      int v; cin >> v;
      f[i / W][i % W] = v == 1;
    }
    cin >> sy >> sx >> gy >> gx >> dir;

    FILL(d, INF);
    int c;
    if(dir == "north") c = 0;
    if(dir == "east") c = 1;
    if(dir == "south") c = 2;
    if(dir == "west") c = 3;
    queue<pair<int, int> > Q;
    Q.push(make_pair(sy * W + sx, c));
    d[sy][sx][c] = 0;
    while(!Q.empty()){
      int y = Q.front().first / W,
	x = Q.front().first % W,
	c = Q.front().second;
      Q.pop();

      if(y == gy && x == gx) break;
      
      for(int k = 1; k <= 3; k++){
	int ny = y + k * dy[c], nx = x + k * dx[c];
	if(in(ny, nx)){
	  if(d[ny][nx][c] == INF){
	    Q.push(make_pair(ny * W + nx, c));
	    d[ny][nx][c] = d[y][x][c] + 1;
	  }
	}else break;
      }
      if(d[y][x][(c + 1) % 4] == INF){
	Q.push(make_pair(y * W + x, (c + 1) % 4));
	d[y][x][(c + 1) % 4] = d[y][x][c] + 1;
      }
      if(d[y][x][(c + 3) % 4] == INF){
	Q.push(make_pair(y * W + x, (c + 3) % 4));
	d[y][x][(c + 3) % 4] = d[y][x][c] + 1;
      }      
    }

    int ans = INF;
    for(int s = 0; s < 4; s++)
      ans = min(ans, d[gy][gx][s]);
    cout << (ans >= INF ? -1 : ans) << endl;
  }
  return 0;
}

2017-10-23

POJ 2341 Spell Checker

POJ 文字列処理

問題

アルファベットと.,;:-!?とスペースからなる文章が与えられる．文の区切りは.か!か?である．文章中の不整合な箇所の個数を答えよ．
ただし不整合なケースとは以下の2種類を指す．

文頭が大文字でない
単語の途中に大文字がある

やりかた

文意を正確に汲み取る必要がある．
注意するのは以下のケース．

文頭が記号である場合はその次の文字を文頭として見る
単語中に記号がある場合は単語の区切りとみなす
改行は単語の区切りとみなす

以下ソース．

int main(){
  char c;
  bool s_start = true, w_start = true;
  int cnt = 0;

  while((c = getchar()) != -1){
    if(c == '\n'){
      w_start = true;
    }
    else if(c == '.' || c == '?' || c == '!'){
      s_start = w_start = true;
    }
    else if(c == ' ' || c == ':' || c == ',' || c == '-'){
      w_start = true;
    }
    else{
      if(s_start && w_start){

	if(!isalpha(c)) continue;
	if(islower(c)) cnt++;
	
	s_start = w_start = false;
	
      }else if(!s_start && w_start){

	s_start = w_start = false;
	
      }else if(!s_start && !w_start){

	if(isupper(c)) cnt++;
	
      }
    }
  }
  cout << cnt << endl;
  return 0;
}