920 – Sunny Mountains

A distance calculation geometry problem.
The illustrative figure makes the problem very simple. Sort the points according to their x co-ordinates. Iterate from left to right. Whenever a point (peak) rises higher than the previous highest peek, then the light is incident to it. What we need to do is accumulate the distance from the current peak to the bottom end of the lit part. So, the problem comes to an intersection task which is trivial given proper tools.
C++ implementation:

#include <iostream>
#include <algorithm>
#include <complex>
#include <limits>
using namespace std;

typedef long double gtype;
#define equ(a, b, e) ( abs((a) - (b)) <= (e) )
const gtype len_before_point_digits = 1E8;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_before_point_digits;
typedef complex<gtype> point;
#define x real()
#define y imag()
#define crs(a, b) ( (conj(a) * (b)).y )
#define intrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define intrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )

bool cmp(point i, point j) {
	return i.x < j.x;
}

int main() {
	int nCase, nP;
	point p[128];
	cin >> nCase;
	cout.precision(2);
	while (nCase--) {
		// input
		cin >> nP;
		for (int i = 0; i < nP; ++i)
			cin >> p[i].x >> p[i].y;
		sort(p, p + nP, cmp);
		// solve
		gtype maxY = 0, red = 0;
		for (int i = nP - 2; i >= 0; --i) {
			// a peak that sees the sun
			if (maxY <= p[i].y) {
				gtype d = intrD(p[i], p[i + 1], point(0, maxY), point(1, maxY));
				if (!equ(d, 0, epsLen)) {
					point r = intrN(p[i], p[i + 1], point(0, maxY), point(1, maxY)) / d;
					red += abs(p[i] - r);
				}
				maxY = p[i].y;
			}
		}
		cout << fixed << red << endl;
	}
	return 0;
}

438 – The Circumference of the Circle

A circumcircle problem.
Given a triangle, find the circumference of its circumcircle. There’s an easy way and a hard way.
The hard way is to find the circle’s centre then calculate the radius. I used this solution for a previous problem. Code:

#include <iostream>
#include <complex>
#include <limits>
using namespace std;

typedef long double gtype;
#define equ(a, b, e) ( abs((a) - (b)) <= (e) )
const gtype len_before_point_digits = 1E6;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_before_point_digits;
const gtype pi = M_PI;
typedef complex<gtype> point;
#define x real()
#define y imag()
#define crs(a, b) ( (conj(a) * (b)).y )
#define perp(p) ( point((p).y, -(p).x) )
#define intrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define intrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )
#define circC(r) (2 * pi * abs(r))
#define plus(x) ((x) > 0 ? " + " : " - ")

int main() {
	point p[3], v[2][2], c;
	gtype r;
	cout.precision(2);
	while (cin >> p[0].x >> p[0].y >> p[1].x >> p[1].y >> p[2].x >> p[2].y) {
		// perp 1
		v[0][0] = 0.5l * (p[0] + p[1]);
		v[0][1] = perp(p[1] - p[0]) + v[0][0];
		// perp 2
		v[1][0] = 0.5l * (p[1] + p[2]);
		v[1][1] = perp(p[2] - p[1]) + v[1][0];
		// intersect
		gtype d = intrD(v[0][0], v[0][1], v[1][0], v[1][1]);
		if (!equ(d, 0, epsLen)) {
			c = intrN(v[0][0], v[0][1], v[1][0], v[1][1]) / d;
			r = abs(c - p[0]);
		} else
			r = 0;
		// out
		cout << fixed << circC(r) << endl;
	}
	return 0;
}

The easy way is to find the circumcircle’s radius through:


And since:



This method is less prone to error. Implementation:

#include <iostream>
#include <complex>
#include <limits>
using namespace std;

typedef long double gtype;
const gtype pi = M_PI;
typedef complex<gtype> point;
#define x real()
#define y imag()
#define crs(a, b) ( (conj(a) * (b)).y )
#define circC(r) (2 * pi * abs(r))
#define area(a, b, c) ( (gtype) 0.5 * crs((b) - (a), (c) - (a)) )
#define plus(x) ((x) > 0 ? " + " : " - ")

int main() {
	point p[3];
	gtype r;
	cout.precision(2);
	while (cin >> p[0].x >> p[0].y >> p[1].x >> p[1].y >> p[2].x >> p[2].y) {
		r = abs(p[0] - p[1]) * abs(p[1] - p[2]) * abs(p[2] - p[0]) / area(p[0], p[1], p[2]) * 0.25;
		cout << fixed << circC(r) << endl;
	}
	return 0;
}

190 – Circle Through Three Points

A geometry problem.
Given three points, find the circle that passes though them. The problem statement illustrates what you have to do:
Given triangle , the circle centre is the intersection of the perpendicular bisectors on the sides. So, we need to find — the mid points of any two sides — and draw two perpendicular lines which intersect in the centre . The distance from to any point is the radius .
We need to transform this process into actual computation with vectors. A mid-point is easily determined. A perpendicular to a vector can be found by swapping the coordinates and negating only one of them. Intersection is also a simple task.






C++ implementation using some geometry tools:

#include <algorithm>
#include <iostream>
#include <complex>
#include <cstdlib>
#include <limits>
using namespace std;

typedef double gtype;
// const
const gtype len_before_point_digits = 1E6;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_before_point_digits;
// point
typedef complex<gtype> point;
#define x real()
#define y imag()
// vector
#define crs(a, b) ( (conj(a) * (b)).y )
#define perp(p) ( point((p).y, -(p).x) )
// line
#define intrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define intrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )
// I/O
#define plus(x) ((x) > 0 ? " + " : " - ")

int main() {
	point p[3], v[2][2], c;
	gtype r;
	cout.precision(3);
	while (cin >> p[0].x >> p[0].y >> p[1].x >> p[1].y >> p[2].x >> p[2].y) {
		// perp 1
		v[0][0] = 0.5 * (p[0] + p[1]);
		v[0][1] = perp(p[1] - p[0]) + v[0][0];
		// perp 2
		v[1][0] = 0.5 * (p[1] + p[2]);
		v[1][1] = perp(p[2] - p[1]) + v[1][0];
		// intersect
		c = intrN(v[0][0], v[0][1], v[1][0], v[1][1]) / intrD(v[0][0], v[0][1], v[1][0], v[1][1]);
		r = abs(c - p[0]);
		// out
		gtype h = c.x, k = c.y;
		cout << fixed << "(x" << plus(-h) << abs(h) << ")^2 + (y" <<
				plus(-k) << abs(k) << ")^2 = " << r << "^2" << endl;
		gtype a = 2 * c.x, b = 2 * c.y, d = c.x * c.x + c.y * c.y - r * r;
		cout << fixed << "x^2 + y^2" << plus(-a) << abs(a) << "x" <<
				plus(-b) << abs(b) << "y" << plus(d) << abs(d) << " = 0" << endl;
		cout << endl;
	}
	return 0;
}

C++ tools for computational geometry

A well-tested C++ library for computational geometry problems. The code has bad style and conventions, but it’ll suffice for hackathons.

Types and constants

The numerical type involved in computation varies depending on the problem. So, for flexibility we’ll name this type gtype that we’ll use throughout the code.
The macros equ and inRange should be self-explanatory, however, note that the comparison they make use some tolerance bound e.
After that, follow some constants that are of type gtype. len_before_point_digits and ang_before_point_digits are also problem dependent and they’re VERY critical. They define how many decimal digits before the decimal point we’re going to use. For example, if we know from our problem that distances will never exceed 100,000, len_before_point_digits should be 1E6. The same applies for ang_before_point_digits.
epsLen, epsAng and epsArea are the tolerance bounds we’re likely to use for distances, radian angles and areas. Yep! It’s that simple.

typedef long double gtype;
#define equ(a, b, e) ( abs((a) - (b)) <= (e) )
#define inRange(a, b, p, e) ( (p) >= min(a, b) - (e) && (p) <= max(a, b) + (e))
// CONST
const gtype len_digits_before_point = 1E6;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_digits_before_point;
const gtype ang_digits_before_point = 1E2;
const gtype epsAng = numeric_limits<gtype>::epsilon() * ang_digits_before_point;
const gtype epsArea = epsLen * epsLen;
const gtype pi = M_PI;

Points

Let’s get down to business. Let’s define a point. What’s a point? It’s an ordered pair of some type. So, instead of writing it from scratch, we’ll use C++‘s std::complex which is also an ordered pair that can be given a numerical type. std::complex spares us a lot of effort in overriding operators and has some very useful methods. We define the type std::complex <gtype> as our point.
To get a point‘s x and y, we get std::complex‘s real and imaginary parts respectively.
Another good feature of std::complex is that it allows both Cartesian and polar constructors. But, the polar constructor must be given the correct type for the magnitude component, that’s why we defined the macro.
The last two macros do simple angle conversion.

typedef complex<gtype> point;
#define x real()
#define y imag()
#define polar(r, t) polar((gtype) (r), (gtype) (t))
bool cmp_point(point const& p1, point const& p2) {
	return p1.x == p2.x ? (p1.y < p2.y) : (p1.x < p2.x);
}
#define deg(a) ( (a) * 180 /pi )
#define rad(a) ( (a) * pi /180 )

Vectors

Vectors are the building blocks of any geometrical structure. Rotations, products and transformations are very important in any problem. Fortunately, we don’t have to implement all of that. norm, abs, conj, arg are provided by std::complex. We’ll use some of these functions and the rest we’ll code on our own. I’ve chosen the fastest and most precise functions I could lay my hands on:

// VECTOR
//#define norm(a) dot(a,a) // fast for [long double]
//#define abs(a) sqrt(norm(a)) // fast for [int]
#define arg(a) ( atan2((a).y, (a).x) )
#define unit(v) ( (v) / abs(v) )
#define dot(a, b) ( (conj(a) * (b)).x )
#define crs(a, b) ( (conj(a) * (b)).y )
#define rot(v, t) ( (v) * polar(1, t) )
#define rot2(v, t, o) ( rot((v) - (o), (t)) + (o) )
#define refelct(v, m) ( conj((v) / (m)) * (m) )
#define perp(v) ( point((v).y, -(v).x) )
#define vecVecProj(u, v) ( dot(u, v) * (v) / norm(v) )

Lines

Lines-point / line-line intersections are very frequently used. So, here are some useful macros:

// LINE
#define pntLinProj(a, b, p) ( (a) + vecVecProj((p) - (a), (b) - (a)) )
#define pntLinDist(a, b, p) ( abs(crs((b)-(a), (p)-(a)) / abs((b)-(a))) )
#define linIntrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define linIntrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )
#define linOnLin(a, b, p) ( abs(crs((b) - (a), (p) - (a))) <= epsLen )
#define linOnRay(a, b, p) ( linOnLin(a, b, p) && dot((b) - (a), (p) - (a)) >= -epsLen )
#define linOnSeg(a, b, p) ( linOnRay(b, a, p) && dot((b) - (a), (p) - (a)) >= -epsLen )
#define linInRect(a, b, p, e) ( inRange((a).x, (b).x, (p).x, e) && inRange((a).y, (b).y, (p).y, e) )

Triangles

Perimeter, surface, and cosine law:

// TRIANGLE
#define triPerim(a, b, c) ( (a) + (b) + (c) )
#define triSurface(a, b, c) ( (gtype) 0.5 * crs((b) - (a), (c) - (a)) )
#define triSurfaceH(s, a, b, c) ( sqrt((s) * ((s) - (a)) * ((s) - (b)) * ((s) - (c))) )
#define cosLaw(a, b, c) ( ((a) * (a) + (b) * (b) - (c) * (c)) / (2 * (a) * (b)) )

Circles

Arc length, chord length, sector surface area, and circle circle intersection:

// CIRCLE
#define cirArcLen(r, a) ((a) * abs(r))
#define cirSectorSurface(r, a) ((a) * norm(r) / 2)
#define cirChordLen(r, a) sqrt(2 * norm(r) * (1 - cos(a)))
int cirCircInter(point & c0, point & c1, gtype r0, gtype r1, vector<point> & out) {
	point cT, v;
	gtype rT, phi;
	// init
	if (r1 > r0)
		rT = r0, r0 = r1, r1 = rT, cT = c0, c0 = point(c1), c1 = point(cT);
	out.clear(), v = c1 - c0, phi = acos(cosLaw(r0, abs(v), r1));
	// same center
	if (equ(abs(v), 0, epsLen))
		equ(r0, r1, epsLen) ? (equ(r0, 0, epsLen) ? out.push_back(c0) : out.resize(3)) : out.clear();
	// none
	else if (abs(abs(v) - r0) > r1 + epsLen)
		out.clear();
	// one point
	else if (equ(abs(abs(v) - r0), r1, epsLen))
		out.push_back(unit(v) * r0 + c0);
	// two points
	else
		out.push_back(rot(unit(v) * r0, -phi) + c0), out.push_back(rot(unit(v) * r0, phi) + c0);
	return out.size();
}

Polygons

Area, convexity check, and rotating calipers (min-width). Also includes convex hull creation:

// POLYGON (must be closed, start & end vertex overlap)
#define polAngle(n) (((n) - 2) * pi / (n))
long double polArea(vector<point> & p) {
	long double area = 0;
	for (size_t i = 0; i < p.size() - 1; ++i)
		area += crs(p[i], p[i + 1]);
	return abs(area * 0.5l);
}
bool polIsConvex(vector<point> & p) {
	int a[2] = { 0, 0 }, n = p.size() - 1;
	for (size_t i = 1; i < p.size(); ++i) {
		long double cross = crs(p[i % n] - p[i - 1], p[(i + 1) % n] - p[i - 1]);
		if (cross != 0)
			++a[cross > 0];
	}
	return !(a[0] && a[1]);
}
// CONVEX HULL
// O(n.h) - jarvis's march
// TODO
// O(n.log(n)) - monotone chain
vector<point> mcH;
void monotoneChain(vector<point> &ps) {
	vector<point> p(ps.begin(), ps.end() - 1);
	int n = p.size(), k = 0;
	mcH = vector<point>(2 * n);
	sort(p.begin(), p.end(), cmp_point);
	for (int i = 0; i < n; i++) {
		while (k >= 2 && crs(mcH[k - 1] - mcH[k - 2], p[i] - mcH[k - 2]) <= 0)
			k--;
		mcH[k++] = p[i];
	}
	for (int i = n - 2, t = k + 1; i >= 0; i--) {
		while (k >= t && crs(mcH[k - 1] - mcH[k - 2], p[i] - mcH[k - 2]) <= 0)
			k--;
		mcH[k++] = p[i];
	}
	mcH.resize(k);
}
// O(n) - rotating calipers (works on a ccw closed convex hull)
gtype rotatingCalipers(vector<point> &ps) {
	int aI = 0, bI = 0;
	for (size_t i = 1; i < ps.size(); ++i)
		aI = (ps[i].y < ps[aI].y ? i : aI), bI = (ps[i].y > ps[bI].y ? i : bI);
	gtype minWidth = ps[bI].y - ps[aI].y, aAng, bAng;
	point aV = point(1, 0), bV = point(-1, 0);
	for (gtype ang = 0; ang < pi; ang += min(aAng, bAng)) {
		aAng = acos(dot(ps[aI + 1] - ps[aI], aV) / abs(aV) / abs(ps[aI + 1] - ps[aI]));
		bAng = acos(dot(ps[bI + 1] - ps[bI], bV) / abs(bV) / abs(ps[bI + 1] - ps[bI]));
		aV = rot(aV, min(aAng, bAng)), bV = rot(bV, min(aAng, bAng));
		if (aAng < bAng)
			minWidth = min(minWidth, pntLinDist(ps[aI], ps[aI] + aV, ps[bI])), aI = (aI + 1) % (ps.size() - 1);
		else
			minWidth = min(minWidth, pntLinDist(ps[bI], ps[bI] + bV, ps[aI])), bI = (bI + 1) % (ps.size() - 1);
	}
	return minWidth;
}