## 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; }

#### 2 comments :

vexorian said...

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 .

Hmnn, incomplete thought is incomplete.

"This speed boost is not as important as other ideas, like cropping the number of directions to check. This approach needs around 2 minutes in my computer using one core and 1 minute using two."

Artyom Razinov said...

There is a way to don't use either floating point operations or even division operation to solve this task! This task is completely discrete. The large input can be processed in 0.2 sec.