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 .
typedef mpz_class big;
const bool ENABLE_MULTI_PROCESS = true;
using namespace std;
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;
}
}
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) {
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;
}
}
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) ) {
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;
}
double y = lineFormula(Y, s.x, X, dy, dx);
if ( (y >= s.y0 + EPS && y <= s.y1 - EPS) && correctChange(y - Y, dy)) {
processEvent( eventType( s.x,y,X,Y, INTER_VERT ,id ) );
}
}
}
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) {
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() {
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] = '.';
}
}
}
x0 = 2*x0 + 1;
y0 = 2*y0 + 1;
D *= 2;
vector<horzSegment> horz;
vector<vertSegment> vert;
vector<point> magicPoints;
vector<square> evilSquares;
for (int i=0; i<W; i++) {
for (int j=0; j<H; j++) {
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)} );
}
}
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) ) );
}
}
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) {
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) {
vert.push_back( vertSegment( (i+1)*2, j*2, k*2) );
} else {
k++;
}
j = k;
}
}
}
for (int j=0; j<H; j++) {
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) {
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) {
horz.push_back( horzSegment( (j+1)*2, i*2, k*2) );
} else {
k++;
}
i = k;
}
}
}
int DEBUG_DX = -8800, DEBUG_DY = -4100;
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);
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");
}
}
}
for_each(evil, evilSquares) {
svg::rectangle(F, evil->x0, evil->y0, 2,2, "gray");
}
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);
}
svg::circle(F, x0,y0, 0.1, "red");
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) {
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;
if (lastkind != INTER_START) {
doPointIntersection(INTER_START, x0, y0, X, Y, dx, dy, 0);
}
for (int i=0; i<magicPoints.size(); i++) {
doPointIntersection(INTER_MAGIC, magicPoints[i].x, magicPoints[i].y, X, Y, dx, dy, i);
}
for (int i=0; i<horz.size(); i++) {
doHorzIntersection( horz[i], X,Y, dx,dy, i);
}
for (int i=0; i<vert.size(); i++) {
doVertIntersection( vert[i], X,Y, dx,dy, i);
}
if (lastkind != INTER_MAGIC) {
for(int i=0; i<evilSquares.size(); i++) {
doSquareIntersection( evilSquares[i], X,Y, dx,dy, i+1);
}
}
if (found == false) {
if (tdx==DEBUG_DX && tdy==DEBUG_DY) {
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;
}
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;
}
}
usedD += sqrt(event.dist2);
if (usedD >= D + EPS) {
if (tdx==DEBUG_DX && tdy==DEBUG_DY) {
cout << "exceed"<<endl;
}
break;
}
assert(event.dist2 > 0);
lastkind = event.kind;
int id = event.id;
if (lastkind == INTER_START) {
canSee = true;
break;
} else if (lastkind == INTER_MAGIC) {
dx = -dx;
dy = -dy;
X = event.x;
Y = event.y;
} else if (lastkind == INTER_HORZ) {
X = event.x; Y = event.y;
dy = -dy;
} else if (lastkind == INTER_VERT) {
X = event.x; Y = event.y;
dx = -dx;
} else if (lastkind == INTER_SQUARE) {
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;
}
}
if (canSee) {
total ++;
}
}
}
}
svg::end(F);
F.close();
return total;
}
inline void init(){}
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") );
mode = 1;
remove(".finished");
cerr<<"--Multi process mode--"<<endl;
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;
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;
int x = solve();
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;
}