AOJ 2537 : Billiard

問題リンク : Billiard | Aizu Online Judge

問題概要 :
(0,0)が左下で(w,h)が右上となるような長方形があり、
その内部にn個の半径rの円がおいてある
それらは1から順番に番号が割り振られている
1番のボールを(vx,vy)方向に打つ
打たれたボールは10000秒間動きつづける
ボールが壁に当たった場合は鏡の反射のようにして反射する
角にあたった場合は来た方向に戻る
打たれたボールが別のボールと衝突するのであれば、一番最初に衝突するボールの番号を答えよ
どのボールとも衝突しない場合は-1と答えよ

解法 :

ボールの数も少ないのでシュミレーションする

( 以下自分用のメモ、日本語ではない )
まずボールが現在の方向に移動した時、どの壁と衝突するのかを確かめる
各壁について、円の中心と(円の中心+壁と垂直なベクトル * 適当な大きな値 )からなる線分と交差する点を求める
その交点方向に円の中心からrだけ移動した点が(vx,vy)に円が移動した際に一番最初にその壁と衝突する場所なので、そこを記録し壁との衝突ポイントを見つける
そこを元に、衝突までの時間を計算し記録しておく
衝突時間がもっとも短い壁と衝突する
これが２つある時、角に衝突しているのでそういう場合は次の方向は(-vx,-vy)となる
１つなら壁と垂直なベクトルをもとに反射させて次の方向を求める
衝突する時間の最小値と残り時間のうち、小さいものが本当の移動時間となるので注意
円と円が衝突するかは解の公式を使う
毎回方向を正規化したら誤差にやられたので、
正規化は必要なタイミングでのみ行うこと !!!
詳しくはノートRendezvous-1を16ページを見ること

コード:

#include<bits/stdc++.h>
 
#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-10)
#define equals(a,b) (fabs((a)-(b)) < EPS)
#define COUNTER_CLOCKWISE 1
#define CLOCKWISE -1 
#define ONLINE_BACK 2
#define ONLINE_FRONT -2
#define ON_SEGMENT 0
 
using namespace std;
 
class Point{
public:
  double x,y;
 
  Point(double x = 0,double y = 0): x(x),y(y){}
 
  Point operator + (Point p){return Point(x+p.x,y+p.y);}
  Point operator - (Point p){return Point(x-p.x,y-p.y);}
  Point operator * (double a){return Point(a*x,a*y);}
  Point operator / (double a){return Point(x/a,y/a);}
  Point operator * (Point a){ return Point(x * a.x - y * a.y, x * a.y + y * a.x); }
 
  bool operator < (const Point& p) const{ return !equals(x,p.x)?x<p.x:(!equals(y,p.y)&&y<p.y); }
 
  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }
 
};
 
struct Segment{
  Point p1,p2;
  Segment(Point p1 = Point(),Point p2 = Point()):p1(p1),p2(p2){}
  bool operator <  (const Segment& s) const { return ( p2 == s.p2 ) ? p1 < s.p1 : p2 < s.p2; }
  bool operator == (const Segment& s) const { return ( s.p1 == p1 && s.p2 == p2 ) || ( s.p1 == p2 && s.p2 == p1 ); }
};
 
typedef Point Vector;
typedef Segment Line;

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 a.x*a.x+a.y*a.y; }
double abs(Point a){ return sqrt(norm(a)); }
 
int ccw(Point p0,Point p1,Point p2){
  Point a = p1-p0;
  Point b = p2-p0;
  if(cross(a,b) > EPS)return COUNTER_CLOCKWISE;
  if(cross(a,b) < -EPS)return CLOCKWISE;
  if(dot(a,b) < -EPS)return ONLINE_BACK;
  if(norm(a) < norm(b))return ONLINE_FRONT;
  return ON_SEGMENT;
}

bool intersectLS(Line l, Line s) {
  return cross(l.p2-l.p1, s.p1-l.p1)*       // s[0] is left of l
         cross(l.p2-l.p1, s.p2-l.p1) < EPS; // s[1] is right of l
}
Point projection(Line l,Point p) {
  double t = dot(p-l.p1, l.p1-l.p2) / norm(l.p1-l.p2);
  return l.p1 + (l.p1-l.p2)*t;
}
Point reflection(Line l,Point p) {
  return p + (projection(l, p) - p) * 2;
}
 
Point crosspoint(Line l,Line m){
  double A = cross(l.p2-l.p1,m.p2-m.p1);
  double B = cross(l.p2-l.p1,l.p2-m.p1);
  if(abs(A) < EPS && abs(B) < EPS){
    vector<Point> vec;
    vec.push_back(l.p1),vec.push_back(l.p2),vec.push_back(m.p1),vec.push_back(m.p2);
    sort(vec.begin(),vec.end());
    assert(vec[1] == vec[2]); //同じセグメントかもよ
    return vec[1];
    //return m.p1;
  }
  if(abs(A) < EPS)assert(false);
  return m.p1 + (m.p2-m.p1)*(B/A);
}
 
typedef pair<double,double> dd;
const double DINF = 1e20;
bool LT(double a ,double b){ return !equals(a,b) && a < b; }
bool LTE(double a,double b){ return  equals(a,b) || a < b; }
#define pow2(a) ((a)*(a))
 
dd calc(double x1,double y1,double vx1,double vy1,
        double x2,double y2,double vx2,double vy2,double r){
  double VX = (vx1-vx2), X = (x1-x2), VY = (vy1-vy2), Y = (y1-y2);
  double a = pow2(VX) + pow2(VY), b = 2*(X*VX+Y*VY), c = pow2(X) + pow2(Y) - pow2(r);
  dd ret = dd(DINF,DINF);
  double D = b*b - 4 * a * c;
 
  if( LT(D,0.0) ) return ret;
 
  if( equals(a,0.0) ) {
    if( equals(b,0.0) ) return ret;
    if( LT(-c/b,0.0)  ) return ret;
    ret.first = - c / b;
    return ret;
  }
 
  if( equals(D,0.0) ) D = 0;
  ret.first  = ( -b - sqrt( D ) ) / ( 2 * a );
  ret.second = ( -b + sqrt( D ) ) / ( 2 * a );
  if( !equals(ret.first,ret.second) && ret.first > ret.second ) swap(ret.first,ret.second);
  return ret;
}
 
const int MAX_N = 20;
const double TIME = 100000;
int n;
double w,h,r;
Point ps[MAX_N];
Line table[4];
Vector es[4] = {Vector(100,0),Vector(0,100),Vector(-100,0),Vector(0,-100)};
 
typedef pair<double,int> di;
bool COMP(const di &a,const di& b){
  if( !equals(a.first,b.first) ) return LT(a.first,b.first);
  return a.second < b.second;
}
 
bool success = false;
 
struct Data {
  int index;
  double time;
 
  Point rp,cp2;
  bool operator < ( const Data &data ) const {
    if( !equals(time,data.time) ) return LT(time,data.time);
    if( index != data.index ) return index < data.index;
    if( !( rp == data.rp ) ) return rp < data.rp;
    return cp2 < data.cp2;
  } 
};
 
vector<Data> getReflection(double remain,double vx,double vy,Point cur,Vector &nv){
  vector<Data> ret,temp;
  Vector v = Vector(vx,vy);
  rep(i,4){
    Vector h = es[i] * Vector(0,1);
    Line line = Line(cur,cur+h);
    Point cp = crosspoint(table[i],line);
    Vector e = ( cp - cur ) / abs( cp - cur );
    Point ep = cur + e * r;
    Segment seg = Segment(ep,ep+v*TIME);
    if( intersectLS(table[i],seg) ) {
      Point cp2 = crosspoint(table[i],seg);
      double timer = abs(cp2-ep)/abs(v);
      if( equals(timer,0.0) ) continue;
 
      Line wall = Line(cp2,cp2+h);
      Point rp = reflection(wall,ep);
      temp.push_back((Data){i,timer,rp,cp2});
    }
  }
  if( !temp.empty() ) {
    sort(temp.begin(),temp.end());
    double mini = temp[0].time;
    if( LT(remain,mini) ) {
      mini = remain;
      ret.push_back((Data){-1,mini});
    }
    rep(i,temp.size()) if( equals(mini,temp[i].time) ) {
      ret.push_back(temp[i]);
      nv = ( temp[i].rp - temp[i].cp2 );
    }
     
  }
  return ret;
}
 
void compute(double remain,double vx,double vy,Point cur){
  if( equals(remain,0) ) return;
  Vector vtemp = Vector(vx,vy) / abs(Vector(vx,vy));
  vx = vtemp.x, vy = vtemp.y;
  Vector nv = Vector(123,456);
  vector<Data> ref = getReflection(remain,vx,vy,cur,nv);
  double timeLimit = ref[0].time;
  vector<Data> collides;
  REP(i,1,n){
    dd temp = calc(cur.x,cur.y,vx,vy,ps[i].x,ps[i].y,0,0,2*r);
    vector<double> vec;
    vec.push_back(temp.first), vec.push_back(temp.second);
    double mini = DINF;
    rep(j,vec.size()) if( LTE(0.0,vec[j]) && LTE(vec[j],timeLimit) ) mini = min(mini,vec[j],LT);
    if( !equals(mini,DINF)) collides.push_back((Data){i,mini});
  }
  if( !collides.empty() ) {
    sort(collides.begin(),collides.end());
    printf("%d\n",collides[0].index+1);
    success = true;
    return;
  }
  double nvx = vx, nvy = vy;
  if( (int)ref.size() == 2 ) nvx *= -1, nvy *= -1;
  else nvx = nv.x, nvy = nv.y;
  Vector v(vx,vy);
  compute(remain-timeLimit,nvx,nvy,cur+v*timeLimit);
}
 
int main(){
  while( cin >> n, n ){
    double vx,vy;
    cin >> w >> h >> r >> vx >> vy;
    table[0] = Line(Point(0,0),Point(w,0));
    table[1] = Line(Point(w,0),Point(w,h));
    table[2] = Line(Point(w,h),Point(0,h));
    table[3] = Line(Point(0,h),Point(0,0));
    rep(i,n) cin >> ps[i].x >> ps[i].y;
    success = false;
    compute(10000,vx,vy,ps[0]);
    if( !success ) puts("-1");
  }
  return 0;
}