11130:Billiard bounces

問題文
http://uva.onlinejudge.org/external/111/11130.html

横a,縦bインチのビリヤード台があって、その中央から玉を反時計周りにA度の角度に速度vで打ち出す。この玉はs秒間動くが、摩擦でスピードが減ってしまう。ビリヤード台の縦の壁と横の壁に玉が当たった回数を出力せよ。

何も思いつかなかったのでシミュレーションした。
本当は割るだけらしい。

#include<cmath>
#include<algorithm>
#include<iostream>
#include<vector>
#include<climits>
#include<cfloat>
#define curr(P, i) P[(i) % P.size()]
#define next(P, i) P[(i+1) % P.size()]
#define prev(P, i) P[(i+P.size()-1) % P.size()]
 
 
using namespace std;
 
double EPS = 1e-10;
const double PI = acos(-1);
 
double add(double a, double b){
  if(abs(a+b) < EPS * (abs(a)+abs(b)))return 0;
  return a+b;
}
 
struct point{
  double x, y;
  point(){}
  point(double x,double y) : x(x) , y(y){}
 
  point operator + (point p){
    return point(add(x,p.x), add(y,p.y));
  }
 
  point operator - (point p){
    return point(add(x,-p.x), add(y,-p.y));
  }
 
  point operator * (double d){
    return point(x*d,y*d);
  }
 
  point operator / (double d){
    return point(x/d,y/d);
  }
   
  bool operator == ( const point &p ) const {
    return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS;
  }
};
 
#define Line pair<point,point>
 
double dot(point a, point b){return (a.x * b.x + a.y * b.y);}
double cross(point a, point b){return (a.x * b.y - a.y * b.x);}
double norm(point a){return sqrt(a.x*a.x+a.y*a.y);}
double dist(point a,point b){return sqrt(pow(a.x-b.x,2)+pow(a.y-b.y,2));}
 
point symmetry(point p1, point p2, point Q){
  double xa,ya,xb,yb,t1,t2;
  point R;
 
  xa=Q.x-p1.x,ya=Q.y-p1.y;
  xb=p2.x-p1.x,yb=p2.y-p1.y;
   
  t1=xa*xb+ya*yb,t2=xb*xb+yb*yb;
   
  R.x=2*(p1.x+xb*t1/t2)-Q.x;
  R.y=2*(p1.y+yb*t1/t2)-Q.y;
  return R;
}
 
point rotate(point t,point p,double r){
  double ta=cos(r)*(t.x-p.x)-sin(r)*(t.y-p.y)+p.x;
  double tb=sin(r)*(t.x-p.x)+cos(r)*(t.y-p.y)+p.y;
  return point(ta,tb);
}
 
int is_point_on_line(point a, point b, point c) {
  return cross(b-a, c-a)==0.0 &&
         (dot(b-a, c-a) > -EPS) &&
         (dot(a-b, c-b) > -EPS);
}
 
int is_intersected_ls(point a1, point a2, point b1, point b2) {
 
  if(cross(a1-a2,b1-b2)==0){
    return is_point_on_line(a1,a2,b1) || is_point_on_line(a1,a2,b2) 
        || is_point_on_line(b1,b2,a1) || is_point_on_line(b1,b2,a2);
  }
  else {
    return ( cross(a2-a1, b1-a1) * cross(a2-a1, b2-a1) < EPS ) &&
           ( cross(b2-b1, a1-b1) * cross(b2-b1, a2-b1) < EPS );
  }
}
 
point intersection_ls(point a1, point a2, point b1, point b2) {
  point b = b2-b1;
  double d1 = abs(cross(b, a1-b1));
  double d2 = abs(cross(b, a2-b1));
  double t = d1 / (d1 + d2);
 
  return a1 + (a2-a1) * t;
}
 
 
int main(void){
 
  point p(0,0),q;
  double v,a,s;
 
  while(true){
    cin >> q.x >> q.y >> v >> a >> s;
    if(q.x==0 && q.y==0 && v==0 && a==0 && s==0)break;
 
    vector<point>pol;
    pol.push_back(point(0,0));
    pol.push_back(point(q.x,0));
    pol.push_back(q);
    pol.push_back(point(0,q.y));
 
    p=q/2;
   
    point e(v*s/2+p.x,p.y);
    e=rotate(e,p,a*(PI/180));
 
    int vrtcl=0,hrzntl=0;
 
    while(dist(p,e)>0){
      bool fg=true;
      for(int i=0;i<pol.size();i++){
    if(is_intersected_ls(curr(pol,i),next(pol,i),p,e) &&
       !is_point_on_line(curr(pol,i),next(pol,i),p)){
      fg=false;
      point ne=symmetry(curr(pol,i),next(pol,i),e);
      p=intersection_ls(curr(pol,i),next(pol,i),p,e);
      e=ne;
 
      for(int j=0;j<pol.size();j++){
        if(p==pol[j]){
          vrtcl++,hrzntl++;
          goto end;
        }
      }
       
      if(i%2==0)hrzntl++;
      else vrtcl++;
 
    end:;
      break;
    }
      }
      if(fg)break;
    }
   
    cout << vrtcl << " " << hrzntl << endl;
  }
  return 0;
}