MAXCIR - Editorial






Sweep line algorithm, Computing square roots


You are given three points A, B and C and N pairs (x[i], y[i]). You can shift coordinates of A by at most K of these pairs. Each pair can be used at most once. What maximal value can get |AB| + |BC| + |CA| in this situation?


Clearly |BC| is constant. So we should consider only |AB| + |AC| and add |BC| in the end. The key idea is to maximize some linear function f(X, Y) = U * X + V * Y instead of maximizing the value |AB| + |AC| that involves square roots. Here (X, Y) are the coordinates of the new position of A and U, V are some fixed parameters. Namely, the following statement holds.

Lemma. There exists some real numbers U and V such that the set of pairs, applied to the point A, that maximizes the function f(X, Y) = U * X + V * Y, also maximizes the sum of the distances |AB| + |AC|.

Proof. Consider the set of at most K pairs that maximizes the sum of distances |AB| + |AC|. Let this maximum be D and the corresponding point A, on which this maximum is achieved, be A’. So for each other set of at most K pairs this sum of distances is not greater than D. Note, that the set of points A for which |AB| + |AC| <= D is a set of points enclosed by an ellipse with foci at B and C that passes through the point A’. Consider the tangent line L to this ellipse at the point A’. Let its equation be U * X + V * Y + W. Since ellipse encloses the convex set and all the points A, corresponding to other sets of at most K pairs, lies inside this set, all these points A will also lie at the same half-plane with respect to the line L. Hence, we can assume, changing signs of U, V and W if necessary, that U * X + V * Y + W <= 0 for all other points A. And since for the point A’ we have equality here it exactly means that the point A’ maximizes the linear function f(X, Y) = U * X + V * Y. We can think on this from another position. Namely, the set of pairs that maximizes considered linear function also maximizes the sum of distances. So we are done.

As we see, the proof do not provide any hints on how to find this U and V. So in order to solve the problem we ideally should try all possible values of U and V, maximize the corresponding linear function and update the actual answer by the value |AB| + |AC| for the optimal point A of the current linear function. Of course, it’s unreal plan, so we should decrease the set of different pairs (U, V), that we will check, to some relatively small set of candidates.

For this consider the problem of maximization of linear function f(X, Y) = U * X + V * Y in detail. At first it is clear that we can remove all completely zero pairs (x[i], y[i]). So now we assume that every given pair is non-zero. Denote by F[i] the value U * x[i] + V * y[i]. If we apply to A the set of pairs with indexes i1, …, iL and obtain the point (X, Y), then by linearity of f we get f(X, Y) = f(Ax, Ay) + F[i<sub>1</sub>] + … + F[i<sub>L</sub>]. It is clear now, how to find the set of pairs that maximizes f(X, Y). If the number of positive values among F[1], …, F[N] is no more than K than choose all pairs with positive value, otherwise choose the K largest by F[i] pairs.

This observation shows that the only essential things that we should know for the given U and V are, which of F[i] are positive, and what their relative order. Clearly F[i] = 0 if and only if (U, V) is parallel to (y[i], -x[i]) and F[i] = F[j] if and only if (U, V) is parallel to (y[i] - y[j], x[j] - x[i]). If for some i and j we have equal pairs then F[i] = F[j] for all U and V and we shouldn’t consider direction produced by this pair. Lines parallel to these directions and passing through the origin divide the plane into several sectors. Clearly the signs of F[i] and their relative order are constant for all (U, V) in each particular sector. So we can choose just one particular pair (U, V) from each sector, find the optimal set of at most K pairs in O(N * log N) time (since we in general need to sort them) and update the actual answer, which is the sum of squares |AB| + |AC|. In this way we get the O(N3 * log N) solution since there are O(N2) sectors.

We can improve this solution to O(N2 * log N) solution. For this we need to maintain the signs and relative order of F[i] during iterations through the sectors without actual values of F[i], as well as, the optimal point.

Let’s sort all critical directions by polar angle (each critical line produce two opposite critical directions). This can be made by comparing at first half-planes of two directions and then by the sign of their cross-product. I suggest you to not use attractive atan2() function in C++ that calculates the polar angle, since in general situation it may leads to precision issues. In this problem it will be enough. But if coordinates of directions can be up to 109 then there can be two different directions that can’t be distinguished by atan2(). For example

**atan2(109, 109 + 1) = 0.78539816389744830936...**; **atan2(109 - 1, 109) = 0.78539816389744830986...**. As you see, they differ only at 19-th digit after the decimal point which is out of `double` precision.

At the beginning we directly calculate relative order and signs distribution of F[i] for some direction (U, V) in the last sector (it is in fact the previous to the first one). For example one can take (U, V) = (2 * maxX + 1, -1), where maxX = 106 is the maximal possible absolute value of x[i] and y[i].

Next we consider critical directions one-by-one. If for the current direction F[i] = 0 then we change the corresponding sign. If for this direction F[i] = F[j] for some i and j we swap their positions. During both these transformations it is clear how the sum of optimal set of pairs is changed. But very special consideration is required when several events occur for the same direction. Currently, please, see tester’s solution as a reference. It currently has no explanatory comments. It will be improved soon, as well as the required explanation will be added here.

Now let’s discuss how to deal with real numbers in this problem with 13 digits after the decimal point. We are very disappointed that nobody figured it out the beautiful way we invent for this. All accepted solutions have cumbersome calculation of square root using some usual algorithm. Common guys! Are you really thought that we want you to add a ton of ugly code just to print crappy 13 digits after the decimal point? If we would really decide to see how you could calculate the square root we would ask to print, say, 100000 digits.

In this problem the only thing that is needed to deal with square roots with required precision is the following simple formula:

**sqrt(A) = M + (A - M * M) / (M + sqrt(A))**,

where M = floor(sqrt(A)). Now the idea is to deal with real numbers here as with pairs (I, F), where I is an integer part and F is a fractional part, so 0 <= F < 1. It is clear how to compare and add such numbers (see tester’s solution as a reference). And the mentioned formula shows that for sqrt(A) we have I = M and F = (A - M * M) / (M + sqrt(A)). Here F should be calculated directly using built-in function sqrt() since we only need F with 13 significant digits.

You should use this not only for printing the answer but also for comparing current answer with some candidate answer. We have tests where very close different choices exist. The simplest such test is the following:
1 1
0 1000000000
1000000000 -1000000000
-1000000000 -1000000000
1 0
We can either use the spell card or not. The answers in these cases are:
So their relative error is 5.5279e-20 and can’t be distinguished by double.

Well, we even have several such tests involving large tests with many different sectors. So we had no doubt that solution with comparison in double will fail. But contestants as usually cheated on us. Check it out this solution and also this one.

Also in the first solution instead of smart iterating over critical directions simple heuristic to try many random directions was used. In fact, we have tests where this heuristic used with correct maximization of f(X, Y) will fail even after 1000000 of random tries. These very special tests were designed by the tester and the probability to guess the direction is lower then 1/11,000,000 in these tests. But look, in this solution for the given direction all prefix sets are tried instead of choosing the optimal set. In fact, after creating those tough test data we develop this hacking solution as well and saw that it easily passes, but we had no time at all to think how to fail it. It was one of the reasons of such delay with this problem. In the second solution some much smarter scheme to crack our test data was used.

After some checks we now realize that by adding the following test
1 1
1 1000000000
1000000000 -1000000000
-1000000000 -1000000000
-1 0
along with the mentioned small test (which presents in test data) we would fail comparison in double.

Finally, let us say a few words about tricky overflows in this problem. The first one can occur when the coordinates of the shifted point A are subtracted from the coordinates of the points B and C, if this calculation is performed in 32-bit signed integer type. This will happen on the following test
500 500
1000000000 1000000000
-1000000000 -1000000000
-1000000000 -1000000000
1000000 1000000 (500 times)
The difference in this case can be 2500000000 = 2.5e9 > 231.

The same test also causes overflow during calculation of the sum of the squares (Bx - X)2 + (By - Y)2, where (X, Y) are the coordinates of the shifted point A, if this calculation performed in 64-bit signed integer type. This sum of squares can equal to 1.25e19 > 263. But using unsigned 64-bit integer type will resolve the issue. We can proudly said that even ACRush was caught by this overflow for some time :slight_smile:


Setter’s solution will be updated soon.
Tester’s solution can be found here.


CPOINT - Codechef April 2011 Long Contest - Central Point.
The tester has invented this hack with square root exactly during participation at that contest (see his solution here).


Really nice problem. Do you have an idea on how many points can there be in the convex hull of all possible uses of the cards?

I suspected it couldn’t be much because how using a card affects all the points, and I couldn’t think of a way of having many more points than n^2, so I tried to look at the problem as a knapsack problem, for each card, try to use it and keep the convex hull since the solution lies in the convex hull, obviously it got TLE.

I’d like to point out that I really dislike problems for which 64 bit unsigned integers are required and 64 bit signed integers overflow. Some languages such as Java do not have proper support for such integers; thus your point about not wanting people to use hacky methods isn’t really valid as that is the only option in Java.


actually i am not able to understand the second part of changing signs and other stuffs.
Can any one explain it explicitly.

If we ignore the k limitation, when a new card is added at most 2 points increased in the new convex hull. So there are at most n*2 points in the final convex hull. But in this problem this method doesn’t work, because the k cards limitation leads to a o(n^3logn) complexity.

We thought on “without limit” version at some point. We have an O(N * log N) solution for such version since we don’t need to consider directions of the form (y[j]-y[i],x[i]-x[j]). But decided that “limit version” is smarter, though it allows some cheating solutions to pass.

I don’t really understand this formula: sqrt(A) = M + (A - M * M) / (M + sqrt(A)) , we have sqrt(A) on both sides, so how does it help?

You can calculate convex hull in O(n), so it’s acutally O(n^3).

If we construct the convex hull by adding cards, we should record how many cards have used by each point. And the points inside the convex hull should not be dropped, because these may have used fewer cards thus may get a better result after adding some new cards.
So…you mean a convex hull corresponding to each number of used cards? O(n^3) also got TLE?

I actually didn’t find a way to bound the number of points in the convex hull, I think this might be the reason I get TLE. But if we have K convex hulls, and each of them has O(n) points, then we have an O(k*n^2) algorithm, for each card, for each convex hull, add all the points of the convex hull with 1 less card, you can merge all the points of two convex hulls in O(n) and calculate the new convex hull in O(n).

“But look, in this solution for the given direction all prefix sets are tried instead of choosing the optimal set.”

But why does it work “better” than finding the optimal set?

@mmaxio It is explained further: M is in integer part and (A - M * M) / (M + sqrt(A)) is a fractional part. Since we need fractional part only with 13 digits we can calculate it directly by this formula using usual sqrt. The main thing is that you should not add M and fractional part since you immediately lost many significant digits in fractional part. See tester solution for details.

@damians Well, it is better at least because more sets of pairs are checked. Most probably some of these sets are optimal for another close sectors.

Definitely backing tripleM on this one. Java is slow enough without being forced to use arbitrary precision arithmetic while the c languages use the standard stuff.

Also, hacky? I dunno what everyone else did but I used binary search to calculate the square root. Binary search should be the most natural thing on Earth to a computer scientist and definitely not hacky.

1 Like

Regarding long long overflow. The constraints was set quite naturally and nobody thought about long long overflow. It was found accidentally and we decided to stay on this version anyway.

I did not use word “hacky” for computation of square roots, that all of you guys did. Your luck that we can’t find tests where a lot of very close possibilities exist. So you all would get TLE with your “natural” way of computing square root on such tests since you should perform this too many times.

@venuswitharms I am curious what this lines in your solution mean?

if(++rtCnt > 4*n+10)

Is this some proved bound or you just hoped that test data are weak?

Also without high percision version I am sure we have a lot of really hacky solutions that tries random angles and so on. This would ruin the problem completely. But precision issues of the current version scares away most of the potential cheaters.

The reason I mention it is that this is the second problem in a few months in which the editorials specifically mention you should be working with 64 bit unsigned integers, as if it were deliberate: CLOSEST - Editorial - editorial - CodeChef Discuss

In that case you could work around it by subtracting 2^63 from each value when comparing, but if this indeed is only happening due to a lack of testing, it’s something that should be carefully checked in future contests. Will always pop up when using distances involving -10^9 to 10^9.

Well, in fact, we could change constraints but I like this overflow very much :P.

In CLOSEST it was me who noticed it, but again since I like it, we did not change the constraints :stuck_out_tongue:

So it is not due to a lack of testing but because of too neat testing :slight_smile:

Well, I only can suggest you to ask Oracle (Sun Microsystems) to add unsigned long long to Java :wink:

But OK, next time I will try to consider you point.

I will not test next two long contests so you can safely participate in them without risk to meet this overflow again :wink:

I really liked the problem btw :slight_smile: