#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>

/*Corrected version of 2-d closest points*/
/*See CLR, p. 908 and Preparata & Shamos, p. 192*/
/*Assumes distinct points*/

#define sqr(x) ((x)*(x))
#define distance(x1,y1,x2,y2) sqrt(sqr(x1-x2) + sqr(y1-y2))

float CPUtime()
{
  struct rusage rusage;
  getrusage(RUSAGE_SELF,&rusage);
  return rusage.ru_utime.tv_sec+rusage.ru_utime.tv_usec/1000000.0
    + rusage.ru_stime.tv_sec+rusage.ru_stime.tv_usec/1000000.0;
}

float elapsedTime()
{
  return (10000.0*times(NULL))/CLOCKS_PER_SEC;
}

void shellsort1(a,b,N)
float *a,*b;
int N;
/*Sorts in ascending order using table a as the key
  and table b as satellite data.  See Sedgewick for details*/
{
int group;
int i,j,h;
int v,vb;
for (h=1;h<=N/9; h=3*h+1) ;
for (;h>0; h /= 3)
  for (group=0; group<h; group++)
     for (i=group+h; i<N; i+=h)
     {
       v=a[i];
       vb=b[i];
       j=i;
       while (j>=h && a[j-h]>v)
       {
         a[j]=a[j-h];
         b[j]=b[j-h];
         j -= h;
       }
       a[j]=v;
       b[j]=vb;
     }
}


void shellsort2(a,b,myRefs,N)
float *a,*b;
int *myRefs;
int N;
/*Sorts in ascending order using table a as the key
  and table b as satellite data.  See Sedgewick for details*/
/*myRefs are subscripts for another table containing
  the same points in a different order.*/
{
int group;
int i,j,h;
int v,vb,vMyRefs;
for (h=1;h<=N/9; h=3*h+1) ;
for (;h>0; h /= 3)
  for (group=0; group<h; group++)
     for (i=group+h; i<N; i+=h)
     {
       v=a[i];
       vb=b[i];
       vMyRefs=myRefs[i];
       j=i;
       while (j>=h && a[j-h]>v)
       {    
         a[j]=a[j-h];
         b[j]=b[j-h];
         myRefs[j]=myRefs[j-h];
         j -= h;
       }
       a[j]=v;
       b[j]=vb;
       myRefs[j]=vMyRefs;
     }
}

void closestPair(Xx,Xy,Yx,Yy,Yrefs,
  numPoints,closestX1,closestY1,closestX2,closestY2,closestDist)
float *Xx,*Xy;  /*Set of points sorted ascending by x*/
float *Yx,*Yy;  /*Set of points sorted ascending by y*/
int *Yrefs;     /*Entry i points to the Xx and Xy entry that
                  is the same as the ith entry of Yx and Yy*/
int numPoints;  /*Input Size*/
float *closestX1,*closestY1,*closestX2,
  *closestY2,*closestDist; /*Used to return the closest pair*/
{
int i,j,k,closestI,closestJ;
float dist;

/*When dividing input,the X arrays are simply handled with pointers
  to the input array*/
float *Xlx,*Xly,*Xrx,*Xry;
/*Y arrays get split/copied according to natural split in X*/
float *Ylx,*Yly,*Yrx,*Yry;
int numPointsL,numPointsR;  /*Sizes for the divided input*/
/*These work just like the input "refs"*/
int *YlRefs,*YrRefs;

int lengthYl,lengthYr;  /*Used for "unmerging" the Y tables*/
/*The next two declarations are for the results of the
  recursive calls*/
float closestLx1,closestLy1,closestLx2,closestLy2,Ldistance;
float closestRx1,closestRy1,closestRx2,closestRy2,Rdistance;
/*These are work variables for the combining "sweep"*/
float lx,rx; /*x-coordinate for dividing lines*/
float lLimit,rLimit;
int rSideSubs[6],nextRpoint;
float sweepDist,divideDist;

/* dumping for the desperate
printf("entering closest\n");
for (i=0;i<numPoints;i++)
  printf("(%f %f) (%f %f)\n",Xx[i],Xy[i],Yx[i],Yy[i]);
*/

if (numPoints<=3)  /*The easy case*/
{
  closestI=0;
  closestJ=1;
  *closestDist=distance(Xx[0],Xy[0],Xx[1],Xy[1]);
  if (numPoints==3)
  {
    dist=distance(Xx[0],Xy[0],Xx[2],Xy[2]);
    if (dist< *closestDist)
    {
      *closestDist=dist;
      closestJ=2;
    }
    dist=distance(Xx[1],Xy[1],Xx[2],Xy[2]);
    if (dist< *closestDist)
    {
      *closestDist=dist;
      closestI=1;
      closestJ=2;
    }
  }
  *closestX1=Xx[closestI];
  *closestY1=Xy[closestI];
  *closestX2=Xx[closestJ];
  *closestY2=Xy[closestJ];
  return;
}

/*divide*/
numPointsL=(numPoints+1)/2;
numPointsR=numPoints/2;
/*This could be replaced by a stack-like allocation scheme*/
Ylx=(float*) malloc(numPointsL*sizeof(float));
Yly=(float*) malloc(numPointsL*sizeof(float));
Yrx=(float*) malloc(numPointsR*sizeof(float));
Yry=(float*) malloc(numPointsR*sizeof(float));
YlRefs=(int*) malloc(numPointsL*sizeof(int));
YrRefs=(int*) malloc(numPointsR*sizeof(int));

if (!Ylx || !Yly || !Yrx || !Yry || !YlRefs || !YrRefs)
{
  printf("malloc failed\n");
  exit(0);
}

/*Compute beginnings of "new" arrays*/
Xlx=Xx;
Xly=Xy;
/*Skip over left side elements (ptr arithmetic)*/
Xrx=Xx+numPointsL;
Xry=Xy+numPointsL;

/*unmerge from p. 911 of CLR*/
lengthYl=lengthYr=(-1);
for (i=0;i<numPoints;i++)
  if (Yrefs[i]<numPointsL) /*Check for left side in X*/
  {
    lengthYl++;
    Ylx[lengthYl]=Yx[i];
    Yly[lengthYl]=Yy[i];
    YlRefs[lengthYl]=Yrefs[i];
  }
  else
  {
    lengthYr++;
    Yrx[lengthYr]=Yx[i];
    Yry[lengthYr]=Yy[i];
    YrRefs[lengthYr]=Yrefs[i]-numPointsL;
  }

/*conquer*/
closestPair(Xlx,Xly,Ylx,Yly,YlRefs,numPointsL,
  &closestLx1,&closestLy1,&closestLx2,&closestLy2,&Ldistance);
closestPair(Xrx,Xry,Yrx,Yry,YrRefs,numPointsR,
  &closestRx1,&closestRy1,&closestRx2,&closestRy2,&Rdistance);
/*Use the closer pair*/
if (Ldistance<Rdistance)
{
  *closestDist=divideDist=Ldistance;
  *closestX1=closestLx1;
  *closestY1=closestLy1;
  *closestX2=closestLx2;
  *closestY2=closestLy2;
}
else
{
  *closestDist=divideDist=Rdistance;
  *closestX1=closestRx1;
  *closestY1=closestRy1;
  *closestX2=closestRx2;
  *closestY2=closestRy2;
}

/*combine - this follows Preparata & Shamos*/
lx=Xlx[numPointsL-1];
rx=Xrx[0];
/*Small improvement here*/
if (rx-lx > divideDist)
  return; /*Impossible to find a closer pair of points*/

/*Adjust coordinates for dividing lines, another improvement*/
lLimit=rx - divideDist;
rLimit=lx + divideDist;

/*Find the six right side points with smallest Y-coodinates*/
j=(-1);
for (nextRpoint=0;nextRpoint<numPointsR;nextRpoint++)
  if (Yrx[nextRpoint]<=rLimit) /*Only take points in Pr*/
  {
    j++;
    rSideSubs[j]=nextRpoint;
    if (j==5)
    {
      nextRpoint++;
      break;
    }
  }

for (j++;j<6;j++)
  rSideSubs[j]=numPointsR;  /*indicates end-of-points when <6 of them*/

/*Check each left side point with (up to) 6 right side points*/
for (i=0;i<numPointsL;i++)
  if (Ylx[i]>=lLimit)  /*Only take points in Pl*/
  {
    /*Move ahead in right side points*/
    while (rSideSubs[0]!=numPointsR
      && Yry[rSideSubs[0]]<Yly[i]-divideDist)
    {
      for (j=0;j<5;j++)
        rSideSubs[j]=rSideSubs[j+1];
      for (;nextRpoint<numPointsR;nextRpoint++)
        if (Yrx[nextRpoint]<=rLimit) /*Only take points in Pr*/
          break;
      if (nextRpoint!=numPointsR)
      {
        rSideSubs[5]=nextRpoint;
        nextRpoint++;
      }
      else
        rSideSubs[5]=numPointsR;
    }
    /*Check for improving closest points*/
    for (j=0;j<6 && rSideSubs[j]!=numPointsR;j++)
    {
      sweepDist=distance(Ylx[i],Yly[i],
        Yrx[rSideSubs[j]],Yry[rSideSubs[j]]);
      if (sweepDist< *closestDist)
      {
        *closestDist=sweepDist;
        *closestX1=Ylx[i];
        *closestY1=Yly[i];
        *closestX2=Yrx[rSideSubs[j]];
        *closestY2=Yry[rSideSubs[j]];
      }
    }
  }

free(Ylx);
free(Yly);
free(Yrx);
free(Yry);
free(YlRefs);
free(YrRefs);

/*SSSSllllooooowwwww output verification*/
/*
for (i=0;i<numPoints;i++)
  for (j=i+1;j<numPoints;j++)
  {
    sweepDist=distance(Xx[i],Xy[i],Xx[j],Xy[j]);
    if (sweepDist < *closestDist)
    {
      printf("verification problem\n");
      for (k=0;k<numPoints;k++)
        printf("(%f %f) (%f %f)\n",Xx[k],Xy[k],Yx[k],Yy[k]);
      printf("(%f,%f) (%f,%f) %f\n(%f,%f) (%f,%f) %f\n",
        Xx[i],Xy[i],Xx[j],Xy[j],sweepDist,
        *closestX1,*closestY1,*closestX2,*closestY2,*closestDist);
    }
  }
*/
}

main()
{
int i;
int numPoints;
float *Xx,*Xy;
float *Yx,*Yy;
/*Yrefs refer to same point in the other array*/
int *Yrefs;
float closestX1,closestY1,closestX2,closestY2,dist;
float startCPU,startTime;

scanf("%d",&numPoints);
Xx=(float*) malloc(numPoints*sizeof(float));
Xy=(float*) malloc(numPoints*sizeof(float));
Yx=(float*) malloc(numPoints*sizeof(float));
Yy=(float*) malloc(numPoints*sizeof(float));
Yrefs=(int*) malloc(numPoints*sizeof(int));
if (!Xx || !Xy || !Yx || !Yy || !Yrefs)
{
  printf("malloc problem\n");
  exit(0);
}
for (i=0;i<numPoints;i++)
  scanf("%f %f",&Xx[i],&Xy[i]);

startCPU=CPUtime();
startTime=elapsedTime();

shellsort1(Xx,Xy,numPoints);
for (i=0;i<numPoints;i++)
{
  Yx[i]=Xx[i];
  Yy[i]=Xy[i];
  Yrefs[i]=i;
}
shellsort2(Yy,Yx,Yrefs,numPoints);

printf("Presort CPU is %f\n",CPUtime()-startCPU);
printf("Presort elapsed time is %f\n",elapsedTime()-startTime);

closestPair(Xx,Xy,Yx,Yy,Yrefs,numPoints,&closestX1,&closestY1,
  &closestX2,&closestY2,&dist);
printf("The closest points are (%f,%f) & (%f,%f)\nat distance %f\n",
  closestX1,closestY1,closestX2,closestY2,dist);
free(Xx);
free(Xy);
free(Yx);
free(Yy);
free(Yrefs);
printf("CPU is %f\n",CPUtime()-startCPU);
printf("Elapsed time is %f\n",elapsedTime()-startTime);
}