#include <alloca.h>
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <stdlib.h>

typedef struct {
  double x;
  double y;
} point;

double sqr(double x)
{
  return x*x;
}

double dist(double x[], double y[], int i, int j) {
  return sqrt(sqr(x[i]-x[j])+
	      sqr(y[i]-y[j]));
}

double DistSqrd(point cities[], int i, int j) {
  return (sqr(cities[i].x-cities[j].x)+
	  sqr(cities[i].y-cities[j].y));
}

void swap(double *p, double *q)
{
  double tmp=*p;
  *p = *q;
  *q = tmp;
}

void tsp(point cities[], double tourx[], double toury[], int ncities)
{
  int i;
  double *ClosePt=tourx;
  double CloseDist;
  double *p;
  
  for (i=1; i<ncities; i++) {
    point c=cities[i-1];
    tourx[i]=c.x;
    toury[i]=c.y;
  }
  tourx[0] = cities[ncities-1].x;
  toury[0] = cities[ncities-1].y;
  
  for (i=1; i<ncities; i++) {
    double ThisX = tourx[i-1];
    double ThisY = toury[i-1];
    CloseDist = DBL_MAX;
    for (p=&tourx[ncities];;) {
      double ThisDist;
      do {
        p--;
        ThisDist = sqr(*p-ThisX);
      } while (ThisDist>CloseDist);
      ThisDist += sqr(toury[p-tourx]-ThisY);
      if (ThisDist <= CloseDist) {
        if (p < tourx+i)
          break;
        CloseDist = ThisDist;
        ClosePt = p;
      }
    }
    swap(&tourx[i],ClosePt);
    swap(&toury[i],&toury[ClosePt-tourx]);
  }
}

int main(int argc, char *argv[])
{
  int i, ncities;
  point *cities;
  double *tourx;
  double *toury;
  FILE *psfile;
  double sumdist = 0.0;

  if (argc!=2) {
    fprintf(stderr, "usage: %s <ncities>", argv[0]);
    exit(1);
  }
  ncities = atoi(argv[1]);
  cities = alloca(ncities*sizeof(point));
  tourx=alloca(ncities*sizeof(double));
  toury=alloca(ncities*sizeof(double));
  for (i=0; i<ncities; i++) {
    cities[i].x = ((double)(random()))/(double)(1U<<31);
    cities[i].y = ((double)(random()))/(double)(1U<<31);
  }
  tsp(cities,tourx,toury,ncities);
  psfile = fopen("tsp.eps","w");
  fprintf(psfile, "%%!PS-Adobe-2.0 EPSF-1.2\n%%%%BoundingBox: 0 0 300 300\n");
  fprintf(psfile, "1 setlinejoin\n0 setlinewidth\n");
  fprintf(psfile, "%f %f moveto\n",
	  300.0*tourx[0], 300.0*toury[0]);
  for (i=1; i<ncities; i++) {
    fprintf(psfile, "%f %f lineto\n",
	    300.0*tourx[i], 300.0*toury[i]);
    sumdist += dist(tourx, toury, i-1, i);
  }
  fprintf(psfile,"stroke\n");
  printf("sumdist = %f\n", sumdist);
  exit(0);
}