AOJ 0284 : Happy End Problem

問題リンク : Happy End Problem | Aizu Online Judge

問題概要 : 略

解法 :
DP
答えるものが多角形の頂点番号なので経路復元もしなければならない
メモ化再帰でも間に合うと思ったが幻想だった
各点を出力に合わせてソートする
下にあるものほど前にくる
同じならば左にあるものが前にくる
各点から偏角ソートしておく
dp[現在の多角形の頂点数][初期位置][最後に追加した点の位置][最後の１つ前に追加した位置] := 面積の最小値

コード :

#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)
#define COUNTER_CLOCKWISE 1
#define CLOCKWISE -1 
#define ONLINE_BACK 2
#define ONLINE_FRONT -2
#define ON_SEGMENT 0

using namespace std;

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

class Point{
public:
  double x,y;
  int index;
  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(y,p.y) ? LT(y,p.y) : LT(x,p.x); }
  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }
};

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)); }

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;
}

inline double getArea(Point A,Point B,Point C){
  double a = abs( A - B );
  double b = abs( B - C );
  double c = abs( C - A );
  double s = ( a + b + c ) / 2.0;
  return sqrt( s * ( s - a ) * ( s - b ) * ( s - c ) );
}

const int MAX_N = 41;
const double DINF = 1e30;
int N;
double dp[MAX_N][MAX_N][MAX_N][MAX_N]; // dp[cnt][sp][ep_0][ep_1] := minimum
int  prev[MAX_N][MAX_N][MAX_N][MAX_N]; // 
double area[MAX_N][MAX_N][MAX_N];
Point ps[MAX_N];
vector<Point> vec[MAX_N];

Point vir;
bool comp_angle(Point a,Point b){
  double rad_a = atan2(a.y-vir.y,a.x-vir.x);
  double rad_b = atan2(b.y-vir.y,b.x-vir.x);
  return LT(rad_a,rad_b);
}

map<Point,int> toIndex; 
double mini[MAX_N];
Point a[MAX_N],b[MAX_N],c[MAX_N];

int main(){
  scanf("%d",&N);
  rep(i,N) scanf("%lf %lf",&ps[i].x,&ps[i].y);

  vector<int> ori(N);
  vector<Point> orip;
  rep(i,N) orip.push_back(ps[i]);
  rep(i,N+1) mini[i] = DINF;
  rep(i,N+1) rep(j,N) rep(k,N) rep(l,N) dp[i][j][k][l] = DINF, prev[i][j][k][l] = -1;
  sort(ps,ps+N);
  toIndex.clear();
  rep(i,N) toIndex[ps[i]] = i;
  rep(i,N) ori[toIndex[orip[i]]] = i;

  rep(sp,N){
    vec[sp].clear();
    REP(i,sp+1,N) vec[sp].push_back(ps[i]);
    if( vec[sp].size() < 2 ) break;
    vir = ps[sp];
    sort(vec[sp].begin(),vec[sp].end(),comp_angle);     
  }

  rep(sp,N) rep(ep_0,vec[sp].size())  REP(ep_1,ep_0+1,vec[sp].size()) {
    dp[3][sp][toIndex[vec[sp][ep_0]]][toIndex[vec[sp][ep_1]]] = getArea(ps[sp],vec[sp][ep_0],vec[sp][ep_1]);
    if( LT(dp[3][sp][toIndex[vec[sp][ep_0]]][toIndex[vec[sp][ep_1]]],mini[3])) {
      mini[3] = dp[3][sp][toIndex[vec[sp][ep_0]]][toIndex[vec[sp][ep_1]]];
      a[3] = ps[sp];
      b[3] = vec[sp][ep_1];
      c[3] = vec[sp][ep_0];
    }
  }

  REP(cnt,3,N) {
    rep(sp,N) {
      rep(ep_0,vec[sp].size()){
        REP(ep_1,ep_0+1,vec[sp].size()) if( !equals(dp[cnt][sp][toIndex[vec[sp][ep_0]]][toIndex[vec[sp][ep_1]]],DINF) ) {
          REP(next,ep_1+1,vec[sp].size()){
            int res = ccw(vec[sp][ep_0],vec[sp][ep_1],vec[sp][next]);
            if( !( res == COUNTER_CLOCKWISE || res == ONLINE_FRONT) ) continue;
            double temp = dp[cnt][sp][toIndex[vec[sp][ep_0]]][toIndex[vec[sp][ep_1]]] + getArea(ps[sp],vec[sp][ep_1],vec[sp][next]);
            if( LT(temp,dp[cnt+1][sp][toIndex[vec[sp][ep_1]]][toIndex[vec[sp][next]]]) ) {
              dp[cnt+1][sp][toIndex[vec[sp][ep_1]]][toIndex[vec[sp][next]]] = temp;
              prev[cnt+1][sp][toIndex[vec[sp][ep_1]]][toIndex[vec[sp][next]]] = toIndex[vec[sp][ep_0]];
              if( LT(temp,mini[cnt+1]) ) {
                mini[cnt+1] = temp;
                a[cnt+1] = ps[sp];
                b[cnt+1] = vec[sp][ep_1];
                c[cnt+1] = vec[sp][next];
              }
            }
          }
        }
      }
    }
  }

  int Q;
  cin >> Q;
  rep(_,Q){
    int use;
    scanf("%d",&use);
    if( equals(mini[use],DINF) ) puts("NA");
    else {
      vector<int> path;
      int ia = toIndex[a[use]];
      int ib = toIndex[b[use]];
      int ic = toIndex[c[use]];
      path.push_back(ori[ic]);
      path.push_back(ori[ib]);
      int remain = use;
      while( remain >= 3 && prev[remain][ia][ib][ic] != -1 ){
        int next = prev[remain][ia][ib][ic];
        path.push_back(ori[next]);
        ic = ib;
        ib = next;
        --remain;
      }

      if( use > 3 ) {
        path.push_back(ori[toIndex[a[use]]]);
        reverse(path.begin(),path.end());
      } else {
        vector<int> t;
        t.push_back(ori[toIndex[a[use]]]);
        t.push_back(ori[toIndex[c[use]]]);
        t.push_back(ori[toIndex[b[use]]]);
        path = t;
      }
      rep(i,path.size()) {
        if( i ) printf(" ");
        printf("%d",path[i]+1);
      } puts("");
    }
  }
  return 0;
}