// Making the 600-cell petrie polygon
// (C)2008 Claudio Rocchini

bool even_permutation( int p[], int n )
{
	int i,ns;
	int * pt = new int[n];
	for(i=0;i<n;++i) pt[i] = p[i];

	ns = 0;
	for(;;) {
		int j,k;
		for(j=0;j<n;++j) if(pt[j]!=j) break;
		if(j==n) break;
		for(k=0;k<n;++k) if(pt[k]==j) break;
		std::swap(pt[k],pt[j]);
		++ns;
	}
	delete[] pt;
	return ns%2==0;
}

void main() {
	const double EPS = 1e-6;
	const int ND = 4; const int NV = 120; const int NE = 720;
	double v[NV][ND]; int e[NE][2];
	double px[NV]; double py[NV];
	
		// Magics! Hard to find projection directions
	double PX[ND] = {  0.33086939364114176, 0.66913060635885813,  0.66913060635885824, 0.33086939364114182 };
	double PY[ND] = { -0.74314482547739424, 0.07032846067420610, -0.07032846067420611, 0.74314482547739424 };

	const double fi = (1.0+sqrt(5.0))/2.0;
	const double le = 1.0/fi;
	
	int i,j,k,l;

	k = 0;
	for(i=0;i<16;++i) {
		for(j=0;j<ND;++j)
			v[k][j] = -0.5 + ( i&(1<<j) ? 0 : 1 );
		++k;
	}
	for(i=0;i<4;++i) for(j=0;j<2;++j) {
		v[k][0] = v[k][1] = v[k][2] = v[k][3] = 0;
		v[k][i] = -1+2*j;
		++k;
	}
	for(j=0;j<8;++j) {
		double val[4];
		int    per[4];

		val[0] = 0;
		val[1] = j&(1<<0) ? -0.5  : 0.5;
		val[2] = j&(1<<1) ? -fi/2 : fi/2;
		val[3] = j&(1<<2) ? -le/2 : le/2;

		for(i=0;i<4;++i) per[i] = i;
		do {
			if(even_permutation(per,4)) {
				for(i=0;i<4;++i) v[k][i] = val[per[i]];
				++k;
			}
		} while(std::next_permutation(per,per+4));
	}
	assert(k==NV);

	k = 0;
	for(i=  0;i<NV-1;++i) for(j=i+1;j<NV  ;++j) {
		double d = 0;
		for(l=0;l<ND;++l) d += (v[i][l]-v[j][l])*(v[i][l]-v[j][l]);
		d = sqrt(d);
		if( fabs(le-d)<EPS ) { e[k][0] = i; e[k][1] = j; ++k; }
	}
	assert(k==NE);

	for(i=0;i<NV;++i){
		px[i] = 0; for(l=0;l<ND;++l) px[i] += v[i][l]*PX[l];
		py[i] = 0; for(l=0;l<ND;++l) py[i] += v[i][l]*PY[l];
	}
	
	const double SX = 800; const double SY = 800;
	const double B  = 32;  const double R  = 6;
	
	const double sca  = std::min((SX-2*B)/2,(SY-2*B)/2);
	for(i=0;i<NV;++i) { px[i] = B+(px[i]+1)*sca; py[i] = B+(py[i]+1)*sca; }

	FILE * fp = fopen("c:\\temp\\New600Cell.svg","w");
	fprintf(fp,
		"<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n"
		"<svg\n"
		"xmlns:svg=\"http://www.w3.org/2000/svg\"\n"
		"xmlns=\"http://www.w3.org/2000/svg\"\n"
		"version=\"1.0\"\n"
		"width=\"%g\"\n" "height=\"%g\"\n"
		"id=\"New600Cell\">\n"
		,SX,SY
	);
	fprintf(fp,"<g style=\"stroke:#000000;stroke-width:2;stroke-opacity:0.5;\">\n");
	for(i=0;i<NE;++i)
		fprintf(fp,
			"<line x1=\"%5.1lf\" y1=\"%5.1lf\" x2=\"%5.1lf\" y2=\"%5.1lf\"/>\n"
			,px[e[i][0]],py[e[i][0]], px[e[i][1]],py[e[i][1]]
		);
	fprintf(fp,"</g>\n");
	fprintf(fp,"<g style=\"stroke:#000000;stroke-width:2;stroke-opacity:0.5;fill:#0080FF\">\n");
	for(i=0;i<NV;++i)
		fprintf(fp,"<circle cx=\"%5.1lf\" cy=\"%5.1lf\" r=\"%5.1lf\"/>\n",px[i],py[i],R);
	fprintf(fp,"</g>\n");
	fprintf(fp,"</svg>\n");
	fclose(fp);
}