## Friday, April 27, 2012

My objective for this match was to be in the first 333. Let us see how it goes. At this moment there are 5 minutes left for the match. But In fact it was 30 minutes ago that I decided to do other things other than paying attention to the contest. C seemed like unapproachable to me under the duration of the contest.

## Problem A - Password problem

Statement

It seems google are aware that this contest makes programming contests mainstream, as this problem included even a tutorial to explain the common trope that is to calculate the minimum expected number of moves.

My solution involves simply devising a function f(x) that returns the probability at least one of the first x keystrokes is wrong. After we make this function, it is possible to see for example that the expected amount of keystrokes you type after deciding to erase 5 keys is: 5 + (B-A+5) + 1 * f(A - 5)*(B+1). Because we erase 5 characters, then type the remaining (B-A+5) characters, press enter and there is a f(A - 5) probability that it is all wrong.

f(x) can be calculated using dynamic programming. Let p(x-1) be the probability (x-1)-th is right. Then either the key is correct or not: f(x) = (1 - p(x-1)) + p(x-1) * f(x-1). - If it is correct, then we need to multiply the probability that the other (x-1) keys are wrong.

double solve(int A, int B, double* p) {
// dp[A]: probability you mistyped at least one of the first A characters.
double dp[A+1];
dp[0] = 0.0;
for (int i=1; i<=A; i++) {
dp[i] = (1 - p[i-1]) + p[i-1] * dp[i-1];
}

// give up typing, just press enter, it will be wrong all the time
double giveup = 2 + B;
// finish it, type the rest, press enter
double finish = (B - A) + 1 + dp[A]*(1 + B);
double backspace = finish;
// press backspace X times (including 0, which is the same as finish).
for (int back = 1; back <= A; back++) {
backspace = std::min(backspace, back + (B - A + back) + 1 + dp[A-back]*(1 + B) );
}
return std::min(giveup, backspace );

}


## Problem B - Kingdom Rush

Statement

At first I was unsure I could solve B-large. Decided that it was best to first do B-small, then use its results to verify B-large should I think of it. Maybe this was a bad idea, because B-large turned to be far easier than I thought, and the delay caused by solving two problems instead of one was the difference between top 100 and pathetic top 200 :).

The idea in B-large is to always pick the level that gives more stars. Let us say you have a list of amounts of stars you already acquired from each level and the current number of stars. According to your current state, some levels can give you 2 stars, some 1 stars and some can't be solved yet - So, pick one that gives you 2 stars. This should be correct, because all 2-stars levels give you the same benefit and will not make you unable to take other levels (the opposite, in fact).

But do not rush into things, what to do if all levels can give you at most one star? Not all the levels are the same here. One of those levels may be better to take later because with more stars you would take 2 stars from that level in one game. It is better to pick the 1-star level with the largest b_i, this way it is less likely that picking another set of games before this game will allow you to take this level and win 2 stars in one try.

#define make_triple(a,b,c) make_pair(a, make_pair(b,c))
int solve(int N, int* a, int* b) {
int stars = 0, games = 0;
int got[N];
fill(got, got+N, 0);
const int INF = 3*N;
while (stars < 2*N) {
// We use STL pairs for tie breaking purposes.
pair<int, pair<int,int> > best = make_triple(0,0,-1);
for (int i=0; i<N; i++) {
if (got[i] == 1) {
if (b[i] <= stars) {
//can do...
best = std::max(best, make_triple(1, INF, i) );
}
} else if (got[i] == 0) {
if (b[i] <= stars) {
best = std::max(best, make_triple(2, INF, i) );
} else if (a[i] <= stars) {
best = std::max(best, make_triple(1, b[i], i) );
}
}
}
int p = best.second.second;
if (p == -1) {
return -1;
}
games++;
//cout << "game "<<games<<" : "<<p<<endl;
if (got[p] == 1) {
got[p] = 2;
stars++;
} else if (got[p] == 0) {
got[p] = best.first;
stars += got[p];
}

}
return games;
}



B-small is harder to explain, but its main advantage is that you do not need to prove anything and there is no risk that the solution is wrong. The maximum number of levels is 10. And the list has for each game, 3 different values (the number of stars you got from it, 0, 1, 2). There are only at most 310 possibilities for a state. And , let us say that for each possible state, you tried all possible level choices and picked the best one. This dynamic programming would solve the problem.

## Monday, April 23, 2012

### Evil SRM 541 is evil

Another day, another topcoder contest.

Oh well. After round 2A and now this SRM my rating was going to suffer. I did so many stupid things today that it is hard to believe. At the end, once I noticed my rating was going to die, I decided that I may as well increase my volatility so I opted to do a wrong challenge

It was a cool problem set by a new talent. The div1 250 was easy, it just required you to think things over to avoid making the implementation too long. But if you do make the implementation short, you got to be careful not to make a very stupid mistake that was not caught by the examples. The story with div1 550 is... actually the same. These two factors combined into the perfect anti-vexorian storm.

### div1 250: The one with the ants

I was first going to do a complicated bruteforce thing finding equations to get the intersections of each line. In the middle of that business, I noticed that there was an easier (wrong) solution.

The maximum time at which ants can meet is after 4000 moves. Now, we can say that each ant has a pair (dx, dy) which describes its direction (The ant's coordinate change after 1 move, for example , an ant that goes north has (dx=0, dy=1), and an ant that goes west (dx=-1, dy=0). So, the position of the ant after t moves will be: (x[i] + t*dx[i], y[i] + t*dy[i]). Great. So, you can just try all times from 0 to 4000. Then for each pair of (not removed yet) ants, find their positions after the given time, and if their positions are the same, take them out.

This solution is wrong, but only slightly wrong. The issue is that ants will not necessarily meet at a point with integer coordinates. The easiest example would be two ants that initially begin next to each other and each goes to the opposite direction of each other. They will meet exactly one step later in the middle point between their cells.

The solution to this bug is simply to notice that if ants find each other at a coordinate, the coordinate will be integer OR it will be a integer plus 0.5. (Can only meet exactly at the middle of between two cells). Thus you can just use time values like t=0, t=0.5, t=1, ... and do the same thing. Alternatively, you can just multiply every starting coordinate by 2 and also the time limit.

int countAnts(vector <int> x, vector <int> y, string direction) {     int n = x.size();     for (int i=0; i<n; i++) {         // multiply everything by 2         x[i] *= 2;         y[i] *= 2;     }     // find directions dx,dy for each ant:     map<char, int> chdy;     chdy['N'] = 1;     chdy['S'] = -1;     map<char, int> chdx;     chdx['E'] = 1;     chdx['W'] = -1;     int dx[n];     int dy[n];     for (int i=0; i<n; i++) {         dx[i] = chdx[direction[i]];         dy[i] = chdy[direction[i]];     }          // simulate everything     bool active[n];     fill(active,active+n, true);         int t = 0;     while(t <= 4000) { //The time limit was 2000, we multiplied everything by 2         for (int i=0; i<n; i++) if (active[i]) {             for (int j=i+1; j<n; j++) if (active[j]) {                 if ( ( y[i] + dy[i]*t == y[j] + dy[j]*t )                    && ( x[i] + dx[i]*t == x[j] + dx[j]*t )                  ) {                     active[i] = active[j] = false;                }             }         }         t++;     }          return count(active, active+n, true); }

### div1 550: The one with the strings

You know, I thought I was going to solve this during the match, because I sort of had the initial idea the second I read the problem. (This problem reminds me of EndlessStringMachine, which I loved although it is a little harder). But I also knew implementation would be tough to do correctly. I also did not optimize my idea well enough before coding it. Mind you, my idea was all right, but it was overcomplicated and I could not finish debugging it in time (I had an hour to code it and no help). It does not help that I needlessly used dynamic programming and matrix power...

Here is the easier solution: For simplicity, I will rename the variables in the input to left, middle, right, S, F and K. So that using the string means concatenating (left + X + middle + X + right). In fact, I will define a string called hack = left+"*"+middle+"*"+right. This string just puts a * in the places where we have to place X.

First of all, let us solve this problem when K is not very large. There is a simple iterative solution that just requires us to tinker some few things. Let us say that the result after (t-1) steps (times the function gets applied) is known: r[t-1]. Then we can find the result after (t) steps by using the following formula:

r[t] = 2*r[t-1] + A(t)


r[t-1] is the number of substrings in f^(t-1){S} that are equal to F., thus the number of substrings inside the * areas in (left*middle*right) is equal to 2*r[t-1]. We just need to know A(t) - The number of substrings equal to F, that may or may not cover characters from the * area, but definitely cover at least one character from left, middle or right. So when calculating A(t) we only care about the instances of F() that begin or finish in one of the characters that were added at time t.

The main trick is to notice that after a quantity of steps, the values of A(t), A(t+1), A(t+2)... all become constant - the same. Let me introduce you to an easy case: a*b*c, x, axb:

0: x
1: axbxc
2: aaxbxcbaxbxcc
3: aaaxbxcbaxbxccbaaxbxcbaxbxccc
4: aaaaxbxcbaxbxccbaaxbxcbaxbxcccbaaaxbxcbaxbxccbaaxbxcbaxbxcccc
...


The idea is that eventually, next to the characters added at time t, there will always be either a repetition of left (in this example aaaaa...), or right (in this example ...ccccc). There will be a time at which the n left repetitions and the n right repetitions will be larger than the length of F. In fact, this will always happen after at most 50 steps (Because that is the length of F). In conclusion, we only need to find A(t) for t <= 50. After that, we can assume that A(k) = ... = A(52) = A(51) = A(50) and be happy, because the iterative solution res[k] = res[k-1]*2 + A can be found easily because k is at most 10000000.

What is left to do is to find the way to solve res[0] and A(t) for t <= 50.

At time 0, the result string is merely S, so res[0] is simply the number of times F appears as a substring of S.

A(t) is the part of the problem that is hard to implement. It helps that we represent the function as left*middle*right though. First of all, we need to find the lengths of the strings returned by the function for all necessary values of t. len[0] = length(S). len[i] = 2*len[i-1] + (characters in left, right and middle). It is very lucky that t is at most 50 in this calculation, because you will see that in the worst case, len[50] is very large and just a couple of steps away of overflowing signed 64 bits integers, thankfully it does not.

At time t, the length of the string is len[t], of all those len[t] positions, we need to extract the positions that would be good for a substring of size len(F) that overlaps with the added left, middle or right. This procedure is not as difficult as you would think. Simple use hack = left*middle*right and iterate through all characters in hack. But for every *, you got len[t-1] starting positions but only some are valid in a way that overlaps with the added parts. Since a string has to begin or end at a character added in time t, there are at most 300 positions to pick. For each of those positions, we need to extract the substring of len(F) elements that begins at that location at time t and compare it with F.

Extracting the substrings

The substrings will always be of size length(F), this means that you will need length(F) steps to extract them if you decide to instead just extract the characters at a given position. The idea is to make a function get(t, p) that returns the p-th character after time t. To implement this function, once again take advantage of hack = left*middle*right. You need to check if p ends at a new character or at a * character. If it ends at a * then you will need to call get(t-1, new position). Since you recurse at most once per function call, you do not need any fancy implementation tricks and can just implement the recursive approach.

Hopefully the code will clear any misunderstanding:

typedef long long int64; #define long int64  const int MOD = 1000000007;   using namespace std;  const int BASIC = 50; struct AkariDaisukiDiv1 {     string hack;     string S, F;          long len[BASIC];      // Get the p-th character of f^t(p)         char get(int t, long p)     {         if (t == 0) {             if (p >= S.length()) {                 return '-';             }             return S[(int)p];         }         long x = len[t-1];         for (int i=0; i<hack.size(); i++) {             if (hack[i] == '*') {                 if (p < x) {                     return get(t-1, p);                 }                 p -= x;             } else if (p == 0) {                 return ( hack[i]  );              } else {                 p--;             }         }         return '-';     }          // At time t, count the number of substrings equal to F, such that they     // overlap with the added left, middle or right parts.     int A(int t)     {         // we got a string left{X}middle{X}righ{X}         int res = 0;         long p = 0, x = len[t-1];         vector<long> valid;         for (int i=0; i<hack.size(); i++) {             if (hack[i] == '*') {                 for (int o=1; o<=x && o<F.length(); o++) {                     valid.push_back(p+x-o);                 }                 p += x;             } else {                 valid.push_back(p);                 p++;             }         }         for (int i=0; i<valid.size(); i++) {             bool ok = true;             for (int j=0; (j<F.length()) && ok; j++) {                 ok = ( get(t, valid[i] + j) == F[j]);             }             res += ok;         }         return res;     }      int countF(string left, string middle, string right, string S, string F, int k)     {         //F = "dc";         hack = left+"*"+middle+"*"+right;         this->S = S;         this->F = F;          len[0] = S.length();         for (int i=1; i<BASIC; i++) {             int a = hack.length() - 2;             len[i] = 2*len[i-1] + a;         }          // Base case: how many times is F a substring of S?         long a = 0;         long r = 0;         for (int i=0; i<S.size(); i++) {             int x = F.length();             r += (S.substr(i,x) == F );         }          for (int i=1; i<=k; i++) {             // a is the number of substrings equal to F that overlap with             // left, middle or right at time i:             //             // After BASIC steps, a will become constant.             if (i < BASIC) {                     a = A(i);             }             r = ( 2*r + a ) % MOD;         }         return r;     } }; #undef long

## Monday, April 16, 2012

### To solve Google Codejam: Hall of Mirrors I

This is probably be the problem I'll remember the most out of those solved this year. (Don't get me wrong, this is not the hardest problem I solved this year nor it is the hardest to implement, but it was cool and I had some breakthroughs while solving it). I initially had the right idea initially (try angles, simulate movement, count correct) but just the right idea is hardly enough to solve this problem. You need to think things througj and it also helps a lot to understand many things about this problem.

## The "easy" D-small

Something the official contest analysis does not say well is ... something that is supposedly explained in the statement, but you do not really understand it until trying to solve and debugging the third example case. What does it mean to count the number of images of yourself you can see? In reality, this problem asks you to count the number of starting angles for the trajectory of the light such that they will return to the starting point before traveling a distance of at most D units. Once this is understood the official analysis is useful enough.

This is supposedly easier than D-large, and in a way it is, but it is not very easy either. There were at least 15000 coders in the qualification round, but the number of people who solved this problem is around 500. And we all had 24 hours to do it.

I supposedly solved this problem during the match, and I did. But not exactly well. I discovered that there was a way to find the possible distances through a formula. But it seems that I discovered the formula in an improvised way. It turns out that the formula based on the number of bounces I explained before was wrong. Yet the overall solution is correct. Curious thing eh? Well, the thing is that although the formula is wrong, what it does is return the wrong final position (dx,dy) for each combination of vertical and horizontal bounces (bx,by), but it does so in a way that it returns the correct (dx,dy) for a different (bx,by) and that it still finds all needed pairs (dx,dy). As you can see, this mistake did not stop me from solving D-small, but in a way, it made it harder for me to solve D-large, you will see why.

The actual way to solve this problem is, again by first taking the number of horizontal bounces (bx,by). The official contest analysis does a great job explaining both how to think of the formula and why it works. A horizontal bounce can be seen as folding the whole board. You can choose to fold left or right, and thus when bx is negative, you can see it as folding left instead of right. Same about folding up or down. So, you try all pairs (bx,by) and for each of them, you find the appropriate (dx,dy). It is important not to take a pair (dx,dy) that is a multiple or divisor of a pair already picked, so that you do not pick a same direction more than once.

If you are interested in knowing what my mistake was , this is the correct formula:

             int dx, dy;
if (bx & 1) {
dx = W * (bx- sign(bx) ) + 2*( (bx >= 0)?(W-x0):(-x0));
} else {
dx = bx*W;
}
if (by & 1) {
dy = H * (by - sign(by) ) + 2*( (by >= 0)?(H-y0):(-y0));
} else {
dy = by * H;
}


Every pair of horizontal folds, move the point by twice W. If there is an odd amount of folds, then you have to increase (W-x0) or decrease (-x0) , depending if the folds are negative or not. For the y axis, the logic is similar.

## What we can learn from D-small to solve D-large

The solution to D-small is important and the analysis made in the official site is also. Not because we will actually use something like that to solve D-large, but because we need to take some conclusions from it.

First conclusion:: When solving D-small, it is best to multiply all values by 2, so that the starting point is at a lattice point instead of fractioned coordinates. It is also useful in the large problem. Thus from now on, we will assume that D is at most 100, not 50.

Second conclusion:From D-small, you can see that the final direction (dx,dy) will always go towards another point that is also at the center of a square. And the distance between (0,0) and (dx,dy) must be less than D. This actually means that (dx,dy) will always be integer and that the integers are between -100 and 100. More so, since the distance has to be at most 100, then we are only interested in lattice (dx,dy) such that their distance to (0,0) is at most 100. This greatly bounds the number of directions we have to try so the idea to just simulate all light directions is viable. We will from now on call (dx,dy) the direction vector. You know, for every dx units we move in the x-axis we will move dy units in the y axis. Note that dx and dy must be co-prime, so because for example, (2*dx, 2*dy) represents the same direction and should not be counted again.

The direction vector determines the slope of our direction, right. But since abs(dx) and abs(dy) are at most 100. This means two different slopes are never very close to each other. This is very useful, this means that precision is not a big deal. When starting this problem, I was very worried about precision, because if we simulate the movement of light until it returns to a point, we need to know when is light intersecting with a point. An error in calculation can make the point not count as intersected. Generally, using floats/dobules to answer a question that has YES or NO as an answer (And the question, [Does going towards this direction makes you return this very point?] is one of such questions) is very risky because floats are not exact numbers. As I mentioned in the previous post, I was about to try fractions from the GMP library so I had arbitrary precision and avoid this risk, but they were too slow for this problem.

You will see that in most of the operations when simulating the trajectory of the light, the most evil floating point operation we do is division. And division will be performed against dx or dy which both are ints with abs(dx) at most 100. This makes division rather stable. Since slopes cannot be very different, then we can have a generous Epsilon when comparing floats/doubles. In fact, the epsilon I ended up using was 1e-4 !, and it works.

Another thing that shall convince you about precision not being a big problem is that the definition from the statement actually asks us not to consider light bounces when they hit the limits of a segment (mirror). There are some exceptions that we need to consider, like those corners that make the light return in the opposite direction and those corners that "destroy" the light. Though.

## Part 2 - the algorithm

After we extract all those ideas from the small version, we can try to solve the large one. I decided to split this post at this place. So that it does not get too large.

## The general idea

We are located at a point that is initially (X,Y), we move towards a direction vector (dx,dy). This movement will generate a straight line. This line will touch a mirror and then the ray of light bounces or gets destroyed. Whenever a light bounces, we update (dx,dy) and (X,y) accordingly. The process finishes when the ray of light touches the original point, gets destroyed or when its total length exceeds D.

From here on we got to clarify some things.

## Light movement towards a direction vector?

The direction vector (dx,dy) means that for every dx units moved in the horizontal axis, we move dy units in the vertical axis or vice versa. When the starting point of the trajectory is (X,Y). Let us rush to mention that (X,Y) and (dx,dy) together form a straight line. However, note that the direction means that only the positive side of the straight line matters to us, so (dx,dy) not only specify the direction but also the sense of it ... The equation of the straight line is either (y = Y + (x - X)*dy/dx) or (x = X + (y - Y)*dx/dy) (We used the point-slope form).

## Light strikes a mirror

Given a mirror, any mirror, imagine it as a square that is surrounded by segments. The segments are what is important. Basically, the effect of bouncing happens when the line from the equation above intersects a side of a mirror. The good thing about our mirrors is that all of their sides are parallel to either the x or y axis. If we have a horizontal segment that goes from coordinates (x=s.x, s.y) to (x= s.x, s.y1) then the segment's line equation is simply (x = s.x). To find the intersection between the segment and the line from the direction is as easy as replacing x=7 in the first equation. If there is a solution to this equation, we will find a value of y such that the segment's line and the line cross at that point. But of course, we also want the intersection to be inside the segment, this just requires our value of y to be inside [s.y0 , s.y1]. After we confirm the segment intersection, we also need to verify that the point lies towards the correct sense as given by (dx,dy).

There are many special cases all specified by the problem statement. For starters, we should not usually consider the intersection if it is touching the very corner of a mirror. Unless the corner is actually adjacent to another mirror (and thus the sides of the mirrors act as a single segment). There is also the two special cases when the corner either destroys the light ray or makes it bounce completely.

There are many ways to deal with these issues, I will do one that is a little odd but does the job. First of all, consider that mirror segments can come from many different mirrors that stand together but also have an empty space in the appropriate side:

A maximal sequence of # characters in the input that are adjacent and have empty spaces (.) above forms a horizontal segment. Likewise, a maximal sequence that has an empty space (.) bellow forms another kind of horizontal segment. Vertical segments are a similar story. One main idea behind this transformation is that I wanted the logic to simply ignore all intersections between the direction and an end point of a mirror segment, so we need multiple mirrors to form larger segments to avoid light getting lost when it should not. The other idea is to reduce the overall amount of segments, this is a good idea to improve execution speed.

Note how this does not generate any segment for the outer area, that is good because light will never bounce there.

We will ignore intersections with the border of a segment. This introduces a new problem, and it is that sometimes we do need to detect when a ray of light hits a corner. There are two cases: a) The light bounces on a corner (because the mirrors make a concave corner) and b) The light gets destroyed because the light hit a convex corner.

## Magic points and evil squares

I happen to use such silly, graphical language to stuff involved in a problem, because it somehow makes it clearer to my brain. The idea is that a magic point is a corner that allows bouncing. When this corner gets hit by the ray of light, it will be returned to the opposite direction. In order to detect this event, you would need to calculate the intersection between the line equation and a single point (If the point is contained by the line equation and in the correct side). Do not worry about precision, since we at the end want to detect if the ray of light hits the starting point, we will need to do this kind of line-point intersection accurately anyway.

Magical points are simply those points that are adjacent to three mirror squares (#) and an empty square (.) so that they form an L-like shape. They are the corners at which light could bounce in some situations. Evil squares are those squares that could destroy light, they possess corners that are surrounded by 3 empty spaces in a similar L-shape. In the image, magical points are green and evil squares are red:

Detecting when an evil square "destroys" the light is tricky. But it is all about the ray of light not only touching one of its corners (Take a look at the image, one of the evil squares has a corner that is also a magical point). (In the statement, there is also a case in which a ray goes though the corners of two squares adjacent diagonally) . It has to be one of the convex corners. A good way to detect this is if the direction line touches a corner AND also intersects the square itself. The square is made of 4 segments so you can detect an intersection against those segments easily.

There is a catch and it is that the evil square's edges may coincide with the mirror's segments. My suggestion is to consider squares that are slightly smaller (subtract an Epsilon from each dimension). Also, note that hitting the corner is necessary, a square will not destroy light if it does not hit a corner (light will bounce),

In the image, note how both crossing a corner of the square and intersecting with the smaller square are necessary conditions to destroy the light. In practice, the square will be smaller but only slightly smaller, just a distance to separate it from the mirror segments but ensuring that every light that is supposed to be destroyed can still intersect with the evil square.

## The movement algorithm

So let us say we have extracted the starting position, the segments, the magic points and the evil squares from the input. Then we kindly multiplied everything by 2 so that all the relevant coordinates are integer. The next step is to just simulate the light trajectory for each (dx,dy) pair. Thus solve this problem:

• The current position is (X,Y).
• The current direction vector is (dx,dy)
• Light has already traveled usedD distance units.
• Does the current direction take us to a) The starting location. b) A mirror's segment. c) A magic point. d) An evil square or e) The remaining distance (D - usedD) is too short and thus the trajectory ends before doing anything else.
• If the light does hit something, the position (X,Y) gets updated to the point at which the something got hit and when there is reflexion, (dx,dy) have to be updated, how to update it? (used D increases by the distance between original (X,Y) and the new value).

Reflection: Try it in some paper, you shall see that bouncing against a horizontal mirror means multiplying dy by -1. A vertical mirror means multiplying by -1 and a magical point multiplying both dx and dy, by -1.

What object gets hit first?: We now how to know if the light intersects one of the objects or not. But there are usually going to be many objects and of different types. Let us say we have a list of all the intersections between the line that starts at (X,Y) and uses the direction vector (dx,dy), there are going to be many of them, we will call each of these intersections an event. The actual object that will be hit is the object that is closest to (X,Y), so let us pick the event with the minimum distance. When we pick the event with the minimum distance, we move to a new (X,Y) and we update usedD and (dx,dy).

There is a literal corner case that is related to evil squares. The location of the intersection with an evil square is the corner (cx,cy). But if we do things recklessly, this corner could also be a magical point. The solution is to either only check against the relevant (convex) corners of the square or you can also check against all corners, but add priorities to the events, so if two events are at the same distance from (X,Y), you can use the kind of event as a tie breaker and give less priority to evil squares.

In the image, you will see a couple of different thins that can happen moving from (X,Y). The direction line intersects with the segment that is even going through (X,Y), but this should not be the first event - The only way for this to happen is if it actually was the previous event - We should ignore any event within 0 units of distance to avoid any infinite loop of picking the same location all the time - . The next closest event is the corner of an evil square, but it does not count - Since the line is not going through the square's edges, there is no intersection to it. Next closest event is a magical point, this is the closest valid event, so the result is to pick this event. There are two remaining events which both are valid but since the magical point happens before them, we have to ignore them - There is an intersection with an evil square, and this one is valid because it crosses the square and touches the intersection - {Even though this intersection is valid, it should not be. Luckily, when this happens there should always be another event at a lower distance when this issue of hitting the square at the wrong corner.

## Part 3 - the implementation

After we extract all those ideas from the small version, we can try to solve the large one. I decided to split this post at this place. So that it does not get too large.

## Debugging tools

The algorithm does not sound very easy to implement. There are possible pitfalls everywhere. Starting from even the part in which we translate the input to segments and other objects. How about this idea? Graphical debugging. In geometry problems like this, it really helps to be able to look at exactly what is going on rather than read variables. The idea is to have some sort of visual feedback about a specific case.

We could use some good graphics environment from which we generate and then output images of what the light trajectory is doing by request. There are many ideas for this. I hesitated to try a library like SDL or openGL - I just did not feel like setting a build environment or remembering the calls. I also did not need a lot of features, I just needed two things: The ability to draw lines, rectangles and circles and ability to save the results in files. Eventually I settled for SVG. SVG is based on XML, so people that already know XML can find out quickly how to generate SVG files. I just used the following link: http://www.svgbasics.com/shapes.html. I reiterate that I had no idea how to make SVG files until Saturday when I was getting ready to code this solution. It was quite easy. Who knows? Perhaps SVG is the best friend of a coder trying a geometry problem for a contest... The images on the top of this post were all generated by the SVG code in my solution, and they show some special cases dealt masterfully.

SVG was very useful, but to debug it you need to first know what case (and starting angle) to try it on. This is where the solution to D-small can come into play. If you solve D-small, you have a simple solution that stays in integer world and also a whole file of test cases and another file with their correct answers. You can reuse the small input on the code you are developing for D-large and make sure the results coincide. You can also just try crafting some test cases.

When there is a disagreement between your current D-large solution and your correct solution for D-small, we can print all the pairs (dx,dy) that were successful in both solutions and then compare the outputs, any pair (dx,dy) that works in one solution but does not work in the other is a good idea for generating images for.

An anecdote about debugging: Do not assume things up. At my last hours trying this problem I assumed that my solution for D-small was correct in its logic. But it was not, as I mentioned in the first part, this solution was correct overall, but the logic had issues. The problem is that the mistake in my solution was making it output the wrong (dx,dy) values as valid. Thus when attempting to debug the third case in D-large, I was trying a value (dx = 5, dy = -6) that was actually not valid, but I tried it under the wrong assumption that it was. In my visual debugging it looked like the light was destined to reach the center, but for some reason it was not, this was confusing and lead me to believe that numerical instability was making the precision hard to maintain. But further debugging (thanks to images) helped me confirm that the case was that D-small was wrong. (As we shown in part I, precision is not a problem).

## Epsilons

We use floating point arithmetic everywhere in this problem. But that means that many times, when doing the line - point intersection sub-problem that is so necessary to detect the moment the ray of light returns to its starting position, the floating point numbers might vary and thus you might end not detecting the intersection even though it was meant to.

The lack of precision in floating points is the reason we will use Epsilon values when comparing coordinates and intervals in this problem. So, instead of wanting an exact match in the line - point intersection result. We will admit points that are EPSILON-units away of the line as contained by it.

We need the tolerance in other places. For example, we had defined that mirror segments should not allow intersections to their corners. Let us say we have a horizontal segment with y coordinate s.y, and x coordinates between s.x0 and s.x1 and we used the line equation to find an x value that intersects with the light direction. x must not only be strictly between s.x0 and s.x1 - You can assume that the borders (part to ignore) of the segment begin at s.x0 and s.x1 an finish at s.x0+EPS and s.x0-EPS.

## Optimizations

The logic described in part 2 works. And it seems that if implemented correctly it takes around 2 minutes even in this Core 2 duo. (The time limit is 8). But if you want to take the solution to the next level, you have to consider that in the current version of the solution needs to try intersections between all objects in the test case. One solution is to use a sweep line algorithm and with the help of a data structure (maybe a binary tree) - So that at every step, we only try against the objects that are within (D - usedD) distance units. But it seems this was not necessary.

## Code

It is rather long this time. Note that I used a trick to use both cores of my core 2 duo. This speed boost is not as .

#ifdef VX_PRECOMPILED #include "precomp.h" typedef mpz_class big; #else #include<string> #include<iostream> #include<sstream> #include<.h> #include<cstdio> #include<map> #include<algorithm> #include<bitset> #include<cmath> #include<queue> #include<functional> #include<set> #include<sys/stat.h> #include<numeric> #include<cstdio> #include<cstdlib> #include<cstring> // http://gmplib.org/ (uncomment if solution does not use big nums) #include "gmpxx.h" #define big mpz_class #endif  #define for_each(q,s) for(typeof(s.begin()) q=s.begin(); q!=s.end(); q++) const bool ENABLE_MULTI_PROCESS = true; using namespace std;  #include <fstream>  //============================================================================== // quick SVG stuff // I coded most of this on Saturday during the qualification round while // chatting in the topcoder arena. But it was way disorganized and mixed with // the rest of the code, so I reorganized it a little. namespace svg {     void start(ofstream &F, int W, int H) {         F  << "<?xml version=\"1.0\" standalone=\"no\"?>" << endl;         F  <<"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << endl;         F << "<svg viewBox = '0 0 "<<W<<" "<<H<<"' version = '1.1'>" << endl;     }     void rectangle(ofstream &F, int x, int y, int width, int height, const char* color) {         F<<"<rect x = '"<<x<<"' y = '"<<y<<"' width = '"<<width<<"' height = '"<<height<<"' fill = '"<<color<<"' stroke = '"<<color<<"' stroke-width = '0'/>"<<endl;     }     void circle(ofstream & F, int cx, int cy, double r, const char* color)     {         F << "<circle cx='"<<cx<<"' cy='"<<cy<<"' r='"<<r<<"' fill='"<<color<<"' stroke='"<<color<<"' stroke-width='0' />" << endl;      }     void line(ofstream &F, double x1, double y1, double x2, double y2, const char* color, double width)     {          F << "<line x1 = '"<<x1<<"' y1 = '"<<y1<<"' x2 = '"<<x2<<"' y2 = '"<<y2<<"' stroke = '"<<color<<"' stroke-width = '"<<width<<"' />"<<endl;     }     void end(ofstream &F) {         F << "</svg>"<<endl;     } }   //========================================================= // program: // int TestCase = 0; const double EPS  = 1e-4;  int W, H, D; char grid[31][31];  int gcd(int a, int b) {     if (a < 0) {         return gcd(-a, b);     } else if (b < 0) {         return gcd(a, -b);     } else if (a < b) {         return gcd(b,a);     }     while (b != 0) {         int c = b;         b = a % b;         a = c;     }     return a; }  struct horzSegment {     int y, x0, x1;          horzSegment(int y, int x0, int x1) {         this->y = y; this->x0 = x0; this->x1 = x1;     } }; struct vertSegment {     int x, y0, y1;          vertSegment(int x, int y0, int y1) {         this->x = x; this->y0 = y0; this->y1 = y1;     } }; struct point {     int x, y;     point(int x0, int y0) {         x = x0; y = y0;     } }; struct square {     int x0,y0, x1, y1;      };  enum interkind {     INTER_START = 0, INTER_MAGIC = 1, INTER_HORZ = 2, INTER_VERT = 3, INTER_SQUARE = 4 };  bool compDouble(double a, double b) {     return (a+EPS >= b && a-EPS <= b); }   double doubleDist2(double x0, double y0, double x1, double y1) {     double dx = x0 - x1;     double dy = y0 - y1;     return dx*dx + dy*dy; }  struct eventType {     double dist2, x, y;     interkind kind;     int id;          eventType( double x1, double y1, double x0, double y0, interkind kind, int id) {         this->x = x1; this->y = y1;         this->dist2 = doubleDist2(x0,y0,x1,y1);         this->kind = kind; this->id = id;     }     eventType() {}; }; eventType event; bool found = false;  void processEvent( const eventType & other ) {     if (other.dist2 < EPS) {         //ignore those events that are too close. (Fixes infinite loops).         return ;     }     if (! found ) {         found = true;         event = other;     } else if (compDouble(other.dist2, event.dist2) ) {         if ( (int)(event.kind) > (int)(other.kind) ) {             event = other;         }     } else if ( event.dist2 > other.dist2 ) {         event = other;     } }   // Equations: y = Y + (x - X)*dy/dx; // Equations: x = X + (y - Y)*dx/dy; double lineFormula(double X, int y, double Y, int dx, int dy) {     double m = dx /(double)dy;     return X + (y - Y)*m; }   bool correctChange(double dif, int dx) {     if (dx == 0) {         return dif >= -EPS && dif <= EPS;     }     if (dx > 0) {         return (dif > -EPS);     }     return (dif < EPS); }  void doPointIntersection(interkind kind, int px, int py, double X, double Y, int dx, int dy, int id) {     if (! correctChange(px - X, dx)) {         return;     }     if (! correctChange(py - Y, dy)) {         return;     }      if (dx == 0) {                  double x = lineFormula(X, py, Y, dx, dy);         if ( compDouble(x,px)) {             processEvent( eventType( px,py,X,Y, kind ,id ) );         }     } else {          double y = lineFormula(Y, px, X, dy, dx);         if ( compDouble(y,py)) {             processEvent( eventType( px,py,X,Y, kind ,id ) );         }     }  }  void doHorzIntersection(horzSegment &s, double X, double Y,  int dx, int dy, int id) {     if (dy != 0) {         if (! correctChange(s.y - Y, dy)) {             return ;         }         double x = lineFormula(X, s.y, Y, dx, dy);         if ( (x >= s.x0+EPS && x <= s.x1 - EPS) && correctChange(x - X, dx) ) {             // intersection!             processEvent( eventType( x,s.y,X,Y, INTER_HORZ ,id ) );         }     } }   void doVertIntersection(vertSegment &s, double X, double Y,  int dx, int dy, int id) {     if (dx != 0) {         if (! correctChange(s.x - X, dx)) {             return;         }          // Y + (s.x - X) * (dy / dx)         double y = lineFormula(Y, s.x, X, dy, dx);         if ( (y >= s.y0 + EPS && y <= s.y1 - EPS) && correctChange(y - Y, dy)) {             // intersection!             processEvent( eventType( s.x,y,X,Y, INTER_VERT ,id ) );         }     } } // This intersection is for the segments in an evil square, so it is different... bool lineHorzIntersection(double y, double x0, double x1, double X, double Y,  int dx, int dy) {     if (dy != 0) {         if (! correctChange(y - Y, dy)) {             return  false;         }         double x = lineFormula(X, y, Y, dx, dy);         return ( (x >= x0 && x <= x1) && correctChange(x - X, dx) );     }     return false; }  const double EPS2 = 1e-4; void doSquareIntersection(square &s, double X, double Y,  int dx, int dy, int id) {     int t = 0;     t += lineHorzIntersection( s.y0+EPS2, s.x0+EPS2, s.x1-EPS2, X,Y,dx,dy);     t += lineHorzIntersection( s.y1-EPS2, s.x0+EPS2, s.x1-EPS2, X,Y,dx,dy);     t += lineHorzIntersection( s.x0+EPS2, s.y0+EPS2, s.y1-EPS2, Y,X,dy,dx);     t += lineHorzIntersection( s.x1-EPS2, s.y0+EPS2, s.y1-EPS2, Y,X,dy,dx);          if (t > 0) {         // intersection!         for (int i=0; i<2; i++) {             for (int j=0; j<2; j++) {                 int px = s.x0 + 2*i;                 int py = s.y0 + 2*j;                 doPointIntersection(INTER_SQUARE, px, py, X, Y, dx, dy, id);             }         }     } }    int solve() {     // Find the starting point:     int x0,y0;     for (int i=0; i<W; i++) {         for (int j=0; j<H; j++) {             if (grid[i][j] == 'X') {                 x0 = i; y0 = j;                 grid[i][j] = '.'; //we will later need all chars to be # or .             }         }     }     //Convert     x0 = 2*x0 + 1;     y0 = 2*y0 + 1;     D *= 2;          vector<horzSegment> horz;     vector<vertSegment> vert;     vector<point> magicPoints;     vector<square> evilSquares; //made of evil too.     for (int i=0; i<W; i++) {         for (int j=0; j<H; j++) {             // Detect "evil" squares, you know, a mirror that has three empty             // squares in at least one of its corners.             if (grid[i][j] == '#') {                 bool evil = false;                 for (int x=i-1; x<=i+1; x+=2) {                     for (int y=j-1; y<=j+1; y+=2) {                         bool valid = (x>=0 && x<W && y>=0 && y<H);                          evil |= ( grid[x][y] == '.' && grid[x][j] == '.' && grid[i][y] == '.');                         }                 }                 if (evil) {                     evilSquares.push_back( (square){2*i,2*j,2*(i+1),2*(j+1)} );                 }             }             // Detect magical points, they are a corner adjacent to three squares.             int c = 0;             for (int x=i; x<W && x<i+2; x++) {                 for (int y=j; y<H && y<j+2; y++) {                     if (grid[x][y] == '#') {                         c++;                     }                 }             }             if (c==3) {                 magicPoints.push_back( point(2*(i+1), 2*(j+1) ) );              }         }         //Vertical segments         if (i>=1) {             int j=0;             while (j < H) {                 int k = j;                 while ( (k < H) && (grid[i][k] == '#') && (grid[i-1][k] == '.') ) {                     k++;                 }                 if (k != j) {                     //a segment                     vert.push_back( vertSegment(i*2, j*2, k*2) );                 } else {                     k++;                 }                 j = k;             }         }         if (i<W-1) {             int j=0;             while (j < H) {                 int k = j;                 while ( (k < H) && (grid[i][k] == '#') && (grid[i+1][k] == '.') ) {                     k++;                 }                 if (k != j) {                     //a segment                     vert.push_back( vertSegment( (i+1)*2, j*2, k*2) );                 } else {                     k++;                 }                 j = k;             }         }     }     for (int j=0; j<H; j++) {         // Horizontal segments         if (j>=1) {             int i=0;             while (i < W) {                 int k = i;                 while ( (k < W) && (grid[k][j] == '#') && (grid[k][j-1] == '.') ) {                     k++;                 }                 if (k != i) {                     //a segment                     horz.push_back( horzSegment(j*2, i*2, k*2) );                 } else {                     k++;                 }                 i = k;             }         }         if (j<H-1) {             int i=0;             while (i < W) {                 int k = i;                 while ( (k < W) && (grid[k][j] == '#') && (grid[k][j+1] == '.') ) {                     k++;                 }                 if (k != i) {                     //a segment                     horz.push_back( horzSegment( (j+1)*2, i*2, k*2) );                 } else {                     k++;                 }                 i = k;             }         }             }     // Debug direction. This is the direction whose traject will be printed to a image.     int DEBUG_DX = -8800, DEBUG_DY = -4100;     //----------------------------------------     // Debug file stuff     string fname = "";     {         char buf[100];         sprintf(buf, "thing%03d.svg", TestCase);         fname = buf;     }     ofstream F;     F.open(fname.c_str() );     svg::start(F, W*2, H*2);      //-------------------------------------     // Draw the background objects. :     //     // Normal squares for reference (they do not interact with the solution but make     // the image easier to interpret     for (int i=0; i<W; i++) {         for (int j=0; j<H; j++) {             if (grid[i][j] == '#') {                 svg::rectangle(F, i*2, j*2, 2,2, "#B0B0B0");             }         }     }     // "evil" squares     for_each(evil, evilSquares) {         svg::rectangle(F, evil->x0, evil->y0, 2,2, "gray");            }     //vertical and horizontal segments:     for_each(seg, vert) {         svg::line(F, seg->x, seg->y0, seg->x, seg->y1, "black", 0.1);     }     for_each(seg, horz) {         svg::line(F, seg->x0, seg->y, seg->x1, seg->y, "black", 0.1);     }     //the start:     svg::circle(F, x0,y0, 0.1, "red");     // "magic" points:     for_each(p, magicPoints) {         svg::circle(F, p->x,p->y, 0.1, "darkgreen");     }          int total = 0;     for (int tdx=-100; tdx<=100; tdx++) {         for (int tdy=-100; tdy<=100; tdy++) {             if (!tdx && !tdy) continue;             if ( tdx*tdx + tdy*tdy >= D*D ) continue;                          int g = gcd(tdx,tdy);             if (g==1) {                 //try this direction!                 int dx = tdx, dy = tdy;                 double X = x0, Y = y0;                                  interkind lastkind = INTER_START;                                  bool canSee = false;                 double usedD = 0;                 while (usedD < D+EPS) {                     found = false;                     //intersection with the start?                     if (lastkind != INTER_START) {                          doPointIntersection(INTER_START, x0, y0, X, Y, dx, dy, 0);                     }                     // intersection with magic points?                     for (int i=0; i<magicPoints.size(); i++) {                         doPointIntersection(INTER_MAGIC, magicPoints[i].x, magicPoints[i].y, X, Y, dx, dy, i);                     }                     // intersection with horizontal segments?                     for (int i=0; i<horz.size(); i++) {                         doHorzIntersection( horz[i], X,Y, dx,dy, i);                     }                     // intersection with vertical segments?                     for (int i=0; i<vert.size(); i++) {                         doVertIntersection( vert[i], X,Y, dx,dy, i);                     }                     if (lastkind != INTER_MAGIC) {                         // intersection with evil squares?                         for(int i=0; i<evilSquares.size(); i++) {                             doSquareIntersection( evilSquares[i], X,Y, dx,dy, i+1);                         }                     }                      if (found == false) {                         //No intersection found? This should never happen -                           // the hall is always a closed area. So this it is likely                          // a bug if we got here. Current solution won't get here                         // under the specified constraints.                         if (tdx==DEBUG_DX && tdy==DEBUG_DY) {                             // Draw a line that gets lost.                             double a = atan2(dy,dx);                             svg::line(F, X, Y, X+(D-usedD)*cos(a), Y+(D-usedD)*sin(a), "#000030", 0.05);                             cout << "not found"<<endl;                         }                         break;                      }                                          // Something to do in case we are doing the debug direction                     // Just draw the trajectory...                     if (tdx==DEBUG_DX && tdy==DEBUG_DY) {                         double dist = min<double>(sqrt(event.dist2), D - usedD );                         double ang = atan2(dy,dx);                         double X2 = X + dist*cos(ang);                         double Y2 = Y + dist*sin(ang);                         svg::line(F, X, Y, X2, Y2, "#000030", 0.05);                         switch ( event.kind )  {                         case INTER_START:                             cout << "posible "<<usedD<<" + "<<dist<<endl;                             break;                         case INTER_VERT:                             cout << "vertical "<<usedD<<" + "<<dist<<endl;                             break;                         case INTER_HORZ:                             cout << "horizontal "<<usedD<<" + "<<dist<<endl;                             break;                         case INTER_MAGIC:                             cout << "magic "<<usedD<<" + "<<dist<<endl;                             break;                          }                       }                     // Update the distance. Seems sqrt is stable enough.                     usedD += sqrt(event.dist2);                      if (usedD >= D + EPS) { //Epsilons are a necesary evil here                         if (tdx==DEBUG_DX && tdy==DEBUG_DY) { //more debug stuff                             cout << "exceed"<<endl;                         }                         break;                     }                     // This assert was useful to find odd situations that lead                     // to infinite loops, like continuously intersecting with the                     // same object.                     assert(event.dist2 > 0);                                          // Update the last kind of intersection.                     lastkind = event.kind;                                          int id = event.id;                     if (lastkind == INTER_START) {                         // We intersect the starting point again? Good stuff.                         canSee = true;                         break;                     } else if (lastkind == INTER_MAGIC) {                         //ooh a magical point! ... reverse the direction                         dx = -dx;                         dy = -dy;                                                  X = event.x;                         Y = event.y;                     } else if (lastkind == INTER_HORZ) {                         // bounce horizontally                         X = event.x; Y = event.y;                         dy = -dy;                     } else if (lastkind == INTER_VERT) {                         // bounce vertically                         X = event.x; Y = event.y;                         dx = -dx;                     } else if (lastkind == INTER_SQUARE) {                         // "destroyed" , this happens when the ray of light gets                         // destroyed by a corner. Some debug stuff, it is a good                         // idea to color the infringing square red, because it turns                         // out these squares caused a lot of bugs in the past...                         if (tdx == DEBUG_DX && tdy == DEBUG_DY) {                             cout << "destroyed by square"<<endl;                             {                                 svg::rectangle(F, evilSquares[id-1].x0, evilSquares[id-1].y0, 2,2, "red" );                             }                         }                         break; // The part that breaks.                     }                  }                 if (canSee) {                     // During development, it was important to print                     // the directions we found to be valid.                     //cout<<"("<<tdx<<", "<<tdy<<") "<<usedD<<endl;                     total ++;                 }             }         }     }     svg::end(F);     F.close();      return total; }  inline void init(){} //========================================================= // I/O: // int main(int argc, const char* argv[]) {     int mode = 0;     if(argc >= 2) sscanf(argv[1],"%d",&mode);     if ( ( mode == 0 ) && ENABLE_MULTI_PROCESS )     {         string inputfile = ".input";         ( system("cat > .input") );         /* I use a dual core CPU, so for long solutions I might use this          multi-process thing, splitting the input in halves effectively          halving execution time of slow solutions. But due to overhead it          increases the time of fast solutions, so it is optional... */         mode = 1;         remove(".finished");         cerr<<"--Multi process mode--"<<endl;         //string inputfile = argv[2];         string command = argv[0];         command += " 2 < "+inputfile+" > .tem &";         ( system(command.c_str()) );         ( freopen(inputfile.c_str(),"r",stdin) );     }          init();     int TestCases;     cin>>TestCases;      for (int _i=1;_i<=TestCases;_i++) {         TestCase = _i;         /* read input goes here */         cin >> W >> H >> D;         for (int i=0; i<W; i++) {             for (int j=0; j<H; j++) {                 cin >> grid[i][j];             }         }                  if( (mode==0) || ( (mode!=2)== (_i*2<=TestCases) ) ) {                         cerr<<"["<<_i<<" / "<<TestCases<<"]"<<endl;             /* program call goes here */             int x = solve();             /* output result goes here */             cout<<"Case #"<<_i<<": "<<x<<endl;                     }         else if(mode!=2) break;                  assert(cin);     }     if(mode==2) {         ( system("echo done > .finished") );     } else if(mode==1) {         struct stat stFileInfo;         while( stat(".finished",&stFileInfo)!=0);         ( system("cat .tem") );     }     return 0; }