#include <random>
#include <cmath>
#include <queue>
#include <iostream>

using namespace std;
using ld = long double;

struct pt {
	ld x = 0, y = 0;
};

constexpr pt zero = pt{0, 0};

pt operator+(pt a, pt b) {
	a.x += b.x;
	a.y += b.y;
	return a;
}

pt operator-(pt a, pt b) {
	a.x -= b.x;
	a.y -= b.y;
	return a;
}

pt operator*(ld s, pt p) {
	p.x *= s;
	p.y *= s;
	return p;
}

ld dot(pt a, pt b) {
	return a.x * b.x + a.y * b.y;
}

ld normsqr(pt a) {
	return dot(a, a);
}

ld lensqr(pt a, pt b) {
	return normsqr(b - a);
}

ld norm(pt a) {
	return sqrt(normsqr(a));
}

struct cone {
	pt p1, d1;
	pt p2, d2;
};

struct cone_d {
	ld d;
	cone c;

	cone_d(cone const & co) : c(co) {
		// TODO: proper point line segment distance
		d = min(normsqr(c.p1), normsqr(c.p2));
	}
};

bool operator<(cone_d const & lhs, cone_d const & rhs) {
	return lhs.d > rhs.d;
}

// two points and a direction from p2
// returns x such that 
ld unit_triangle(pt p1, pt p2, pt d2) {
	const auto d1 = p1 - p2;
	// from https://en.wikipedia.org/wiki/Triangle#Using_coordinates
	// 1 = area = 1/2 * abs(x * d2.x * d1.y - x * d1.x * d2.y)
	// 2 = x * abs(d2.x * d1.y - d1.x * d2.y)
	// 2 / abs(d2.x * d1.y - d1.x * d2.y) = x
	const auto x = 2.0 / abs(d2.x * d1.y - d1.x * d2.y);
	return x;
}

int main() {
	mt19937 gen(random_device{}());
	
	// This is used for the starting triangles
	uniform_real_distribution<ld> dist(0.8, 1.6);

	uniform_real_distribution<ld> unif01(0, 1);

	// the three starting rays
	pt c1 = {0, 1};
	pt c2 = { sqrt(3) / 2.0, -0.5};
	pt c3 = {-sqrt(3) / 2.0, -0.5};

	// and the random triangles attached to them
	pt p12a = dist(gen) * c1;
	pt p12b = unit_triangle(p12a, {0, 0}, c2) * c2;
	cone c12{p12a, c1, p12b, c2};

	pt p23a = dist(gen) * c2;
	pt p23b = unit_triangle(p23a, {0, 0}, c3) * c3;
	cone c23{p23a, c2, p23b, c3};

	pt p31a = dist(gen) * c3;
	pt p31b = unit_triangle(p31a, {0, 0}, c1) * c1;
	cone c31{p31a, c3, p31b, c1};

	// we keep track of all lines in the construction
	vector<pair<pt, pt>> rays;
	rays.emplace_back(zero, c1);
	rays.emplace_back(zero, c2);
	rays.emplace_back(zero, c3);

	vector<pair<pt, pt>> lines;
	lines.emplace_back(p12a, p12b);
	lines.emplace_back(p23a, p23b);
	lines.emplace_back(p31a, p31b);

	// I use a priority queue to pick the closest cone
	priority_queue<cone_d> work;
	work.push(c12);
	work.push(c23);
	work.push(c31);

	// we pick the closest cone to expand
	while(!work.empty() && lines.size() < 2000) {
		auto c = work.top().c;
		work.pop();

		// one quadrant is enough for me
		if((c.p1.x < -4 && c.p2.x < -4) || (c.p1.y < -4 && c.p2.y < -4))
			continue;

		// split or not?
		if (lensqr(c.p1, c.p2) > 16) {
			// choose midpoint on c.p2 - c.p1
			// such that both are length >= 2
			// choose direction such that we get two cones
			// The two pow(-) functions skew the distribution,
			// you can artistically choose something here.
			const auto l = sqrt(lensqr(c.p1, c.p2));
			const auto wiggle = l - 4;
			const auto r0 = unif01(gen);
			const auto r = 2.0 / l + pow(r0, 1.5) * wiggle / l;
			const auto pm = (1.0 - r) * c.p1 + r * c.p2;

			const auto r20 = unif01(gen);
			const auto r2 = 0.001 + 0.998 * pow(r20, 1.2);
			const auto dm = r2 * c.d1 + (1.0 - r2) * c.d2;

			rays.emplace_back(pm, dm);

			cone c1{c.p1, c.d1, pm, dm};
			cone c2{pm, dm, c.p2, c.d2};
			work.push(c1);
			work.push(c2);
		} else {
			// try the two triangles and see which one fits the bill		
			const auto x1 = unit_triangle(c.p2, c.p1, c.d1);
			const auto x2 = unit_triangle(c.p1, c.p2, c.d2);

			const auto p1b = c.p1 + x1 * c.d1;
			const auto p2b = c.p2 + x2 * c.d2;
			
			// Here I do not follow the paper correctly. I should pick the
			// new cone with angles at most bla bla bla. Instead, I choose
			// the one with a shorter edge, sounds fine to me. I reintroduce
			// some randomness for artistic reasons. Increase 20 for wild
			// results.
			if (lensqr(c.p1, p2b) + 20*unif01(gen) < lensqr(p1b, c.p2) + 20*unif01(gen)) {
				cone c2{c.p1, c.d1, p2b, c.d2};
				work.push(c2);
				lines.emplace_back(c.p1, p2b);
			} else {
				cone c2{p1b, c.d1, c.p2, c.d2};
				work.push(c2);
				lines.emplace_back(p1b, c.p2);
			}
		}
	}

	// the print an svg attribute
	const auto p = [](auto str, auto v) {
		cout << ' ' << str << "=\"" << v << "\"";
	};

	// I roughly estimate the maximal distance from zero
	ld max_l = 0;

	shuffle(lines.begin(), lines.end(), gen);

	cout << "<svg>\n";
	for(auto && l : lines) {
		max_l = 0.5 * max_l + 0.5 * max(max_l, max(norm(l.first), norm(l.second)));
		cout << "<line";
		p("x1", l.first.x);
		p("y1", l.first.y);
		p("x2", l.second.x);
		p("y2", l.second.y);
		p("style", "stroke:#000;stroke-width:0.1");
		cout << "/>\n";
	}
	for(auto && l : rays) {
		const auto p1 = l.first;
		const auto l_left = max(ld(1.0), max_l - norm(p1));
		const auto p2 = p1 + (l_left / norm(l.second)) * l.second;
		cout << "<line";
		p("x1", p1.x);
		p("y1", p1.y);
		p("x2", p2.x);
		p("y2", p2.y);
		p("style", "stroke:#000;stroke-width:0.1");
		cout << "/>\n";
	}
	cout << "</svg>" << endl;
}