Entries

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

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

コメント

コメントの投稿

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

トラックバック

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

April 2011 Challenge 1問目 - Central Point (この記事を編集する[管理者用])

Source

codechef April 2011 Challenge 1問目
Problem Statement
April 2011 Challengeの参加記録

問題概要

2次元平面上のN点 (2000以下) が与えられる.
そのN点までのユークリッド距離の和が最小となる格子点を選んだ時,その距離の和を求める問題.
与えられる点の座標の絶対値は10^9以下.

解法

最初,格子点でなく,整数座標以外の点も許して,山登り法で適当に和が小さくなる方に動かしていく.
後は,その周辺の格子点を全部調べる.
想定解法は,徐々に近づくていくのではなく,3分探索を使うみたい.

C言語のスパゲッティなコード
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#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;}

#define RED 0.85

int main(){
  int i,j,k,l,m,n;
  double sum, tmp, res;
  double go;
  double dx[8]={1,-1,0,0,0.2,0.7,-0.4,-0.9};
  double dy[8]={0,0,1,-1,-0.7,0.2,0.9,-0.4};
  pnt p[2100], now;

  for(;;){
    scanf("%d",&n); if(!n)break;
    rep(i,n) scanf("%lf%lf",&p[i].x,&p[i].y);
    res = 1e100;

    now = pntGenerator(0,0);
    rep(i,n) now = pntPlus(now,p[i]);
    rep(i,n) now = pntMultipleDouble(now,1.0/n);

    go = 1e9;
    sum = 0; rep(k,n) sum += pntDistance(now,p[k]);
    while(go > 1e-5){
      rep(m,8){
        now.x += dx[m] * go;
        now.y += dy[m] * go;

        tmp = 0;
        rep(k,n) tmp += pntDistance(now,p[k]);

        if(tmp > sum){
          now.x -= dx[m] * go;
          now.y -= dy[m] * go;
        } else {
          sum = tmp;
        }
      }
      go *= RED;
    }

    now.x = (int) now.x;
    now.y = (int) now.y;
    REP(i,-5,6) REP(j,-5,6){
      tmp = 0;
      now.x += i; now.y += j;
      rep(k,n) tmp += pntDistance(now,p[k]);
      now.x -= i; now.y -= j;

      if(res > tmp) res = tmp;
    }

    printf("%.6lf\n",res);
  }


  return 0;
}

コメント

コメントの投稿

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

トラックバック

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

Appendix

Recent Articles

ブログ内検索

Ads


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