POJ2826 An Easy Problem?!（与心机婊出题人的对决）
这道题看到名字应该就发现有多“难”了，我做完之后只好这样评价“计算几何 == WA很久 + dcmp + EPS + 这TM就一水题！”。
自认为代码没问题的看客们可以先给答案加上一个EPS再提交一遍试试，如果AC了请看本文最后一个部分。
本文的各部分请根据需求选择阅读，不喜勿喷，欢迎提出文中错误。
一、解决本题的基础知识
1.1 叉积的基础知识
学过数学的都知道有个叫向量(vector)的东西，不过这个不是STL里面的那个vector模板，而是几何概念。
为了方便表示，我们把所有的点都看做从原点到该点的向量，这样点之间的减法也有了几何意义。
我们都知道向量点积的定义：(x1,y1) · (x2,y1) == x1 * x2 + y1 * y2，显然满足交换律；代数意义是模之积乘上cos<v1,v2>（v1,v2表示向量），几何意义是一向量对另一向量的投影与其模的的乘积；据此我们可以判断两线段（直线）是否平行（共线）或者垂直，其他的都会出现两种位置关系，这时我们就可以用到叉积了！
叉积的定义：(x1,y1) × (x2,y2) == x1 * y2 - x2 * y1，显然满足反交换律；当然，从公式中是看不出什么的（除非你作死去证明），所以用好结论就够了：令t = v1 × v2，若t > 0则v1在v2的顺时针
180°范围内，若t < 0则v1在v2的逆时针180°范围内，若t == 0则v1与v2同向或反向。（如图）
1.2 两线段判交
在平面上，两直线只有两种关系，平行(parallel)和相交(intersect)，在判断两线段的位置关系时，若这两条线段所在的直线相互平行，它们显然不会相交；
由上图，我们可以得到两线段P1P2和Q1Q2相交的必要条件：(P2 - P1) × (Q1 - P1) * (P2 - P1) × (Q2 - P1) <= 0；
然而这或许还不足以判断它们相交，如下图：
我们发现，在左图中只满足了一个“跨立”条件，即Q1和Q2在直线P1P2的两侧，于是我们应该做两次跨立实验才能保证它们是相交的；而右图则是边界条件的特殊情况，不方便用跨立实验来判断，不过通过观察两张图后，我们发现：当两线段不相交时，以这两条线段为对角线的格点矩形没有交集！于是我们可以借助这一点来增加判交条件以保证算法的正确性，我们称之为快速排斥实验，如下图：
（这张图网上好像到处都是-_-||）
1.3 求交点
首先它们得相交=。=；
可以通过叉积得到的有向面积的比值来判断交点到某一端点的向量，然后加上那个端点就得到交点坐标了。
（懒得画图了= =）
二、AC代码
1 #include<cstdio> 2 #include<iostream> 3 #include<algorithm> 4 #include<ctime> 5 #include<cstring> 6 #include<string> 7 #include<queue> 8 #include<cmath> 9 10 #define file(f) freopen(f".in","r",stdin); freopen(f".out","w",stdout) 11 12 #define min(a,b) ((a)<(b)? (a):(b)) 13 #define max(a,b) ((a)>(b)? (a):(b)) 14 15 using namespace std; 16 17 const double EPS = 1e-8; 18 const double INF = 10000 + 10; 19 20 int dcmp(double a){ 21 if(fabs(a) < EPS) return 0; 22 else return a < 0 ? -1 : 1; 23 } 24 25 struct R{ 26 double dat; 27 void r(){ 28 scanf("%lf",&dat); 29 } 30 void w(){ 31 printf("%.2lf\n",dat); 32 } 33 bool operator == (const R c){ 34 return !dcmp(dat - c.dat);//dat-EPS < c.dat && dat+EPS > c.dat; 35 } 36 bool operator < (const R c){ 37 return dcmp(dat - c.dat) == -1; 38 } 39 bool operator > (const R c){ 40 return dcmp(dat - c.dat) == 1; 41 } 42 bool operator >= (const R c){ 43 return (*this).operator==(c) || (*this).operator>(c); 44 } 45 bool operator <= (const R c){ 46 return (*this).operator==(c) || (*this).operator<(c); 47 } 48 R operator + (const R c){ 49 R a; 50 a.dat = dat + c.dat; 51 return a; 52 } 53 R operator - (const R c){ 54 R a; 55 a.dat = dat - c.dat; 56 return a; 57 } 58 R operator * (const R c) const{ 59 R a; 60 a.dat = dat * c.dat; 61 return a; 62 } 63 R operator / (const R c){ 64 R a; 65 a.dat = dat / c.dat; 66 return a; 67 } 68 }; 69 R zero,two,minu;//zero == (double)0,two == (double)2,minu == (double)-1; 70 71 struct P{ 72 R x,y; 73 void r(){ 74 x.r(),y.r(); 75 } 76 P operator - (const P c){ 77 P a; 78 a.x = x - c.x,a.y = y - c.y; 79 return a; 80 } 81 P operator + (const P c){ 82 P a; 83 a.x = x + c.x,a.y = y + c.y; 84 return a; 85 } 86 P operator * (const R c){ 87 P a; 88 a.x = x * c,a.y = y * c; 89 return a; 90 } 91 P operator / (const R c){ 92 P a; 93 a.x = x / c,a.y = y / c; 94 return a; 95 } 96 bool operator < (const P c){ 97 return mo() < c.mo(); 98 } 99 bool operator > (const P c){ 100 return mo() > c.mo(); 101 } 102 R CP(P c){ 103 return (x*c.y) - (y*c.x); 104 } 105 R DP(P c){ 106 return (x*c.x) + (y*c.y); 107 } 108 R mo() const{ 109 R ret; 110 ret.dat = sqrt((x.dat * x.dat) + (y.dat * y.dat)); 111 return ret; 112 } 113 }; 114 115 struct L{ 116 P f,t;//from,to 117 R k; 118 bool p; 119 void getk(){ 120 P v = f+t; 121 if(t.mo() == zero) k = zero; 122 else if(v.x == f.x){ 123 p = 1; 124 k.dat = INF; 125 } 126 else{ 127 p = 0; 128 k = (v.y-f.y) / (v.x-f.x); 129 } 130 } 131 void r(){ 132 P to; 133 f.r(),to.r(); 134 if(to.y < f.y) swap(f,to); 135 t = to - f; 136 getk(); 137 } 138 bool cr(L c){//判交 139 P A1 = f,A2 = f + t,B1 = c.f,B2 = c.f + c.t; 140 if(!(min(A1.x,A2.x) <= max(B1.x,B2.x) && 141 min(B1.x,B2.x) <= max(A1.x,A2.x) && 142 min(A1.y,A2.y) <= max(B1.y,B2.y) && 143 min(B1.y,B2.y) <= max(A1.y,A2.y) ))//快速排斥实验 144 return false; 145 if((A1-B1).CP(A2-B1) * (A1-B2).CP(A2-B2) > zero || 146 (B1-A1).CP(B2-A1) * (B1-A2).CP(B2-A2) > zero)//跨立实验 147 return false; 148 else 149 return true; 150 } 151 P node(L c){//求交点 152 P A1 = f,A2 = f + t,B1 = c.f,B2 = c.f + c.t; 153 R p = (B1-A1).CP(B2-A1) / ((A1-B2).CP(A2-B2) - (A1-B1).CP(A2-B1)); 154 return (t * p) + f; 155 } 156 }; 157 L X,Y; 158 159 void init(L &l1,L &l2); 160 void work(L l1,L l2); 161 R work2(L l1,L l2); 162 R work3(L l1,L l2); 163 164 int main(){ 165 #ifndef ONLINE_JUDGE 166 file("2826"); 167 #endif 168 L l1,l2; 169 int T; 170 scanf("%d",&T); 171 while(T --){ 172 init(l1,l2); 173 work(l1,l2); 174 } 175 return 0; 176 } 177 178 void init(L &l1,L &l2){ 179 l1.r(); 180 l2.r(); 181 zero.dat = 0.00; 182 two.dat = 2.00; 183 minu.dat = -1.00; 184 Y.t.y.dat = INF,Y.t.x.dat = 0; 185 } 186 187 void work(L l1,L l2){ 188 R ans; 189 if(l1.k == zero || l2.k == zero || (!l1.cr(l2)))//有一条线段平行于x轴或互不相交 190 ans = zero; 191 else if(l1.p || l2.p)//有一条线段平行于y轴(已经排除了平行和不相交) 192 ans = work2(l1,l2); 193 else//还剩三种一般情况：同时左(右)斜或一左一右 194 ans = work3(l1,l2); 195 ans.w(); 196 } 197 198 R work2(L l1,L l2){ 199 if(l2.p) swap(l1,l2);//l1为平行于y轴的线段 200 P N = l1.node(l2),v1,v2; 201 v1 = l2.f + l2.t - N; 202 v2 = l1.f + l1.t - N; 203 if(v1.mo() == zero || v2.mo() == zero) return zero; 204 R cos = v1.DP(v2) / (v1.mo() * v2.mo()); 205 P h = (Y.t / Y.t.mo()) * min((v1 * cos).mo(),v2.mo()); 206 v1 = (v1/v1.mo()) * (h / cos).mo(),v2 = h; 207 return max(v1.CP(v2),v1.CP(v2) * minu) / two; 208 } 209 210 R work3(L l1,L l2){ 211 P N = l1.node(l2); 212 if(abs(l1.k.dat) < abs(l2.k.dat)) swap(l1,l2);//l1相比l2的更靠近y轴 213 P v1,v2; 214 L l3 = Y; 215 l3.f = N; 216 v1 = l1.f + l1.t - N; 217 v2 = l2.f + l2.t - N; 218 if(v1.mo() == zero || v2.mo() == zero) return zero; 219 R cos1,cos2; 220 cos1 = v1.DP(l3.t) / (v1.mo() * l3.t.mo()); 221 cos2 = v2.DP(l3.t) / (v2.mo() * l3.t.mo()); 222 P h1,h2; 223 h1 = (Y.t / Y.t.mo()) * (v1.mo() * cos1); 224 h2 = (Y.t / Y.t.mo()) * (v2.mo() * cos2); 225 if((v1-h1).mo() >= (v2-h2).mo() && l1.k * l2.k > zero) return zero; 226 else{ 227 l1.t = (l1.t / l1.t.mo()) * min(v1.mo(),(v2.mo() * cos2) / cos1); 228 l2.t = (l2.t / l2.t.mo()) * min(v2.mo(),(v1.mo() * cos1) / cos2); 229 l3.t = (l3.t / l3.t.mo()) * l1.t.mo() * cos1; 230 R r1 = l3.t.CP(l1.t) / two,r2 = l3.t.CP(l2.t) / two,ret; 231 ret = r1 - r2; 232 return max(ret,ret * minu);//即abs(ret) 233 } 234 }
三、爬坑指南
1、记住dcmp（三态函数）的作用，做double类型答案的题都要注意精度误差；
2、EPS的选择通常为1e-6,1e-8,1e-12,1e-19，根据题意和需要选取；
3、有时测评机开了严格比较时，答案为-0.00是无法通过的，这时给答案加上一个EPS就变成(+)0.00了（这也是开头那句话的原理）。
四、后记
这只是一道非常简单的计算几何题，但是还是有很多东西需要注意和总结，希望本文对看过的朋友们有所帮助。
