CIRKILL - Non-official editorial describing my approach

Well, I shall begin by say that this was by far, the most interesting problem I have managed to solve so far on Codechef and I loved everyday of work on it!! And as I liked it so much, I feel like I should describe my solution here and state some things that are now fully clear. (Note: this is a non-official solution nor tutorial and my solution might be very more hard-working than the presented solution, but, this was my geometric approach using divisions and written before editorial was published.)

  • Formulating the algorithmic approach and program setup

This problem’s constrains were low enough such that a brute-force geometric approach was possible. That is because we can have as maximum number of points, only 30, and as it was stated on problem statement that in terms of the ordering of TR they are indistinguishable between them, this essentially means that there are 30C3 maximum arrangements to consider, when nCr stands for N choose R. It is important to mention that on this problem R always equals 3, since we are dealing with triangles only (degenerate or non-degenerate, they are still triangles).

NC3 is quite small number (however, calculating it directly from its formula would be complicated), and using properties of factorials we can reduce it to (N(N-1)(N-2))/6, which will be the number of triangles (again, degenerate or not, they are still called triangles on my tutorial) in a set of N points.

As we understood from beginning, since we are dealing with geometric objects, it’s useful to have suitable representations of these entities which can be easily manipulated, i.e., we want to be able to compute distance between points, areas of triangles, etc. so it helps if we abstract these concepts to a higher level (yes, this is Object-Oriented Programming - lesson 0 eheh) so what we will do is to represent both the concept of point and triangle as new data structures, such that we can write in C++:

struct point
{
	double x,y;
};

struct triangle
{
	point a,b,c;
}; 

Note here that, despite the coordinates of the points being integers, I am defining them as double already, as I am aware that all precision will be very important and key on solving this problem efficiently.

  • Reading the input data and filling in the necessary data structures

As we have seen on program setup section, we are now ready to create and fill in our necessary data structures, which is where the brute-force approach is more noticeable.

For each test case, we simply read N pairs of integers (which notice, are read as double data type), or points, from stdin and place them on a vector of points with size N.

Then, we create vector of triangles with size (N(N-1)(N-2))/6 and we use triple nested for loop to fill this vector.

Note that both vectors are declared & initialized, which essentially makes them work as arrays here. The reason why they are initialized, is mostly because I knew their size beforehand and growing vectors during runtime can be quite expensive in terms of time comsumption.

In C++, this can be stated as:

int t;
	scanf("%d",&t);
	while(t--)
	{
		int n;
		scanf("%d",&n);
		vector<point> v(n);
		vector<triangle> t((n*(n-1)*(n-2))/6);
		for(int np = 0; np < n; np++)
		{
			point p;
			scanf("%lf %lf",&p.x, &p.y);
			v[np] = p;
		}
		int idx=0;
		for(int i = 0; i < n-2; i++)
		{
			for(int j = i+1; j < n-1; j++)
			{
				for(int k = j+1; k < n; k++)
				{
					t[idx].a = v[i];
					t[idx].b = v[j];
					t[idx].c = v[k];
					idx++;
				}
			}
		}

Now we have a vector of points and a vector of triangles which we will use to process all data and perform all the remaining calculations.

Now, we will turn ourselves to the most tricky parts of this problem, which comprises the writing and correctness of all “helper” functions we will need for the rest of the problem.

  • Making our way out of precision-land without getting lost

Now comes, at last, the bottleneck of this problem: geometry and precision issues.

After we have elegantly got rid of how to represent triangles and points and how to store them (yes, in the process we even got a lovely closed formula for NC3) now comes the time to manipulate them and to introduce yet another new additional complication… Circles!!

Although problem is absolutely trivial, the Geometry concepts it involves, are not so trivial at all and in fact, we are faced with a multitude of challenges now, which is what makes this problem so absolutely awesome:

  1. Given 2 points, find their distance;
  2. Given 3 points, efficiently determine if they form a triangle or not;
  3. Given 2 points, check if they are equal or not;
  4. Given 3 points that are known to form a non-degenerate triangle, calculate the circumcenter and radius of the respective unique circumcircle they belong to;
  5. Given 3 points that are known to form triangle(degenerate or not), calculate the area of that triangle efficiently;
  6. Given a point and the properties of a circle (center coordinates and/or radius), efficiently determine if a point lies inside or on the boundary of that circle;

From the above list, points 4. and 6. are undoubtedly the hardest ones and possibly the ones which were responsible for such low accept rate and for so many WA veredicts.

The main reason for it, imho, is that there are several algorithms to do both tasks. Ones are taught in Linear Algebra classes, while others can be easily looked-up on the Internet, or even deduced by oneself if you like playing with geometric proofs (I don’t! :p).

Some points are trivial, while others will rely heavily on Linear Algebra results and more specifically on the results that follow a Geometrical Interpretation of the Determinant of a Matrix.

For those of you who don’t know what a determinant is, basically, it’s a property that all square matrices have that allows one to infer other useful properties about it, like if the matrix is invertible or if given system of linear equations expressed in matrix form has solution or not, etc. They are very easily computed for matrices 2x2 and 3x3, and you can read their computation on wikipedia :slight_smile:

However, it turns out that historically, Linear Algebra was designed to tackle complex Analytic Geometry problems!! So, it’s quite amazing that a given determinant can also represent the signed area of a given geometric figure, or, even cooler, the signed volume of a given solid.

For instance, let us consider for a moment, the triangle presented below:

![alt text][1]

Now, consider the areas of the trapezoids with top P1P2, with top P2P3 and with top P3P1.

They are respectively:

  • (1/2)(y1 +
    y2)(x1 -
    x2)

  • (1/2)(y2 +
    y3)(x2 -
    x3)

and the one with top P3P1 is the negation of:

  • (1/2)(y3 +
    y1)(x3 -
    x1)

Adding all of the above together, we arrive at the formula:

(1/2)(x1 y2 - x2 y1 + x2 y3 - x3 y2 + x3 y1 - x1 y3)

For anyone who has the least amount of experience with determinants manipulation, it’s now trivial to see that we can re-write the above as:

(1/2)det(A), where A is 3x3 square matrix defined as:

A = |x2  y1  1|
    |x2  y2  1|
    |x3  y3  1|

and det() stands for the determinant of the matrix.

Some mathematical properties of the invariance of determinants over line exchanges and related properties make this value of the “area” to be positive or negative accordingly we choose points in anti-clockwise order or clockwise order respectively, but, we can easily take the absolute value of the area and we are done.

This area-determinant formulation has a multitude of advantages over any other method we could think of, namely, because with it we can immediately cross out from our list the points 5 and 2.

5 is trivially cut out since we just got what we wanted and managed to compute the area of the triangle given the coordinates of its vertexes with just multiplications and subtractions on numbers between -50 and 50 and a single division by 2, which will always be either exact or leave a remainder of 1 (this happens because points have integer coordinates and can be proved) so, we have a safe and accurate way of computing area of a triangle.

As a direct consequence of the above, we see that we have also solved point 2 as if this determinant value evaluates to 0, then it means that the 3 points are on a straight line. We can translate this to code almost directly:

//(signed) Area function is assumed to be 100% correct
inline double area(triangle t)
{
	return (t.a.x*(t.b.y-t.c.y) - t.a.y*(t.b.x-t.c.x) + t.b.x*t.c.y-t.c.x*t.b.y)*0.500000;
}

inline int is_triangle(triangle T)
{
	if(fabs(area(T)) > 0)
		return 1;
	else
		return 0;
}

Points 1 and 3 are also very trivial to implement, as there is a very well known formula for Euclidean Distance of 2 points and to check if two points are equal or not, we simply compare both of their coordinates:

inline double dist(point a, point b)
{
	return sqrt(pow((a.x-b.x),2.0) + pow((a.y-b.y),2.0));
}

inline int points_diff(point a, point b)
{
	if(a.x==b.x and a.y==b.y)
		return 0;
	else
		return 1;
}

With the above code, we also solve points 1 and 3 and we are left with the hardest tasks of solving points 4 and 6 efficiently.

  • Implementing the hardest points with two clever ideas

Both of these last two points were quite hard to implement properly to avoid WA due to the required precision, and it turns out that once again, Mathematics can be amazing and can come to our rescue to help us, lets see how.

Implementing Point 6 using the concept of [Power of a Point][2]

Instead of computing the distance of any given point to the center of circle and checking for the error tolerance using +/- eps or increasing tolerance or some other possibly faulty hack, if we use the idea given on the link, we can simply compare a subtraction result with <= eps and we are assured that the result is well computed.

I believe this to be a more efficient and possibly less error-prone way of dealing with the condition of a point being inside or outside the circle and also especially on its boundary, which was a clear bottleneck during the contest.

[Sidenote: In fact, only after I have written this, I have realized that the cases where a point could lie exactly on the boundary of the circle are pretty much, non-existent, in fact, testing my code below with < eps instead of <= eps, which implies that I am not testing for points on the boundary of the circle, still gets AC, which means that possibly there are no cases like this on the test case set (also, keep in mind that points can’t be coincident, that is, Ash can’t be placed on any of the points where are TR members, which, for a fortunate coincidence (or maybe not), are integer points that lie on the boundary of the circle. Nontheless, if I remove this restriction of points being able to be coincident and do the check for points on the boundary, this gets WA, which can be seen as a sort of guarantee on how this method was actually pretty accurate for this problem).]

So, using this idea we have also solved point 6 efficiently. We can pass it to C++ code simply as:

double s_squared = pow(dist(p[i],center_coord(t)),2.0);
double radius_squared = pow(radius(t),2.0);
if(s_squared - radius_squared <= eps)
// do stuff

where eps = 10-6 is required tolerance.

We are only left with the final point, 4, and we will be done with the code :slight_smile:

  • Implementing point 4:

**

  1. Radius:

The formula for the radius of the circumcircle formed by the triangle with vertexes with known coordinates could be found at wikipedia and it proved to be robust enough for the purpose of this problem.

Implementing it directly in C++ yields:

//radius function is assumed to be 100% correct
inline double radius(triangle t)
{
	return (dist(t.a,t.b)*dist(t.b,t.c)*dist(t.c,t.a))/(4.0000*area(t));
}

**

Now comes the final and real challenge of computing the coordinates of the center of the circumcircle efficiently. Unlike what many people thought, the concept of a circle being unique, was a bit misleading for this problem, as there is always (I repeat, ALWAYS) a single and uniquely determined circle passing trough 3 non-collinear points. One of the challenges of this problem was actually finding a good way of computing the coordinates of the center of such circle as there is absolutely no restriction on those coordinates being integer or not. I will describe my approach for this last step below.

  • The end: Computing the coordinates of the center of the circumcircle using coordinate translation method

There are again, many formulas for this computation… Most direct ones, simply exploit the fact that the center, let’s call it C, must lie within the same distance from the three points that are the vertexes of the triangle.

With those 3 distance equations, we can setup a system of 3 equations and we can solve it to find an expression for the coordinates x and y of the center separately.

As all of the above is already done, some more hours of digging trough Algebra books and the Internet led me to these two formulas:

![alt text][3]

where h and k stand for x and y respectively.

As you shall be able to see from all my WA submissions, this was the formula I used until almost the very end before getting AC, even when I had changed almost all other formulas for its most efficient implementations and used all data types I could remember, I was still getting WA… And, thanks to a comment posted on the problem page, I was inspired to sit back and question everything I had implemented… Especially all the formulas I had taken for granted.

Then, all of a sudden, I remembered that floating point operations in matrices or using matrices properties never scale well and absolute values can be computed with very severe precision errors…

Absolute values was the key here!!

The formula I had been using was computing the coordinates in an absolute manner (that is, it was not using anything as reference and was dealing with squares of distances directly), so I understood I needed to find a new frame of reference, something that would serve the purpose of an anchor which holds a ship in place, this new frame of reference would be my “anchor”.

It was clear, that on my case, my “anchor” could be any of the three points of the triangle (I chose point A), and then, all I needed to do was to compute the “translation factors” and sum them algebraically to the already known coordinates of point A and I would obtain a more accurate answer. This yields true because often these translation factors are, in a sense, intrinsic of the new frame of reference we are using. If we are standing on a point which lies in the same direction as the center of the circumference, then the “translation factors” we need to compute are more numerically stable and can usually be computed with more accuracy then computing the values directly in an absolute manner. This was the missing piece on my implementation for over 40 WA submissions.

So, after googling for numerically stable formulas that would serve the purpose I wanted, I simply switched this WA formula:

//center_coord function is assumed to be 100% correct
point center_coord(triangle t)
{
	point ans;
	long double a,b,c,d,e,f,aa,bb,cc,dd,ee,ff;
	a=t.a.x*t.a.x + t.a.y*t.a.y;
	b=t.a.y;
	c=t.b.x*t.b.x + t.b.y*t.b.y;
	d=t.b.y;
	e=t.c.x*t.c.x + t.c.y*t.c.y;
	f=t.c.y;
	//cout << "Vals from center_coord x: " <<endl << a << " " << b << endl << c << " " << d << endl << e << " " << f << endl;
	ans.x = (a*(d-f) - b*(c-e) + c*f -d*e)/(2.00000*area(t));
	//cout << "X coord= " << ans.x << endl;
	aa = t.a.x;
	bb = a;
	cc = t.b.x;
	dd = c;
	ee = t.c.x;
	ff = e;
	ans.y = (aa*(dd-ff) - bb*(cc-ee) + cc*ff -dd*ee)/(2.00000*area(t));
	return ans;
}

with this new and more numerically stable one, that uses as “anchor”, the vertex A of triangle and got AC without further changes on my Algorithm:

//center_coord function is assumed to be 100% correct
inline point center_coord(triangle t)
{
	point ans;
	double a,b,c,d,e,f,aa,bb,cc,dd,ee,ff;
	
	ans.x = t.a.x - ((((t.b.y-t.a.y)*dist(t.c,t.a)*dist(t.c,t.a)) -  ((t.c.y-t.a.y)*dist(t.b,t.a)*dist(t.b,t.a)))/(2.0*((t.b.x-t.a.x)*(t.c.y-t.a.y)  - (t.c.x-t.a.x)*(t.b.y-t.a.y))));
	ans.y = t.a.y + (((t.b.x-t.a.x)*(dist(t.c,t.a)*dist(t.c,t.a)) - (t.c.x-t.a.x)*(dist(t.b,t.a)*dist(t.b,t.a)))/(2.0* (((t.b.x-t.a.x)*(t.c.y-t.a.y))-((t.c.x-t.a.x)*(t.b.y-t.a.y)))));
	return ans;
}

Now, having the coordinates of the center of the circle correctly computed, all there is left is to calculate answer.

For that, we simply update our answer variable which will store for each triangle, the number of points which lie inside its circumcircle. Updating this variable for all triangles will yield the number of positions where Ash gets killed.

Dividing this value by the total of possible positions for TR (on our case, it’s only size of vector of triangles multiplied by N-3), yields, at last, correct answer. :slight_smile:

  • Final remarks, conclusions and some questions

Some solutions during the contest managed to get down to 0.00sec and “dodged” all precision issues that most people faced, presumably because they managed to avoid doing divisions all together… I am curious to how such method can be actually implemented, as I saw no other possible way besides the one I followed (it’s obvious that there must be very clever and less error-prone and faster ways), but, I am here to learn new things :slight_smile:

Possibly, the general form of the equation of the circle could be used but, after searching over the Internet, I still think we would need to compute the center… So, I hope editorial will clear all my doubts.

Also, I know this can be an easy problem, and possibly many people won’t bother reading this whole thing, but, the great north-american physicist Richard P. Feynman once said:

“No problem is too small or too trivial if you can actually do something about it.”

Richard P. Feynman

This has been my motto for almost everything I have done ever since and this problem, clearly was no exception :smiley:

Also, on a final note about precision issues and use of double vs float, what I learnt is that it’s always better and safer to use double, as it is more accurate than float. In fact, computing the answer for the 2nd test case using float and 15 decimal places of precision, yields:

0.400000005960464

while the same answer computed with the double data type yields the correct answer of:

0.400000000000000

so, as we can see data type can really also play an important role on this problem (and in general in all geometry-related problems I think).

To end this text, all there is left to say is to congratulate both Setter and Tester for having prepared such amazing problem with such well designed corner cases! I have learnt immensely while solving it :slight_smile:

I hope you have enjoyed the explanation of my adventure with this problem and sorry for such long text.

Best regards,

Bruno
[1]: http://discuss.codechef.com/upfiles/cirkill_fig1.jpg
[2]: Power of a point - Wikipedia
[3]: http://discuss.codechef.com/upfiles/cirkill_img2.jpg

26 Likes

Every time there is a problem that requires precision, many users find themselves puzzled. There is a very basic piece of information that all programmers should know.

I am sharing the link. “What Every Computer Scientist Should Know About Floating-Point Arithmetic”

1 Like

Nice one kuruma!

I will add the official editorial shortly! I request you to tag your editorial with the “july13” tag (it is “jul13”) right now, so that everyones able to find it.

I have to refer to this post :smiley:

1 Like

@kuruma : I did this problem without any floating point computations . The approach was similarly based on determinants . You may like to have a look at my solution : CodeChef: Practical coding for everyone .

1 Like

Hi, instead of finding the parameters of the circumcircle of a triangle, you can directly find whether three points form a unique circle or not in the following way:

Equation of a circle (expanded):

x^2 + y^2 + gx + fy + c = 0

After substituting the three points (xi,yi) in the above equation we get three linear equations in g,f and c. If they have a unique solution, then following determinant is not 0 and there exists a unique circle which passes through the three points (xi,yi) :

| x1 y1 1 |

| x2 y2 1 |

| x3 y3 1 |

If indeed it is not 0, then you can find the values of g,f and c by Cramer’s rule. Now you have all the parameters required for the equation of the circle.

Finally to check if a point (x_ash,y_ash) lies inside or on the circle check if p(x) <= 0 (use an epsilon value)

p(x) = x^2 + y^2 + gx + fy + c

P.S. p(x) = 0 when point lies on the circle and p(x) > 0 when point lies outside the circle

My solution implements this idea : CodeChef: Practical coding for everyone

2 Likes

We dont need any floating point calculation:
for each possible 4 points do the following

.
    |x1 y1 1 |
    |x2 y2 1 | = D1
    |x3 y3 1 |
.

if the determinant is 0 the points are collinear, so quit.
else the sign of the determinant tells you if they are arranged CW/anti-CW
Now check the determinant

.
    |x1 y1 x1*x1+y1*y1 1|
    |x2 y2 x2*x2+y2*y2 1|
    |x3 y3 x3*x3+y3*y3 1| = D2
    |x4 y3 x4*x4+y4*y4 1|
.

if D1*D2<=0 then the agent is killed.
see my implementation here

5 Likes

@Kuruma I really appreciate your passion. Its truly inspiring !!

1 Like

Great work kuruma! I’ve always felt you try to give your very best to this CodeChef community. Keep it up! :slight_smile:

I did this problem by almost exactly the same concept as described above by kriateive and it turned out to be a rather simple approach. It uses the final result of the long geometrical proof as described in the official editorial.

1 Like

Learned new approaches from your Work, thank you for your work :slight_smile:
I starting solving this problem using the properties of triangle, concept of circum-circle then I ended up with the equation of circle using 3 points.
Here is my approach

1 Like

@bruno: good work man, keep it up :slight_smile:

1 Like

This solution looks good.

I would like to show my approach.

Equn. of circle passing through 3 points using System of Circles:

When P(x1, y1), Q(x2, y2) and R(x3, y3) are the points:

Equn. of PQ: L: (y-y2 / x-x2) - (y1-y2 / x1-x2) = 0

Equn of Circle with PQ as diameter: C: (x-x1)(x-x2) + (y-y1)(y-y2) = 0

System of circles through the intersection of C and L:

C + kL = 0

=> k = -C/L (if L equals 0 then k = C)

substitute R in C and L to get k.

Now, you get the equation of Circle as C + kL = 0

For any other value D(x, y) (call this guy Ash).

put D in C + kL = 0
if the value is negative (inside the circle) or zero (on the circle) Ash gets killed.

For this approach, we need 4 points P, Q, R => Team Rocket, D => Ash

For each P, Q, R iterate D (find positions for Ash that gets him killed).

Got this idea form here http://www.qc.edu.hk/math/Advanced%20Level/circle%20given%203%20points.htm (method 5).

Let me know your thoughts :slight_smile:

1 Like

Good Work! You should remove the tag “jul13”, which will be removed anyways by the admins some time later. The tag is reserved for editorials.

I think that the tag jul13 is reserved for all questions related to the contest… But, I can also easily remove it :slight_smile:

Hello @shilp_adm I have added the requested tag :smiley: And thanks!!

I can recognize some mathematical expressions there very similar to the ones I have used yes :slight_smile: It seems you also used the 3 points of the triangle and possibly exploited some mathematical and/or geometrical facts that I haven’t used! Nontheless, thanks for showing me a solution without any divisions :smiley: I was curious during contest to see how this was done and now I know

@vineetpaliwal : I saw your solution and I myself employed a similar kind of approach but the center coordinates were being computed differently if the order of points was changed.

Hello @sanchit_h, Yes, after contest I was also told that general equation of circle could be used in the solution, but, I couldnt remember Cramer’s rule at the time :frowning: So I ended up using this more geometric approach :(( I still have much, much to learn… That’s why I am here for :slight_smile: Thanks for you explanation!

1 Like

That’s what I did. Great & simple approach, isn’t it? :slight_smile:

@kuruma Awesome work . Really appreciate your time and dedication to help fellow coders. Keep up the good work buddy !!

2 Likes

@ani94, As I said before, it’s really a pleasure to be able to learn and grow up as a coder with this community and I hope I can stick around for many years to come :slight_smile:

1 Like