Entries

スポンサーサイト (この記事を編集する[管理者用])

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

コメント

コメントの投稿

コメントの投稿
管理者にだけ表示を許可する

トラックバック

トラックバック URL
http://rsujskf.blog32.fc2.com/tb.php/2039-d37771e7
この記事にトラックバックする(FC2ブログユーザー)

SRM508 DIV1 HARD - BarbarianInvasion2 (この記事を編集する[管理者用])

Source

TopCoder SRM508 DIV1 HARD (1000pt)
Problem Statement
SRM508 DIV1 自分の参加記録

問題概要

厳密に凸なN角形 (50以下) と,その厳密に内部の点M個 (5以下) が与えられる.
N角形の周上に,等間隔に1000000!人 (凄い沢山) の人が並んでいる.
それらの人が,内部のM箇所の点に,1000000! / M人ずつ移動する.
各人の移動スピードは1.
必要時間の最小値を求める問題.

解法

答えに対して2分探索する.
M個の点を中心とする円と,N角形の各辺との交点で分割して,それぞれの線分がどの点に移動可能かを求める.
後は,最大流を流してやり,全て対応関係が取れるかどうかを調べる.
分割された線分の数は,最大で,N*(2M+1)ぐらいになるが,どの点に移動可能かで分別したら2^Mで抑えることができる.(が,やらなくても十分大丈夫)

C++によるスパゲッティなソースコード
// #includeとusing namespace std;は略

#define REP(i,a,b) for(i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)

#define EPS 1e-10
#define PI 3.141592653589793238462

#define MAX(a,b) (((a)>(b))?(a):(b))
#define MIN(a,b) (((a)<(b))?(a):(b))
typedef struct struct_point{double x,y;}pnt;
typedef struct struct_line{pnt a,b;}line;
typedef struct struct_circle{pnt c; double r;}circle;

int doubleSign(double a){if(a>EPS) return 1; if(a<-EPS) return -1; return 0;}
int doubleSignR(double a,double b){ if(doubleSign(b)!=0 && doubleSign(a/b-1)==0) return 0; return doubleSign(a-b);}
int pntSign(pnt a){int i=doubleSign(a.x); if(i) return i; return doubleSign(a.y);}

pnt pntPlus(pnt a,pnt b){a.x+=b.x; a.y+=b.y; return a;}
pnt pntMinus(pnt a,pnt b){a.x-=b.x; a.y-=b.y; return a;}
pnt pntMultiple(pnt a,pnt b){pnt c; c.x=a.x*b.x-a.y*b.y; c.y=a.x*b.y+a.y*b.x; return c;}
pnt pntMultipleDouble(pnt a,double k){a.x*=k; a.y*=k; return a;}

int pntIsEqual(pnt a,pnt b){return pntSign(pntMinus(a,b))==0;}

double pntLength(pnt a){return sqrt(a.x*a.x+a.y*a.y);}
double pntLength2(pnt a){return a.x*a.x+a.y*a.y;}
double pntDistance(pnt a,pnt b){return pntLength(pntMinus(a,b));}
double pntDistance2(pnt a,pnt b){return pntLength2(pntMinus(a,b));}
double pntArgument(pnt a){return atan2(a.y,a.x);}

pnt pntNormalize(pnt a){double n=pntLength(a); a.x/=n; a.y/=n; return a;}
double pntInnerProduct(pnt a,pnt b){return a.x*b.x+a.y*b.y;}
double pntOuterProduct(pnt a,pnt b){return a.x*b.y-a.y*b.x;}

pnt pntGenerator(double x,double y){pnt res; res.x=x; res.y=y; return res;}
line pntToLine(pnt a,pnt b){line res; res.a=a; res.b=b; return res;}
circle pntToCircle(pnt c,double r){circle res; res.c=c; res.r=r; return res;}

int isPntInCircle(pnt p,circle c){return doubleSign(c.r-pntDistance(p,c.c))+1;}
int isPntOnCircle(pnt p,circle c){if(doubleSign(pntDistance(p,c.c)-c.r)==0) return 1; return 0;}

/* a x^2 + b x + c = 0 を解く.戻り値は解の数.解の数が2ならres1 < res2. */
/* 解の数が無限なら return 3; */
int doubleSolveSecondDegreeEquation(double a,double b,double c,double *res1,double *res2){
  if(doubleSign(a)){
    double d=b*b-4*a*c; int m=doubleSign(d)+1;
    if(m==1) *res1=-b/(2*a);
    if(m==2){d=sqrt(d); *res1=(-b-d)/(2*a); *res2=(-b+d)/(2*a);}
    return m;
  }
  if(doubleSign(b)){*res1 = -b/c; return 1;}
  if(doubleSign(c)) return 0; return 3;
}

/* 線分 t (端は含む) と円 s の交点を求める.戻り値は交点の数.(3次元での球-直線の交点とほぼ同じコード) */
/* 2点で交わるなら,res1はt.aに近い側,res2はt.bに近い側 */
int pntIntersectionSegmentCircle(line t,circle s,double *res1,double *res2){
  int m; double a,b,c,t1,t2,k1,k2; pnt ab=pntMinus(t.b, t.a);
  a=b=0; c=-s.r*s.r;
  t1=t.a.x-s.c.x; t2=t.b.x-t.a.x; a+=t2*t2; b+=2*t1*t2; c+=t1*t1;
  t1=t.a.y-s.c.y; t2=t.b.y-t.a.y; a+=t2*t2; b+=2*t1*t2; c+=t1*t1;
  m=doubleSolveSecondDegreeEquation(a,b,c,&k1,&k2);
  if(m>=2) if(doubleSign(k2)==-1 || doubleSign(k2-1)==1) m--;
  if(m>=1) if(doubleSign(k1)==-1 || doubleSign(k1-1)==1) k1=k2, m--;
  if(m>=1) *res1 = k1;
  if(m>=2) *res2 = k2;
  return m;
}

void doubleSort(double d[],int s){int i,j;double k1,k2,t;if(s<=1)return;k1=(d[0]+d[s-1])/2.0;k2=k1+EPS;k1-=EPS;i=-1;j=s;for(;;){while(d[++i]<k1);while(d[--j]>k2);if(i>=j)break;t=d[i];d[i]=d[j];d[j]=t;}doubleSort(d,i);doubleSort(d+j+1,s-j-1);}



#define INT_INF 1000000000
#define DOUBLE_INF 1e300

#define DOUBLE_LIST_MAX_FLOW_NODE 2002
int *doubleLmfEdge[DOUBLE_LIST_MAX_FLOW_NODE];
double *doubleLmfFlow[DOUBLE_LIST_MAX_FLOW_NODE];
int *doubleLmfInv[DOUBLE_LIST_MAX_FLOW_NODE];
int doubleLmfEdgeSize[DOUBLE_LIST_MAX_FLOW_NODE];
int doubleLmfEdgeMemory[DOUBLE_LIST_MAX_FLOW_NODE];
int doubleLmfLevel[DOUBLE_LIST_MAX_FLOW_NODE];
int doubleLmfQueue[DOUBLE_LIST_MAX_FLOW_NODE];

void doubleListMaxFlowLevelize(int n,int st,int ed){
  int i,j,k,t;
  int q_st,q_end;
  rep(i,n) doubleLmfLevel[i]=-1; doubleLmfLevel[st]=0; q_st=0; q_end=1; doubleLmfQueue[0]=st;
  while(q_st!=q_end){
    i=doubleLmfQueue[q_st++]; t=doubleLmfLevel[i]+1;
    rep(j,doubleLmfEdgeSize[i]) if(doubleLmfFlow[i][j] > EPS){
      k=doubleLmfEdge[i][j]; if(doubleLmfLevel[k]!=-1) continue;
      doubleLmfLevel[k]=t; doubleLmfQueue[q_end++]=k; if(k==ed) return;
    }
  }
}

double doubleListMaxFlowFlow(int i,int ed,double lim){
  int j,k,ji; double s, t, ret = 0;
  if(i==ed) return lim;
  rep(j,doubleLmfEdgeSize[i]) if(doubleLmfFlow[i][j] > EPS){
    k=doubleLmfEdge[i][j]; if(doubleLmfLevel[k]<=doubleLmfLevel[i]) continue;
    s=lim; if(s>doubleLmfFlow[i][j]) s=doubleLmfFlow[i][j];
    t=doubleListMaxFlowFlow(k,ed,s); if(t < EPS) continue;
    ret+=t; lim-=t; ji=doubleLmfInv[i][j];
    doubleLmfFlow[i][j]-=t; doubleLmfFlow[k][ji]+=t;
    if(lim < EPS) break;
  }
  if(lim) doubleLmfLevel[i]=-1;
  return ret;
}

/* Dinic(Dinitz) */
/* 制約: n<=DOUBLE_LIST_MAX_FLOW_NODE, st!=ed */
double doubleListMaxFlow(int n,int st,int ed){
  double ret=0;
  for(;;){
    doubleListMaxFlowLevelize(n,st,ed); if(doubleLmfLevel[ed]==-1) break;
    ret += doubleListMaxFlowFlow(st,ed,DOUBLE_INF);
  }
  return ret;
}

void doubleListMaxFlowAllocInit(int s){
  int i;
  rep(i,DOUBLE_LIST_MAX_FLOW_NODE){
    doubleLmfEdgeMemory[i]=s;
    if(s){
      doubleLmfEdge[i] = (int*) malloc(s*sizeof(int));
      doubleLmfFlow[i] = (double*) malloc(s*sizeof(double));
      doubleLmfInv[i] = (int*) malloc(s*sizeof(int));
    }
  }
}

void doubleListMaxFlowAlloc(int i,int s){
  if(doubleLmfEdgeMemory[i]>=s) return;
  doubleLmfEdgeMemory[i]=s;
  free(doubleLmfEdge[i]); free(doubleLmfFlow[i]); free(doubleLmfInv[i]);
  doubleLmfEdge[i] = (int*) malloc(s*sizeof(int));
  doubleLmfFlow[i] = (double*) malloc(s*sizeof(double));
  doubleLmfInv[i] = (int*) malloc(s*sizeof(int));
}

void doubleListMaxFlowMemoryExpand(int i){
  int k, nw = doubleLmfEdgeMemory[i]*2;
  int *m1,*m3;
  double *m2;

  if(nw<3) nw=3;
  m1=(int*)malloc(nw*sizeof(int));
  m2=(double*)malloc(nw*sizeof(double));
  m3=(int*)malloc(nw*sizeof(int));

  if(doubleLmfEdgeMemory[i]){
    rep(k,doubleLmfEdgeSize[i]){
      m1[k]=doubleLmfEdge[i][k];
      m2[k]=doubleLmfFlow[i][k];
      m3[k]=doubleLmfInv[i][k];
    }
    free(doubleLmfEdge[i]); free(doubleLmfFlow[i]); free(doubleLmfInv[i]);
  }
  doubleLmfEdgeMemory[i]=nw;
  doubleLmfEdge[i]=m1;
  doubleLmfFlow[i]=m2;
  doubleLmfInv[i]=m3;
}

void doubleListMaxFlowEdgeInitAdv(int n){
  int i;
  rep(i,n) doubleLmfEdgeSize[i]=0;
}

void doubleListMaxFlowEdgeInit(){
  int i;
  rep(i,DOUBLE_LIST_MAX_FLOW_NODE) doubleLmfEdgeSize[i]=0;
}

/* 制約: n1!=n2 */
void doubleListMaxFlowEdgeAdd(int n1,int n2,double f1,double f2){
  int s1=doubleLmfEdgeSize[n1]++, s2=doubleLmfEdgeSize[n2]++;
  doubleLmfEdge[n1][s1]=n2; doubleLmfEdge[n2][s2]=n1;
  doubleLmfFlow[n1][s1]=f1; doubleLmfFlow[n2][s2]=f2;
  doubleLmfInv[n1][s1]=s2;  doubleLmfInv[n2][s2]=s1;
}

/* 制約: n1!=n2 */
void doubleListMaxFlowEdgeAddAdv(int n1,int n2,double f1,double f2){
  int s1=doubleLmfEdgeSize[n1]++, s2=doubleLmfEdgeSize[n2]++;
  if(s1 >= doubleLmfEdgeMemory[n1]) doubleListMaxFlowMemoryExpand(n1);
  if(s2 >= doubleLmfEdgeMemory[n2]) doubleListMaxFlowMemoryExpand(n2);
  doubleLmfEdge[n1][s1]=n2; doubleLmfEdge[n2][s2]=n1;
  doubleLmfFlow[n1][s1]=f1; doubleLmfFlow[n2][s2]=f2;
  doubleLmfInv[n1][s1]=s2;  doubleLmfInv[n2][s2]=s1;
}



class BarbarianInvasion2 {
public:
double minimumTime(vector <int> segX, vector <int> segY, vector <int> cityX, vector <int> cityY) {
  int i,j,k;
  double a,b,c;

  pnt city[10]; int city_num;
  line seg[60]; int seg_num;
  pnt p_tmp1, p_tmp2;
  double tmp1, tmp2;

  double split[66][666]; int split_size[66];
  int per_size[66], per_sum[66];

  int node, st, ed;
  double flow;
  double sum_len;

  doubleListMaxFlowAllocInit(700);

  city_num = cityX.size();
  rep(i,city_num) city[i] = pntGenerator( cityX[i], cityY[i] );

  seg_num = segX.size();
  rep(i,seg_num){
    pnt p1 = pntGenerator( segX[i],             segY[i] );
    pnt p2 = pntGenerator( segX[(i+1)%seg_num], segY[(i+1)%seg_num] );
    seg[i] = pntToLine( p1, p2 );
  }

  sum_len = 0;
  rep(i,seg_num) sum_len += pntDistance( seg[i].a, seg[i].b );

  a = 0; b = 1e6;

  while(b-a > EPS){
    c = (a+b)/2;

    rep(i,seg_num) split_size[i] = 0;
    rep(i,seg_num) rep(j,city_num){
      circle en;
      en.r = c; en.c = city[j];
      k = pntIntersectionSegmentCircle(seg[i],en,&tmp1,&tmp2);
      if(k>=1) split[i][split_size[i]++] = tmp1;
      if(k>=2) split[i][split_size[i]++] = tmp2;
    }

    rep(i,seg_num) split[i][split_size[i]++] = 0;
    rep(i,seg_num) split[i][split_size[i]++] = 1;
    
    rep(i,seg_num) doubleSort(split[i], split_size[i]);

    rep(i,seg_num){
      k = 1;
      REP(j,1,split_size[i]){
        if(split[i][k-1] < split[i][j]-EPS) split[i][k++] = split[i][j];
      }
      split_size[i] = k;
    }

    per_sum[0] = 0;
    rep(i,seg_num){
      per_size[i] = split_size[i] - 1;
      per_sum[i+1] = per_sum[i] + per_size[i];
    }

    node = city_num + per_sum[seg_num];
    st = node++; ed = node++;

    doubleListMaxFlowEdgeInitAdv(node);

    rep(i,city_num) doubleListMaxFlowEdgeAddAdv(st,i,sum_len/city_num,0);

    rep(i,seg_num) rep(j,per_size[i]){
      pnt oa = seg[i].a;
      pnt ob = seg[i].b;
      pnt pa = pntPlus( oa, pntMultipleDouble(pntMinus(ob,oa),split[i][j]  ) );
      pnt pb = pntPlus( oa, pntMultipleDouble(pntMinus(ob,oa),split[i][j+1]) );
      pnt pc = pntMultipleDouble( pntPlus(pa,pb), 0.5 );
      int ns = per_sum[i] + j;
      double len_seg = pntDistance(pa,pb);

      doubleListMaxFlowEdgeAddAdv(city_num+ns, ed, len_seg, 0);
      
      rep(k,city_num){
        double dist = pntDistance(city[k], pc);
        if(dist < c+EPS) doubleListMaxFlowEdgeAddAdv(k,city_num+ns,len_seg,0);
      }
    }

    flow = doubleListMaxFlow(node,st,ed);
    if(flow < sum_len-EPS) a=c; else b=c;
  }

  return (a+b)/2;
}

};

コメント

コメントの投稿

コメントの投稿
管理者にだけ表示を許可する

トラックバック

トラックバック URL
http://rsujskf.blog32.fc2.com/tb.php/2039-d37771e7
この記事にトラックバックする(FC2ブログユーザー)

Appendix

Recent Articles

ブログ内検索

Ads


(プライバシーポリシー)
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。