/**************************
 * Where Am I? C++ Edition
 */

#include <sstream>
#include <set>
#include <cmath>
#include <numeric>
#include <cstdlib>
#include <algorithm>
#include <iterator>
#include "WhereAmI.h"
using namespace std;

static const int MIN_DISTANCES = 2;
static const double THRESHOLD = -0.00001;

/* Checks for a precondition and reports an error otherwise. */
#define Assert(n) CheckAssertion((n), (#n))
static void CheckAssertion(bool assertion, const char *message)
{
	if(assertion)
		return;

	cerr << endl;
	cerr << "Runtime Error: Function precondition not satisfied. " << endl;
	cerr << message << endl;
	exit(-1);
}

bool circleT::operator< (const circleT &other) const
{
	if(radius < other.radius)
		return true;
	if(radius > other.radius)
		return false;
	return pt < other.pt;
}

bool pointT::operator <(const pointT &other) const
{
	if(x < other.x)
		return true;
	if(x > other.x)
		return false;
	return y < other.y;
}

/* Utility stream insertion operator for the point class. */
ostream& operator <<(ostream &out, const pointT &pt)
{
	return (out << '(' << pt.x << ", " << pt.y << ')');

}

/* Utility stream insertion operator for the circle class. */
ostream& operator <<(ostream &out, const circleT &circle)
{
	return (out << '{' << circle.pt << ", " << circle.radius << '}');
}

/* MakeCircle: Utility function that returns a new circleT by invoking the
 * class constructor.
 */
static circleT MakeCircle(double newRadius, const pointT &newPt)
{
	return circleT(newRadius, newPt);
}

/* Returns the distance between two points. */
double ComputeDistance(pointT one, pointT two)
{
	const double dx = one.x - two.x;
	const double dy = one.y - two.y;
	return sqrt(dx * dx + dy * dy);
}

/* Use getline to read from cin and return the read string. */
static string GetLine()
{
	string result;
	getline(cin, result);
	return result;
}

/* Read an unsigned integer (size_t), prompting until the user enters valid data. */
static unsigned GetPositiveInteger()
{
	while(true)
	{
		stringstream converter;
		converter << GetLine(); // Read data and load into the stream.
		
		/* Try to read an int. */
		unsigned result;
		converter >> result;
		
		/* If we read an int successfully, see if anything's left over. */
		if(!converter.fail())
		{
			/* Try to read a single char. */
			char remaining;
			converter >> remaining;
			
			/* If the stream broke, we ran out of characters and are done. */
			if(converter.fail())
			{
				return result;
			}
			else
				cout << "Unexpected character: " << remaining << endl;
		}
		else
			cout << "Please enter a real number." << endl;
		
		cout << "Retry: ";
	}
}

/* Read a double, prompting until the user enters valid data. */
static double GetReal()
{
	while(true)
	{
		stringstream converter;
		converter << GetLine(); // Read data and load into the stream.
		
		/* Try to read a double. */
		double result;
		converter >> result;
		
		/* If we read an int successfully, see if anything's left over. */
		if(!converter.fail())
		{
			/* Try to read a single char. */
			char remaining;
			converter >> remaining;
			
			/* If the stream broke, we ran out of characters and are done. */
			if(converter.fail())
				return result;
			else
				cout << "Unexpected character: " << remaining << endl;
		}
		else
			cout << "Please enter a real number." << endl;
		
		cout << "Retry: ";
	}
}

/* Read a positive double, prompting until the user enters valid data. */
static double GetPositiveReal()
{
	while(true)
	{
		stringstream converter;
		converter << GetLine(); // Read data and load into the stream.
		
		/* Try to read a double. */
		double result;
		converter >> result;
		
		/* If we read an int successfully, see if anything's left over. */
		if(!converter.fail())
		{
			/* Try to read a single char. */
			char remaining;
			converter >> remaining;
			
			/* If the stream broke, we ran out of characters and are done. */
			if(converter.fail())
			{
				if(result <= 0)
					cout << "Please enter a positive real number." << endl;
				else
					return result;
			}
			else
				cout << "Unexpected character: " << remaining << endl;
		}
		else
			cout << "Please enter a positive real number." << endl;
		
		cout << "Retry: ";
	}
}

/* Converts a set<circleT> to a vector<circleT> */
static vector<circleT> SetToVector(const set<circleT> &input)
{
	return vector<circleT>(input.begin(), input.end());
}

/* Recursively generate guesses. */
static void RecGenerateGuesses(const vector<pointT> &points, const vector<double> &distances, int index, vector<pointT> &soFar, set<set<circleT> > &result)
{
	/* If we've made a subset of the right size, generate all permutations and record them. */
	if(soFar.size() == distances.size())
	{
		do
		{
			set<circleT> thisResult;
			transform(distances.begin(), distances.end(), soFar.begin(), inserter(thisResult, thisResult.begin()), MakeCircle);
			result.insert(thisResult);			
		}
		/* This uses the next_permutation algorithm, which generates the next permuation
		 * of a data set and returns true if successful.  Once all permutations have
		 * been generated, it sorts the data and returns false.  Nifty, huh?
		 */
		while(next_permutation(soFar.begin(), soFar.end()));		
		return;
	}

	/* If we're out of guesses, abort. */
	if(index >= points.size())
		return;

	/* Otherwise, try without, then try with. */
	RecGenerateGuesses(points, distances, index + 1, soFar, result);
	
	soFar.push_back(points[index]);
	RecGenerateGuesses(points, distances, index + 1, soFar, result);
	soFar.pop_back();
}

/* Function: GenerateAllGuesses(points, distances)
 * Accepts a set of points representing the stars and a set of distances to
 * some set of stars, then returns a vector<vector<circleT > representing
 * all possible permutations of guesses.  This function assumes no duplicate
 * entries.  We assume that there are at least as many points as there are
 * distances.  We also assume that there is at least one point.
 */
vector<vector<circleT> > GenerateAllGuesses(const vector<pointT> &points, const vector<double> &distances)
{
	Assert(!points.empty());
	Assert(points.size() >= distances.size());

	/* All of the temporary outputs are stored in sets to eliminate duplicates. */
	set<set<circleT > > result;
	vector<pointT> localPoints;
	sort(localPoints.begin(), localPoints.end());

	/* "APPLY RECURSIVE ALGORITHM" to solve the problem. */
	vector<pointT> empty;
	RecGenerateGuesses(points, distances, 0, empty, result);

	/* Now, transform the set<set<..>> into a vector<vector<..>> */
	vector<vector<circleT> > vectorResult;
	transform(result.begin(), result.end(), back_inserter(vectorResult), SetToVector);

	return vectorResult;
}



/* Returns the number of distances to stars to read in. */
static unsigned GetNumDistances()
{
	while(true)
	{
		cout << "Enter the number of distances to read (" << MIN_DISTANCES << " minimum)" << endl;
		unsigned result = GetPositiveInteger();
		if(result >= MIN_DISTANCES)
			return result;

		cout << "Sorry, that's not enough distances." << endl;
	}
}

/* This function reads in a number of distances between the current point
 * and the stars.
 */
vector<double> ReadDistances()
{
	vector<double> result(GetNumDistances());
	for(int i = 0; i < result.size(); i++)
	{
		cout << "Enter distance #" << (i + 1) << ": ";
		result[i] = GetPositiveReal();
	}

	return result;
}

/* Reads in a set of points, reprompting until the user enters the specified
 * number of unique points.
 */
vector<pointT> ReadPoints(unsigned numPoints)
{
	set<pointT> result;
	for(unsigned i = 0; i < numPoints; ++i)
	{
		while(true)
		{
			pointT newPoint;
			cout << "Please enter X coordinate of point #" << (i + 1) << ": ";
			newPoint.x = GetReal();
			cout << "Please enter Y coordinate of point #" << (i + 1) << ": ";
			newPoint.y = GetReal();
	
			/* Try inserting the point, reporting an error if it's a duplicate. */
			if(result.insert(newPoint).second)
				break;

			cout << "Point already exists!" << endl;
		}
	}

	/* Construct return value using iterator range constructor. */
	return vector<pointT>(result.begin(), result.end());
}

/* Returns the number of stars we'll end up reading. */
unsigned GetNumPoints(unsigned minPoints)
{
	while(true)
	{
		cout << "Enter number of stars (" << minPoints << " required): ";
		unsigned result = GetPositiveInteger();
		if(result >= minPoints)
			return result;

		cout << "Sorry, that's not enough points." << endl;
	}
}

/* ComputeIntersectionPoints(const circleT &one, const circleT &two)
 * Returns a vector<pointT> of the intersection points of the two circles.
 * If the circles do not intersect, this function will return an empty vector.
 * If this function is called on two identical circles, it will cause a runtime
 * error, so be sure to check your input before calling this function!
 */
vector<pointT> ComputeIntersectionPoints(const circleT &one, const circleT &two)
{
	/* Cause a runtime error if the input is invalid. */
	Assert(one.pt.x != two.pt.x || one.pt.y != two.pt.y || one.radius != two.radius);

	/* We now want to apply a linear transformation to rotate the two points so
	 * that they're both on the X axis and so that c0 is at (0, 0) and c1 is at (1, 0).
	 * This transformation is represented by the matrix
	 * | dx   dy |     1
	 * |         | ---------
	 * | -dy  dx | dx^2 + dy^2
	 * Since this is a rotation/scaling operation, we must update the radii of the
	 * two circles by dividing by sqrt(x^2 + y^2).
	 */
	const double dx = two.pt.x - one.pt.x;
	const double dy = two.pt.y - one.pt.y;
	const double  r = sqrt(dx * dx + dy * dy);

	/* Circle radii after scaling. */
	const double r0 = one.radius / r;
	const double r1 = two.radius / r;

	/* Otherwise, they intersect.  The transformed X coordinate of the intersection
	 * is given by the formula below.
	 */
	const double xIntersect = (1 - r1*r1 + r0*r0) / 2.0;

	/* Now we want to find the transformed Y coordinate of the intersections.  This
	 * is given by the formula below:
	 */
	const double yIntersectSquared = r0*r0 - xIntersect*xIntersect;

	/* If this is negative, we've hit some unspecified "this is invalid" case, so
	 * bail out.
	 */
	if(yIntersectSquared < THRESHOLD)
		return vector<pointT>();

	const double yIntersect = yIntersectSquared < 0 ? 0 : sqrt(yIntersectSquared);

	/* Otherwise, the two intersection points have transformed coordinates
	 * (xIntersect, yIntersect) and (xIntersect, -yIntersect).  We need to
	 * invert the transformation by multiplying by the following matrix:
	 *     | dx  -dy |
	 * T = | dy   dx |
	 * And translating by one.x and one.y (since we assumed that the first circle was
	 * the origin).
	 */

	pointT pt1;
	pt1.x = xIntersect * dx - yIntersect * dy + one.pt.x;
	pt1.y = xIntersect * dy + yIntersect * dx + one.pt.y;

	pointT pt2;
	pt2.x = xIntersect * dx + yIntersect * dy + one.pt.x;
	pt2.y = xIntersect * dy - yIntersect * dx + one.pt.y;

	vector<pointT> result;
	result.push_back(pt1);
	result.push_back(pt2);

	return result;
}