#include <stdio.h>
#include <math.h>

class point
{
public:
	double x,y;
	point() {}
	point( double nx, double ny ) : x(nx),y(ny) {}
	point  operator+ ( const point & p ) const { return point(x+p.x,y+p.y); }
	point  operator- ( const point & p ) const { return point(x-p.x,y-p.y); }
	point  operator* ( const double s  ) const { return point(x*s,y*s); }
	double operator* ( const point & p ) const { return x*p.x+p.y*y; }	// dot   product
	double operator^ ( const point & p ) const { return x*p.y-p.x*y; }	// cross product
	double angle() const { return atan2(y,x); }
	point& at_angle( double a ) { x = cos(a); y = sin(a); return *this; }
	point  perp() const { return point(y,-x); }
	double norm () const { return sqrt(x*x+y*y); }
	point& normalize() { double n = norm(); if(n!=0) { x/=n; y/=n; } return *this; }
};

class line
{
public:
	point orig;
	point dire;
	line() {}
	line( const point & no, const point & nd ) : orig(no),dire(nd) {}
	point param( double t ) const { return orig+dire*t; }
	point intersect( const line & l ) { return param( (l.dire^(orig-l.orig))/(dire^l.dire) ); }
};

void main()
{
	const double SX = 800; const double SY = 800;
	const double S  = 700; const double Q = 75;
	const int N = 3;
	int i;

	point tri[N] = { point(80,50) ,point(680,750), point(165,715) };
	line  median   [N];
	line  bisector [N];
	line  symmedian[N];
	point inter[N][2];

	for(i=0;i<N;++i) {
		median[i].orig = tri[i];
		median[i].dire = ((tri[(i+1)%N]+tri[(i+2)%N])*0.5 - tri[i]).normalize();

		bisector[i].orig = tri[i];
		bisector[i].dire = (tri[(i+1)%N]-tri[i]).normalize() + (tri[(i+2)%N]-tri[i]).normalize();
		bisector[i].dire.normalize();

		symmedian[i].orig = tri[i];
		symmedian[i].dire.at_angle( bisector[i].dire.angle()*2 - median[i].dire.angle() );
	}

	point symmedian_p = symmedian[0].intersect(symmedian[1]);

	for(i=0;i<N;++i) {
		line pa(symmedian_p, (tri[(i+2)%N]-tri[(i+1)%N]).normalize() );
		for(int j=0;j<2;++j) {
			line la(tri[i], (tri[(i+1+j)%N]-tri[i]).normalize() );
			inter[i][j] = pa.intersect(la);
		}
	}

	line l1( (inter[0][0]+inter[0][1])*0.5, (inter[0][0]-inter[0][1]).perp().normalize() );
	line l2( (inter[1][0]+inter[1][1])*0.5, (inter[1][0]-inter[1][1]).perp().normalize() );
	point  first_lem_c = l1.intersect(l2);
	double first_lem_r = (first_lem_c-inter[0][0]).norm();

	FILE * fp = fopen("c:\\temp\\Lemoine_Hexagon.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=\"Lemoine_Hexagon\">\n"
		,SX,SY
	);

	fprintf(fp,"<g style=\"stroke:#0000C0;stroke-width:1;stroke-opacity:1;stroke-dasharray:6,4;fill:none\">\n");
	for(i=0;i<N;++i)
		fprintf(fp,
			"<line x1=\"%5.1lf\" y1=\"%5.1lf\" x2=\"%5.1lf\" y2=\"%5.1lf\"/>\n"
			,median[i].orig.x,median[i].orig.y
			,(median[i].orig+median[i].dire*S).x,(median[i].orig+median[i].dire*S).y
		);
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:#00C000;stroke-width:1;stroke-opacity:1;stroke-dasharray:6,4;fill:none\">\n");
	for(i=0;i<N;++i)
		fprintf(fp,
			"<line x1=\"%5.1lf\" y1=\"%5.1lf\" x2=\"%5.1lf\" y2=\"%5.1lf\"/>\n"
			,bisector[i].orig.x,bisector[i].orig.y
			,(bisector[i].orig+bisector[i].dire*S).x,(bisector[i].orig+bisector[i].dire*S).y
		);
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:#C00000;stroke-width:1;stroke-opacity:1;stroke-dasharray:6,4;fill:none\">\n");
	for(i=0;i<N;++i)
		fprintf(fp,
			"<line x1=\"%5.1lf\" y1=\"%5.1lf\" x2=\"%5.1lf\" y2=\"%5.1lf\"/>\n"
			,symmedian[i].orig.x,symmedian[i].orig.y
			,(symmedian[i].orig+symmedian[i].dire*S).x,(symmedian[i].orig+symmedian[i].dire*S).y
		);
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:#C0C0C0;stroke-width:1;stroke-opacity:1;stroke-dasharray:6,4;fill:none\">\n");
	for(i=0;i<N;++i) {
		point d = (inter[i][1] - inter[i][0]).normalize();
		point p1 = inter[i][0] - d*Q; point p2 = inter[i][1] + d*Q;
		fprintf(fp,
			"<line x1=\"%5.1lf\" y1=\"%5.1lf\" x2=\"%5.1lf\" y2=\"%5.1lf\"/>\n"
			,p1.x, p1.y, p2.x, p2.y
		);
	}
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:#00C0C0;stroke-width:1;stroke-opacity:1;fill:none\">\n");
	fprintf(fp,"<circle cx=\"%5.1lf\" cy=\"%5.1lf\" r=\"%5.1lf\"/>\n"
		,first_lem_c.x, first_lem_c.y
		,first_lem_r
	);
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:#C000C0;stroke-width:2;stroke-opacity:1;fill:none\">\n");
	fprintf(fp,"<path d=\"");
	for(i=0;i<3;++i){
		if(i==0) fprintf(fp,"M ");
		else     fprintf(fp,"L ");
		fprintf(fp,"%5.1lf,%5.1lf ",tri[i].x,tri[i].y);
	}
	fprintf(fp,"z\"/>\n");
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:#000;stroke-width:2;stroke-opacity:1;fill:#C0C000;fill-opacity:0.5\">\n");
	fprintf(fp,"<path d=\"");
	for(i=0;i<3;++i) {
		for(int j=0;j<2;++j) {
			if(i==0 && j==0) fprintf(fp,"M ");
			else             fprintf(fp,"L ");
			fprintf(fp,"%5.1lf,%5.1lf ",inter[i][j]);
		}
	}
	fprintf(fp,"z\"/>\n");
	fprintf(fp,"</g>\n");

	fprintf(fp,"<g style=\"stroke:none;fill:#000\">\n");
	fprintf(fp,"<circle cx=\"%5.1lf\" cy=\"%5.1lf\" r=\"%5.1lf\"/>\n"
		,symmedian_p.x, symmedian_p.y
		,6.0
	);
	for(i=0;i<N;++i) for(int j=0;j<2;++j)
			fprintf(fp,"<circle cx=\"%5.1lf\" cy=\"%5.1lf\" r=\"%5.1lf\"/>\n"
				,inter[i][j].x, inter[i][j].y
				,6.0
			);
	fprintf(fp,"</g>\n");

	fprintf(fp,"</svg>\n");
	fclose(fp);
}