typedef struct {
	double	x;
	double	y;
	int		group;
} POINT;



POINT * gen_xy(int num_pts, double radius)
{
	int		i;
	double	ang, r;
	POINT * pts;

	pts = (POINT*) malloc(sizeof(POINT) * num_pts);

	for ( i = 0; i < num_pts; i++ ) {
		ang = 2.0 * M_PI * rand() / (RAND_MAX - 1.);
		r = radius * rand() / (RAND_MAX - 1.);
		pts[i].x = r * cos(ang);
		pts[i].y = r * sin(ang);
	}
	return pts;	
}


double dist2(POINT * a, POINT * b)
{
	double x = a->x - b->x;
	double y = a->y - b->y;
	return x*x + y*y;
}


int nearest(POINT * pt, POINT * cent, int n_cluster)
{
	int i, clusterIndex;
	double d, min_d;

	min_d = HUGE_VAL;
	clusterIndex = pt->group;	
	for (i = 0; i < n_cluster; i++) {
		d = dist2(&cent[i], pt);
		if ( d < min_d ) {
			min_d = d;
			clusterIndex = i;
		}
	}	
	return clusterIndex;
}


double nearestDistance(POINT * pt, POINT * cent, int n_cluster)
{
	int i;
	double d, min_d;

	min_d = HUGE_VAL;
	for (i = 0; i < n_cluster; i++) {
		d = dist2(&cent[i], pt);
		if ( d < min_d ) {
			min_d = d;
		}
	}

	return min_d;
}


int bisectionSearch(double *x, int n, double v)
{
	int il, ir, i;


	if (n < 1) {
		return 0;
	}
	
	if (v < x[0]) {
		return 0;
	}
	else if (v > x[n-1]) {
		return n - 1;
	}

	
	il = 0;
	ir = n - 1;

	i = (il + ir) / 2;
	while ( i != il ) {
		if (x[i] <= v) {
			il = i;
		} else {
			ir = i;
		}
		i = (il + ir) / 2;		
	}		

	if (x[i] <= v)
		i = ir;
	return i;
} 


void kppAllinger(POINT * pts, int num_pts, POINT * centroids,
		 int num_clusters)
{
	int j;
	int selectedIndex;
	int cluster;
	double sum;
	double d;
	double random;	
	double * cumulativeDistances;
	double * shortestDistance;


	cumulativeDistances = (double*) malloc(sizeof(double) * num_pts);
	shortestDistance = (double*) malloc(sizeof(double) * num_pts);	


	
	selectedIndex = rand() % num_pts;
	centroids[0] = pts[ selectedIndex ];

	for (j = 0; j < num_pts; ++j)
		shortestDistance[j] = HUGE_VAL;	

	
	for (cluster = 1; cluster < num_clusters; cluster++) {

		
		for ( j = 0; j < num_pts; j++ ) {
			d = dist2(&pts[j], &centroids[cluster-1] );

			if (d < shortestDistance[j])
				shortestDistance[j] = d;
		}

		
		sum = 0.0;
		for (j = 0; j < num_pts; j++) {
			sum += shortestDistance[j];
			cumulativeDistances[j] = sum;
		}

		
		random = (float) rand() / (float) RAND_MAX * sum;
		selectedIndex = bisectionSearch(cumulativeDistances, num_pts, random);

		
		centroids[cluster] = pts[selectedIndex];
	}

	
	for (j = 0; j < num_pts; j++)
		pts[j].group = nearest(&pts[j], centroids, num_clusters);

	free(shortestDistance);
	free(cumulativeDistances);

	return;
}	


void kpp(POINT * pts, int num_pts, POINT * centroids,
		 int num_clusters)
{
	int j;
	int cluster;
	double sum;
	double * distances;


	distances = (double*) malloc(sizeof(double) * num_pts);

	
	centroids[0] = pts[ rand() % num_pts ];


	
	for (cluster = 1; cluster < num_clusters; cluster++) {

		
		sum = 0.0;
		for ( j = 0; j < num_pts; j++ ) {
			distances[j] = 
				nearestDistance(&pts[j], centroids, cluster);
			sum += distances[j];
		}

		
		sum = sum * rand() / (RAND_MAX - 1);

		
		for (j = 0; j < num_pts; j++ ) {
			sum -= distances[j];
			if ( sum <= 0)
			{
				centroids[cluster] = pts[j];
				break;
			}
		}
	}

	
	for (j = 0; j < num_pts; j++)
		pts[j].group = nearest(&pts[j], centroids, num_clusters);

	free(distances);

	return;
}	



POINT * lloyd(POINT * pts, int num_pts, int num_clusters, int maxTimes)
{
	int i, clusterIndex;
	int changes;
	int acceptable = num_pts / 1000;	


	if (num_clusters == 1 || num_pts <= 0 || num_clusters > num_pts )
		return 0;


	POINT * centroids = (POINT *)malloc(sizeof(POINT) * num_clusters);

	if ( maxTimes < 1 )
		maxTimes = 1;



	


	
	kppAllinger(pts, num_pts, centroids, num_clusters);

	do {
		

		
		for ( i = 0; i < num_clusters; i++ ) {
			centroids[i].group = 0;		
			centroids[i].x = 0;			
			centroids[i].y = 0;			
		}

		
		for (i = 0; i < num_pts; i++) {
			clusterIndex = pts[i].group;
			centroids[clusterIndex].group++;
			centroids[clusterIndex].x += pts[i].x;
			centroids[clusterIndex].y += pts[i].y;
		}

		
		for ( i = 0; i < num_clusters; i++ ) {
			centroids[i].x /= centroids[i].group;
			centroids[i].y /= centroids[i].group;
		}

		
		changes = 0;
		for ( i = 0; i < num_pts; i++ ) {
			clusterIndex = nearest(&pts[i], centroids, num_clusters);
			if (clusterIndex != pts[i].group) {
				pts[i].group = clusterIndex;
				changes++;
			}
		}

		maxTimes--;
	} while ((changes > acceptable) && (maxTimes > 0));

	
	for ( i = 0; i < num_clusters; i++ )
		centroids[i].group = i;

	return centroids;
}	


void print_eps(POINT * pts, int num_pts, POINT * centroids, int num_clusters)
{



	int i, j;
	double min_x, max_x, min_y, max_y, scale, cx, cy;
	double *colors = (double *) malloc(sizeof(double) * num_clusters * 3);

	for (i = 0; i < num_clusters; i++) {
		colors[3*i + 0] = (3 * (i + 1) % 11)/11.;
		colors[3*i + 1] = (7 * i % 11)/11.;
		colors[3*i + 2] = (9 * i % 11)/11.;
	}

	max_x = max_y = - HUGE_VAL;
	min_x = min_y = HUGE_VAL;
	for (j = 0; j < num_pts; j++) {
		if (max_x < pts[j].x) max_x = pts[j].x;
		if (min_x > pts[j].x) min_x = pts[j].x;
		if (max_y < pts[j].y) max_y = pts[j].y;
		if (min_y > pts[j].y) min_y = pts[j].y;
	}

	scale = W / (max_x - min_x);
	if (scale > H / (max_y - min_y))
		scale = H / (max_y - min_y);
	cx = (max_x + min_x) / 2;
	cy = (max_y + min_y) / 2;

	printf(, W + 10, H + 10);
	printf( 
		
		
		
		
	);


	for (i = 0; i < num_clusters; i++) {
		printf(,
			colors[3*i], colors[3*i + 1], colors[3*i + 2]);

		for (j = 0; j < num_pts; j++) {
			if (pts[j].group != i) continue;
			printf(,
				(pts[j].x - cx) * scale + W / 2,
				(pts[j].y - cy) * scale + H / 2);
		}
		printf(,
			(centroids[i].x - cx) * scale + W / 2,
			(centroids[i].y - cy) * scale + H / 2);
	}
	printf();

	free(colors);

	return;
}	


int main()
{
	int		num_pts = NUMBER_OF_POINTS;
	int		num_clusters = NUMBER_OF_CLUSTERS;
	int		maxTimes = MAXIMUM_ITERATIONS;
	double	radius = RADIUS;
	POINT * pts;
	POINT * centroids;

	
	pts = gen_xy(num_pts, radius);

	
	centroids = lloyd(pts, num_pts, num_clusters, maxTimes);

	
	print_eps(pts, num_pts, centroids, num_clusters);

	free(pts);
	free(centroids);

	return 0;
}