Cod sursa(job #196867)

Utilizator nusmaibunkeleviprofesor cicalescu nusmaibunkelevi Data 29 iunie 2008 19:16:33
Problema Adapost 2 Scor 80
Compilator cpp Status done
Runda Arhiva de probleme Marime 5.76 kb
#include<stdio.h>
#include<math.h>
#define LIM 0.001
#define NMAX 50000
#define APR 30
#define ex LIM
#define ey LIM

struct pct{double x,y;};
struct cmapct{pct p;int nr;double d,sd;};
pct v[NMAX];
cmapct w[NMAX/2];
int n,nrpa;

double dist(double x,double y){
int i;
double rez=0.0,patr,rad,term;
for(i=0;i<n;++i){
	patr=(x-v[i].x)*(x-v[i].x)+(y-v[i].y)*(y-v[i].y);
	if(patr<0) patr=-patr;
	rad=sqrt(patr);
	term=rad;
	rez+=term;
	}
return rez;
}

double dd(double x1,double y1,double x2,double y2){
double dx=x1-x2,dy=y1-y2;
return sqrt(dx*dx+dy*dy);
}

double fx(double x,double y){
int i;
double rez=0.0,patr,rad,term;
for(i=0;i<n;++i){
	patr=(x-v[i].x)*(x-v[i].x)+(y-v[i].y)*(y-v[i].y);
	rad=sqrt(patr);
	if(rad) term=(x-v[i].x)/rad;
	else term=0.0;
	rez+=term;
	}
return rez;
}

double fy(double x,double y){
int i;
double rez=0.0,patr,rad,term;
for(i=0;i<n;++i){
	patr=(x-v[i].x)*(x-v[i].x)+(y-v[i].y)*(y-v[i].y);
	rad=sqrt(patr);
	if(rad) term=(y-v[i].y)/rad;
	else term=0.0;
	rez+=term;
	}
return rez;
}
double dfx(double x,double y){
int i;
double rez=0.0,patr,rad,term;
for(i=0;i<n;++i){
	patr=(x-v[i].x)*(x-v[i].x)+(y-v[i].y)*(y-v[i].y);
	rad=sqrt(patr);
	if(rad) term=(y-v[i].y)*(y-v[i].y)/rad/patr;
	else term=0.0;
	rez+=term;
	}
return rez;
}

double dfy(double x,double y){
int i;
double rez=0.0,patr,rad,term;
for(i=0;i<n;++i){
	patr=(x-v[i].x)*(x-v[i].x)+(y-v[i].y)*(y-v[i].y);
	rad=sqrt(patr);
	if(rad) term=(x-v[i].x)*(x-v[i].x)/rad/patr;
	else term=0.0;
	rez+=term;
	}
return rez;
}

double dfxy(double x,double y){
int i;
double rez=0.0,patr,rad,term;
for(i=0;i<n;++i){
	patr=(x-v[i].x)*(x-v[i].x)+(y-v[i].y)*(y-v[i].y);
	rad=sqrt(patr);
	if(rad) term=-(x-v[i].x)*(y-v[i].y)/rad/patr;
	else term=0.0;
	rez+=term;
	}
return rez;
}

void f2x(double x,double y,double &fx1,double &fx2){
int i;
double x1,x2,patr1,patr2,rad1,rad2,term1,term2;
fx1=fx2=0;
x1=x-ex;x2=x+ex;
for(i=0;i<n;++i){
	patr1=(x1-v[i].x)*(x1-v[i].x)+(y-v[i].y)*(y-v[i].y);
	patr2=(x2-v[i].x)*(x2-v[i].x)+(y-v[i].y)*(y-v[i].y);
	if(patr1<0) patr1=-patr1;
	if(patr2<0) patr2=-patr2;
	rad1=sqrt(patr1);
	rad2=sqrt(patr2);
	if(rad1)
		term1=(x1-v[i].x)/rad1;
	else term1=0.0;
	if(rad2)
		term2=(x2-v[i].x)/rad2;
	else term2=0.0;
	fx1+=term1;fx2+=term2;
	}
}
void f2y(double x,double y,double &fy1,double &fy2){
int i;
double y1,y2,patr1,patr2,rad1,rad2,term1,term2;
fy1=fy2=0;
y1=y-ex;y2=y+ex;
for(i=0;i<n;++i){
	patr1=(y1-v[i].y)*(y1-v[i].y)+(x-v[i].x)*(x-v[i].x);
	patr2=(y2-v[i].y)*(y2-v[i].y)+(x-v[i].x)*(x-v[i].x);
	if(patr1<0) patr1=-patr1;
	if(patr2<0) patr2=-patr2;
	rad1=sqrt(patr1);
	rad2=sqrt(patr2);
	if(rad1)
		term1=(y1-v[i].y)/rad1;
	else term1=0.0;
	if(rad2)
		term2=(y2-v[i].y)/rad2;
	else term2=0.0;
	fy1+=term1;fy2+=term2;
	}
}

void poz(int li,int ls,int &piv){
int i=li,j=ls,d=0;
cmapct t;
while(i<j){
	if(w[i].d>w[j].d){
		t=w[i];w[i]=w[j];w[j]=t;d=1-d;
		}
	i+=d;j-=1-d;
	}
piv=i;
}
void qsrt(int st,int dr){
int piv;
if(st<dr){
	poz(st,dr,piv);
	qsrt(st,piv-1);
	qsrt(piv+1,dr);
	}
}

int solutie(double x,double y){
double fx1,fy1,fx2,fy2;
/*fx1=fx(x-ex,y);
fx2=fx(x+ex,y);
fy1=fy(x,y-ey);
fy2=fy(x,y+ey);
*/

f2x(x,y,fx1,fx2);
f2y(x,y,fy1,fy2);
//return 1;
return fx1*fx2<0&&fy1*fy2<0;
}

int main(){
FILE *f=fopen("adapost2.in","r");
freopen("adapost2.out","w",stdout);
fscanf(f,"%d",&n);
double sx=0.0,sy=0.0,xmed,ymed,sdmed;
double xi,xs,yi,ys,mx,my,xmax,ymax,xmin,ymin;
double dmin,_dmin,dxmax,dymax,k;
int i,j;
xmax=ymax=0.0;
xmin=ymin=10000.0;
for(i=0;i<n;++i){
	fscanf(f,"%lf%lf",&v[i].x,&v[i].y);
	sx+=v[i].x;
	sy+=v[i].y;
	if(xmax<v[i].x) xmax=v[i].x;
	if(ymax<v[i].y) ymax=v[i].y;
	if(xmin>v[i].x) xmin=v[i].x;
	if(ymin>v[i].y) ymin=v[i].y;
	}
xmed=sx/n;ymed=sy/n;
xi=xmed;yi=ymed;
if (solutie(xi,yi)) goto finish;
sdmed=dist(xi,yi);

dxmax=xmax-xmin;
dymax=ymax-ymin;
dmin=dxmax>dymax?dxmax:dymax;
k=0.2;
do{
	_dmin=dmin*k;
	nrpa=0;
	for(i=0;i<n;++i)
		if(v[i].x>=xmed-_dmin&&v[i].x<=xmed+_dmin&&
		   v[i].y>=ymed-_dmin&&v[i].y<=ymed+_dmin){
			w[nrpa].nr=i;
			w[nrpa].p=v[i];
			nrpa++;
			}
	k+=0.2;
	}while(nrpa<n&&nrpa<APR);

for(i=0;i<nrpa;++i) w[i].d=dd(w[i].p.x,w[i].p.y,xmed,ymed);
qsrt(0,nrpa-1);
if(nrpa>APR) nrpa=APR;
for(i=0;i<nrpa;++i) w[i].sd=dist(w[i].p.x,w[i].p.y);
cmapct t;
for(i=0;i<nrpa-1;++i)
	for(j=i+1;j<nrpa;++j)
		if(w[i].sd>w[j].sd){
			t=w[i];w[i]=w[j];w[j]=t;
			}
dmin=dist(w[0].p.x,w[0].p.y);
//if(dmin>sdmed) goto finish;

xi=w[0].p.x;
yi=w[0].p.y;
if(solutie(xi,yi))	goto finish;

double x2med,y2med,ddd,newd,oldd;
double _fx,_fy,_dfxy,_dfx,_dfy,delta,deltax,deltay;
int nrp;
nrp=0;
x2med=y2med=0;
for(i=0;i<3;++i) {
	x2med+=w[i].p.x;
	y2med+=w[i].p.y;
	}
x2med/=3;
y2med/=3;
xi=x2med;
yi=y2med;
ddd=0.1*dist(w[0].p.x,w[0].p.y);
xi=w[0].p.x;
yi=w[0].p.y;
oldd=w[0].sd;
do{
	//ddd=ddd*0.1;
	_fx=fx(xi,yi);
	_fy=fy(xi,yi);
	_dfxy=dfxy(xi,yi);
	_dfx=dfx(xi,yi);
	_dfy=dfy(xi,yi);
	delta=_dfx*_dfy-_dfxy*_dfxy;
	deltax=_dfy*(_fx-_dfx*xi-_dfxy*yi-ddd)-_dfxy*(_fy-_dfxy*xi-_dfy*yi-ddd);
	deltay=_dfx*(_fy-_dfxy*xi-_dfy*yi-ddd)-_dfxy*(_fx-_dfx*xi-_dfxy*yi-ddd);
	if(delta!=0){
		xs=-deltax/delta;
		ys=-deltay/delta;
		}
	else goto finish;
	newd=dist(xs,ys);
	mx=xs-xi;
	if(mx<0) mx=-mx;
	my=ys-yi;
	if(my<0) my=-my;
	nrp++;
	if(mx<=ex&&my<=ey)
		if(solutie(xi,yi)) break;
		else {mx=my=1;ddd=ddd*0.2;continue;}
	if(newd>oldd) ddd=ddd*0.2;
	else {
		ddd=newd*0.1;
		xi=xs;
		yi=ys;
		oldd=newd;
		}
/*
printf("%d   ",nrp);
printf("%0.4lf %0.4lf   ",xi,yi);
printf("%0.6lf %0.6lf %0.6lf %0.6lf %0.6lf\n",
		_fx,_fy,oldd,newd,ddd);
	   */
	}while(mx>ex||my>ey);


finish:
printf("%0.3lf %0.3lf",xi,yi);
//printf(" %0.6lf %0.6lf\n",xmed,ymed);
return 0;
}