[SHOI2015]激光发生器,计算几何 直线相交

https://www.lydsy.com/JudgeOnline/problem.php?id=4595

算是比较明显的线段相交应用题,

涉及到:方向向量,直线相交,点在线段上的判断,直线的方向判断,向量求夹角,向量的旋转

感觉acm赛制下很容易错几次,但是只要写对了总能ac,

而oi赛制想直接拿到满分就很难很难了...

坑点的话,用atan2求向量夹角会出大事,基本过不去

#include<iostream>
#include<streambuf>
#include<string>
#include<cmath>
#include<math.h>
#include<vector>
#include<algorithm>
#include<map>
#include<deque>
#include<stdio.h>
#include<iomanip>
#include<queue>
#define ll long long
#define rep(ii,a,b) for(int ii=a;ii<=b;++ii)
#define per(ii,a,b) for(int ii=b;ii>=a;--ii)
#define forn(i,x,g,e) for(int i=g[x];i;i=e[i].next)
#define IO ios::sync_with_stdio(false);cin.tie(0);cout.tie(0)
#define ull unsigned long long
#define fi first
#define se second
#define mp make_pair
#define db double
#define pb push_back
#define pii pair<db,int>
#define all(x) x.begin(),x.end()
#define show(x) cout<<#x<<"="<<x<<endl
#define show2(x,y) cout<<#x<<"="<<x<<" "<<#y<<"="<<y<<endl
#define show3(x,y,z) cout<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl
#define show4(w,x,y,z) cout<<#w<<"="<<w<<" "<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl
#define show5(v,w,x,y,z) cout<<#v<<"="<<v<<" "<<#w<<"="<<w<<" "<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl
#define showa(x,a,b) cout<<#x<<": ";rep(i,a,b) cout<<x[i]<<' ';cout<<endl
using namespace std;//head
const int maxn=4e5+10,maxm=2e6+10;
const ll INF=0x3f3f3f3f,mod=1e9+7;
int casn,m,k,n;
//@点类与基础@
const db eps = 1e-8;
const db pi = acos(-1);
const db inf = 1e9;
int sign(db x){ if(fabs(x) < eps) return 0;  return x > 0 ? 1 : -1; }
//#define sign(x) x
int cmp(db k1, db k2){ return sign(k1-k2); }
//@k1在k2、k3之间@:
bool inmid(db k1, db k2, db k3){ return sign(k2-k1)*sign(k3-k1) <= 0; }
//@区间相交判定,区间1在区间2前@:
bool intersect(db l1,db r1,db l2,db r2){
  if(l1>r1) swap(l1,r1); if(l2>r2) swap(l2,r2); return cmp(r1,l2)!=-1 && cmp(r2,l1)!=-1;
}
struct point{
	db x, y;
	point(){}
	point(db k1, db k2){ x = k1, y = k2; }
	//@向量加法、点+向量=点:@*/
	point operator + (const point &k1) const { return point(x+k1.x, y+k1.y); }
	//@向量减法、点-点=向量:@*/
	point operator - (const point &k1) const { return point(x-k1.x, y-k1.y); }
	//@向量数乘:@*/
	point operator * (db k1) const { return (point){x*k1, y*k1}; }
  //@向量数除:@*/
  point operator / (db k1) const { return (point){x/k1, y/k1}; }
  //@比较两个点(向量)是否相同:@*/
  bool operator == (const point &k1) const { return cmp(x,k1.x)==0 && cmp(y,k1.y)==0; }
  //@逆时针旋转:@*/
  point turn(db k1){return (point){x*cos(k1)-y*sin(k1), x*sin(k1)+y*cos(k1)};}
  //@逆时针旋转90度:@*/
  point turn90(){return (point){-y, x};}
  //@比较两个点(向量)的大小:@
  //@x越小则点越小,若x相等,则y越小点越小.可以实现按点的坐标排序@*/
  bool operator < (const point k1) const{
    int a = cmp(x, k1.x);
    if(a == -1) return 1; else if(a == 1) return 0; else return cmp(y,k1.y)==-1;
  }
  //@向量模长:@
  db len(){ return sqrt(x*x+y*y); }
  //@向量模长的平方:@
  db len2(){ return x*x+y*y; }
  //@单位向量:@
  point unit(){ return (*this)/(*this).len(); }
  //@向量的极角:@
  db angle() { return atan2(y, x); }
  //@将点放入第一象限:@
  //@当横坐标为负时,或横坐标为0纵坐标为负时,将点按原点做对称 角度是$[-π/2, π/2]$@
  point getdel(){
    if (sign(x)==-1||(sign(x)==0&&sign(y)==-1)) return (*this)*(-1);
    else return (*this);
  }
  //@判断点是否在1 2象限,或者在x的负半轴上 角度是(0, π]@
  bool getp() const {return sign(y)==1 || (sign(y)==0&&sign(x)==-1); }
  void scan(){cin>>x>>y;}
  void print(){cout<<x<<' '<<y<<'\n'; }
};
//@判断k1 在 [k2,k3] 内:@
bool inmid(point k1, point k2, point k3){
  return inmid(k1.x,k2.x,k3.x) && inmid(k1.y,k2.y,k3.y);
}
//@得到两点中点:@
point midpo(point k1, point k2){ return (k1+k2)/2; }
//@两点距离的平方@
db dis2(point k1, point k2){ return (k1.x-k2.x)*(k1.x-k2.x) + (k1.y-k2.y)*(k1.y-k2.y); }
//@两点距离@
db dis(point k1, point k2){ return sqrt(dis2(k1, k2)); }
//@叉乘:@
db cross(point k1, point k2){ return k1.x*k2.y - k1.y*k2.x; }
//@点乘:@
db dot(point k1, point k2){ return k1.x*k2.x + k1.y*k2.y; }
//@向量夹角:@
db rad(point k1, point k2){ 
  return acos(dot(k1,k2)/k1.len()/k2.len()); 
}
//@极角排序,[-π, π]:@
bool compareangle (point k1,point k2){
  return k1.getp()<k2.getp() || (k1.getp()==k2.getp() && sign(cross(k1,k2))>0);
}
bool cmp2(point a,point b){
  if(sign(a.y*b.y)==0) return a.y<b.y;
  if(sign(a.x*b.x)==0) return (a.x<b.x)^(a.y>-eps);
  return a.y/a.x<b.y/b.x;
}
//@k1 k2 k3 逆时针1 顺时针-1 否则0:@
int clockwise(point k1,point k2,point k3){return sign(cross(k2-k1,k3-k1));}
//@直线与线段@
//@直线类@
struct line{ 
//@方向为p[0]->p[1]@
  point p[2];
  line(){}
  line(db x1,db y1,db x2,db y2){p[0]=point(x1,y1),p[1]=point(x2,y2);}
  line(point k1,point k2){p[0]=k1; p[1]=k2;}
  point& operator [] (int k){return p[k];}
//@点在直线左侧的判定:@
//@沿着p0->p1的左侧为1,右侧为0@
  bool include(point k){ 
    return sign(cross(p[0]-k,p[1]-k))>0; 
  } 
//@方向向量:@
  point dir(){return p[1]-p[0];} 
//@向外(左)平移eps@
  line push(){ 
    point delta=(p[1]-p[0]).turn90().unit()*eps;
    return {p[0]-delta, p[1]-delta};
  }
};
//@线段类:@
struct segment{
  point p[2];
  segment(){}
  segment(db x1,db y1,db x2,db y2){p[0]=point(x1,y1),p[1]=point(x2,y2);}
  segment(point a, point b){ p[0]= a, p[1] = b; }
  point dir(){return p[1]-p[0];} 
  point& operator [] (int k){ return p[k]; }
};
//@q 到直线 k1,k2 的投影:@
point proj(point q, point k1, point k2){
  point k=k2-k1; 
  return k1+k*(dot(q-k1,k)/k.len2()); 
} 
//@q 关于直线 k1, k2 的对称点:@
point reflect(point q, point k1, point k2){ 
  return proj(q,k1,k2)*2-q; 
} 
//@点在线段上的判定:@
bool checkons(point q,point k1,point k2){ 
  return inmid(q,k1,k2) && sign(cross(k1-q, k2-k1))==0; 
} 
//@点在直线上的判定:@
bool checkonl(point q,point k1,point k2){ 
  return sign(cross(k1-q, k2-k1))==0; 
} 
//@点在射线k1->k2上的判定:@
bool checkonr(point q, point k1, point k2){ 
  return sign(cross(q-k1, k2-k1)) == 0 && sign(dot(q-k1, k2-k1)) >= 0;
}
//@直线平行判定,可以重合:@
bool parallel(line k1,line k2){ return sign(cross(k1.dir(),k2.dir()))==0; }
//@直线同向判定:@
bool samedir(line k1,line k2){ 
  return parallel(k1,k2)&&sign(dot(k1.dir(),k2.dir()))==1; 
} 
//@直线的比较,极角排序,范围是[-π, π]:@
bool operator < (line k1,line k2){ 
  if (samedir(k1,k2)) return k2.include(k1[0]);
  return compareangle(k1.dir(),k2.dir());
}
//@直线相交判定:@
//@叉积计算面积,两直线不平行必相交(除去重合的情况),平行时,三角形面积相等:@
bool checkll(point k1,point k2,point k3,point k4){ 
  return cmp(cross(k3-k1,k4-k1),cross(k3-k2,k4-k2))!=0;
}
//@直线相交判定:@
bool checkll(line k1,line k2){ 
  return checkll(k1[0],k1[1],k2[0],k2[1]); 
} 
//@直线交点:@
point getll(point k1,point k2,point k3,point k4){ 
  db w1=cross(k1-k3,k4-k3),w2=cross(k4-k3,k2-k3); 
  return (k1*w2+k2*w1)/(w1+w2);
}
//@直线交点:@
point getll(line k1,line k2){
  return getll(k1[0],k1[1],k2[0],k2[1]); 
} 
//@直线与线段相交判定:@
//@线段的两端点在直线的两侧@
bool checkls(point k1, point k2, point k3, point k4){ 
  return sign(cross(k1-k3, k2-k3)) * sign(cross(k1-k4, k2-k4)) <= 0;
}
//@ 线段相交判定:@
bool checkss(point k1,point k2,point k3,point k4){ 
  return intersect(k1.x,k2.x,k3.x,k4.x)&&intersect(k1.y,k2.y,k3.y,k4.y) &&
  sign(cross(k3-k1,k4-k1))*sign(cross(k3-k2,k4-k2))<=0 &&
  sign(cross(k1-k3,k2-k3))*sign(cross(k1-k4,k2-k4))<=0;
}
//@ 线段相交判定:@
bool checkss(segment k1, segment k2){ 
  return checkss(k1[0], k1[1], k2[0], k2[1]); 
}
//@线段规范相交判定:@
//@端点相交不算@
bool strictcheckss(point k1, point k2, point k3, point k4){ 
  return sign(cross(k3-k1,k4-k1))*sign(cross(k3-k2,k4-k2))<0 &&
         sign(cross(k1-k3,k2-k3))*sign(cross(k1-k4,k2-k4))<0;
}
//@ 线段规范相交判定:@
bool strictcheckss(segment k1, segment k2){ 
  return strictcheckss(k1[0], k1[1], k2[0], k2[1]); 
}
//@点到直线的距离:@
db displ(point q, point k1, point k2){ 
  if(k1 == k2) return dis(q, k1);
  return fabs(cross(k2-k1, q-k1)) / (k2-k1).len();
}
//@点到直线的距离:@
db displ(point q, line l){ 
  return displ(q, l[0], l[1]); 
}
//@点到线段的距离:@
db disps(point q,point k1,point k2){
  point k3 = proj(q,k1,k2);
  if (inmid(k3,k1,k2)) return dis(q, k3); 
  else return min(dis(q, k1),dis(q, k2));  
}
//@点到线段的距离:@
db disps(point q, segment k1){ 
  return disps(q, k1[0], k1[0]); 
}
//@线段到线段间的距离:@
db disss(point k1,point k2,point k3,point k4){ 
  if (checkss(k1,k2,k3,k4)) return 0;
  else return min(min(disps(k3,k1,k2),disps(k4,k1,k2)),
                  min(disps(k1,k3,k4),disps(k2,k3,k4)));
}
//@线段到线段间的距离:@
db disss(segment k1, segment k2){ 
  return disss(k1[0], k1[1], k2[0], k2[1]); 
}
//@多边形函数@
//@三角形面积:@
db tarea(point a,point b,point c){
  return fabs((b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x))/2;
}
//@多边形面积:@
//@多边形用 vector<point> 表示 , 逆时针@
db polyarea(vector<point>A){
  db ans = 0;
  sort(all(A),compareangle);
  for(int i=0;i<A.size();i++) ans += cross(A[i],A[(i+1)%A.size()]);
  return fabs(ans/2);
}
//@多边形周长:@
db polyperimeter(vector<point>&A){
  db ans = 0;
  for(int i = 0; i < A.size(); i++) ans += dis(A[i], A[(i+1)%A.size()]);
  return ans;
}
//@多边形重心:@
point polyfocus(vector<point>&A){
  int n = A.size();
  db sumx= 0, sumy = 0, sumarea = 0, area;
  for(int i = 1; i+1 < n; i++){
    area = cross(A[i]-A[0], A[i+1]-A[0])/2.0;
    sumarea += area;
    sumx += (A[0].x+A[i].x+A[i+1].x)*area;
    sumy += (A[0].y+A[i].y+A[i+1].y)*area;
  }
  return point(sumx/sumarea/3.0, sumy/sumarea/3.0);
}
//@点与多边形的位置关系:@
//@ 2 内部 1 边界 0 外部@
int contain(vector<point>&A, point q){
  int pd=0; A.push_back(A[0]);
  for (int i=1;i<A.size();i++){
    point u=A[i-1], v=A[i];
    if (checkons(q,u,v)) return 1;
    if (cmp(u.y,v.y)>0) swap(u,v);
    if (cmp(u.y,q.y)>=0||cmp(v.y,q.y)<0) continue;
    if (sign(cross(u-v,q-v))<0) pd^=1;
  }
  return pd<<1;
  int wn = 0;
  int n = A.size();
  for(int i = 0; i < n; i++){
    if(checkons(q, A[i], A[(i+1)%n])); return -1;//onside
    int k = sign(cross(A[(i+1)%n]-A[i], q-A[i]));
    int d1 = sign(A[i].y-q.y);
    int d2 = sign(A[(i+1)%n].y-q.y);
    if(k > 0 && d1 <= 0 && d2 > 0) wn++;
    if(k < 0 && d2 <= 0 && d1 > 0) wn--;
  }
  if(wn != 0) return 1;//inside
  return 0;//outside
}
//@逆时针凸包判定:@
int checkconvex(vector<point>&A){
  int n=A.size(); A.pb(A[0]); A.pb(A[1]);
  for (int i=0;i<n;i++) if(sign(cross(A[i+1]-A[i],A[i+2]-A[i]))==-1) return 0;
  return 1;
}
//@求凸包:@
//@flag=0 不严格 flag=1 严格@
vector<point> convexhull(vector<point>A, int flag=1){
  int n=A.size(); vector<point>ans(n*2);
  sort(A.begin(), A.end());
  int now=-1;
//@下凸壳@
  for(int i=0;i<n;i++){
    while(now>0 && sign(cross(ans[now]-ans[now-1], A[i]-ans[now-1])) < flag) now--;
    ans[++now]=A[i];
  }
  int pre=now;
//@上凸壳@
  for(int i=n-2;i>=0;i--){
    while(now>pre && sign(cross(ans[now]-ans[now-1], A[i]-ans[now-1])) < flag) now--;
    ans[++now]=A[i];
  }
//@因为A[0]会被算两次,所以舍弃最后一次的A[0]@
  ans.resize(now);
  return ans;
}
//@切割凸包:@
//@保留直线左边的所有点@
vector<point> convexcut(vector<point>A,point k1,point k2){
  int n=A.size(); A.push_back(A[0]); vector<point>ans;
  for(int i=0;i<n;i++){
    int w1=clockwise(k1,k2,A[i]), w2=clockwise(k1,k2,A[i+1]);
    if (w1>=0) ans.push_back(A[i]);
    if (w1*w2<0) ans.push_back(getll(k1,k2,A[i],A[i+1]));
  }
  return ans;
}
//@求半平面交:@
//@半平面是逆时针方向 , 输出按照逆时针@
vector<point> gethalf(vector<line> L){
  int n = L.size();
  sort(L.begin(), L.end());
  int first = 0, last = 0;
//@双端队列指针@
  line *q = new line[n];
//@双端队列@
  point *p = new point[n];
//@p[i]为l[i]和l[i+1]的交点@
  q[last] = L[0];
//@初始化为一个半平面@
  for(int i = 0; i < n; i++){
    while(first < last && !L[i].include(p[last-1])) last--;
    while(first < last && !L[i].include(p[first])) first++;
    q[++last] = L[i];
    if(samedir(q[last], q[last-1])) last--;
    if(first < last) p[last-1] = getll(q[last], q[last-1]);
  }
  while(first < last && !q[first].include(p[last-1])) last--;
  vector<point>ans;
  if(last - first <= 1) return ans;
  p[last] = getll(q[last], q[first]);
  for(int i = first; i <= last; i++) ans.pb(p[i]);
  return ans;
}
int checkpos(line k1,line k2,line k3){return k3.include(getll(k1,k2));}
//@求半平面交:@
//@半平面是逆时针方向 , 输出按照逆时针@
vector<line> gethl(vector<line> L){
  sort(L.begin(),L.end()); deque<line> q;
  for (int i=0;i<(int)L.size();i++){
    if (i&&samedir(L[i],L[i-1])) continue;
    while (q.size()>1&&!checkpos(q[q.size()-2],q[q.size()-1],L[i])) q.pop_back();
    while (q.size()>1&&!checkpos(q[1],q[0],L[i])) q.pop_front();
    q.push_back(L[i]);
  }
  while (q.size()>2&&!checkpos(q[q.size()-2],q[q.size()-1],q[0])) q.pop_back();
  while (q.size()>2&&!checkpos(q[1],q[0],q[q.size()-1])) q.pop_front();
  vector<line>ans; for (int i=0;i<q.size();i++) ans.push_back(q[i]);
  return ans;
}
//@凸包最近点对:@
//@先要按照 x 坐标排序@
bool _cmp(point k1,point k2){return k1.y<k2.y;}
db closestpoint(vector<point>&A,int l,int r){
  if (r-l<=5){
//@当点数小于等于5时,暴力计算:@
    db ans=1e20;
    for (int i=l;i<=r;i++) for (int j=i+1;j<=r;j++) ans=min(ans,dis(A[i],A[j]));
    return ans;
  }
  int mid=l+r>>1; db ans=min(closestpoint(A,l,mid),closestpoint(A,mid+1,r));
  vector<point>B; for (int i=l;i<=r;i++) if (abs(A[i].x-A[mid].x)<=ans) B.push_back(A[i]);
  sort(B.begin(),B.end(),_cmp);
  for (int i=0;i<B.size();i++) for (int j=i+1;j<B.size()&&B[j].y-B[i].y<ans;j++) ans=min(ans,dis(B[i],B[j]));
  return ans;
}
//@凸包的直径(最远点对):@
db convexdiameter(vector<point>&A){
  int n = A.size();
  int now = 1;
  db res = 0;
  for(int i = 0; i < n; i++){
    while(cross(A[i]-A[(i+1)%n], A[i]-A[(now+1)%n]) > cross(A[i]-A[(i+1)%n], A[i]-A[now])){
      now = (now+1)%n;
    }
    res = max(res, dis2(A[now], A[i]));
  }
  return res;
}
//@点集中的最大三角形:@
db maxtriangle(vector<point>&A){
  int m = A.size();
  int a = 1, b = 2;
  db res = 0;
  for(int i = 0; i < m; i++){
		while(cross(A[a]-A[i], A[(b+1)%m]-A[i]) > cross(A[a]-A[i], A[b]-A[i]))
			b = (b + 1) % m;
		res = max(res, cross(A[a]-A[i], A[b]-A[i]) / 2.0);
		while(cross(A[(a+1)%m]-A[i], A[b]-A[i]) > cross(A[a]-A[i], A[b]-A[i]))
			a = (a + 1) % m;
		res = max(res, cross(A[a]-A[i], A[b]-A[i]) / 2.0);
	}
  return res;
}
//@凸包间的最小距离:@
db mindisbetconvex(vector<point>&A, vector<point>&B){
  int n = A.size(), m = B.size();
  if(n < 3 && m < 3){
    if(n == 1){
      if(m == 1) return dis(A[0], B[0]);
      else return disps(A[0], B[0], B[1]);
    }
    else{
      if(m == 1) return disps(B[0], A[0], A[1]);
      else return disss(A[0], A[1], B[0], B[1]);
    }
  }
  int ai = 0, bi = 0;
  for(int i = 0; i < n; i++) if(A[i].y < A[ai].y){ ai = i; }
  for(int i = 0; i < m; i++) if(B[i].y > A[bi].y){ bi = i; }
  db ans = 1e18;
  for(int i  = 0; i < n; i++){
    db ck;
    while(ck = sign(cross(B[(bi+1)%m]-B[bi], A[(ai+1)%n]-A[ai])) < 0) bi = (bi+1)%m;
    if(ck == 0) ans = min(ans, disss(A[(ai+1)%n], A[ai], B[(bi+1)%m], B[bi]));
    else ans = min(ans, disps(B[bi], A[(ai+1)%n], A[ai]));
    ai = (ai+1)%n;
  }
  return ans;
}
//@最小正方形覆盖:@
db minsquarecover(vector<point>&A, db rad){
  db minx = inf, maxx = -inf, miny = inf, maxy = -inf;
  for(int i = 0; i < A.size(); i++){
    point p = A[i].turn(rad);
    minx = min(minx, p.x);
    miny = min(miny, p.y);
    maxx = max(maxx, p.x);
    maxy = max(maxy, p.y);
  }
  return max(maxx-minx, maxy-miny);
}
//@三分--最小正方形覆盖:@
db t_divide(vector<point>&A, db l, db r){
  db m, rm, eps = 1e-8;
  while(r-l >= eps){
    m = l+(r-l)/3;
    rm = r-(r-l)/3;
    if(minsquarecover(A, m) > minsquarecover(A, rm)) l = m;
    else r = rm;
  }
  return minsquarecover(A, (m+rm)/2);
}
vector<line> seg;
vector<int> ans;
db c[maxn];
int main() {IO;
  cout<<fixed<<setprecision(4);
  db x,y,dx,dy;
  cin>>x>>y>>dx>>dy>>n;
  point now(x,y),v0(dx,dy);
  line s(now,now+v0);
  rep(i,1,n){
    db x1,x2,y1,y2,a,b;
    cin>>x1>>y1>>x2>>y2>>a>>b;
    c[i-1]=a/b;
    seg.push_back(line(x1,y1,x2,y2));
  }
  rep(_,1,10){
    int flag=0,id;
    db mn=1e100;
    rep(j,0,n-1){
      if(parallel(s,seg[j])) continue;
      point x=getll(s,seg[j]);
      if(checkons(x,seg[j][0],seg[j][1])&&sign(dot(s.dir(),x-now))>0){
        db d=(x-now).len();
        if(d<mn)mn=d,id=j;
        flag=1;
      }
    }
    if(!flag) break;
    else {
      cout<<id+1<<' ';
      ans.push_back(id+1);
      now=getll(s,seg[id]);
      if(sign(dot(seg[id].dir(),s.dir()))==0){
        s=line(now,now-s.dir());
        continue;
      }
      point v1=seg[id].dir();
      if(sign(dot(v1,s.dir()))<0) v1=v1*-1;
      db al=pi/2-rad(v1,s.dir());
      int bt=sign(cross(v1,s.dir()));
      s=line(now,now+v1.turn(bt*(al*c[id]-pi/2)));
    }
  }
  if(ans.size()==0)cout<<"NONE";
  return 0;
}