poj 1755 （半平面交解不等式）

转自http://blog.csdn.net/non_cease/article/details/7820361

题目链接：http://poj.org/problem?id=1755

题意：铁人三项比赛，给出ｎ个人进行每一项的速度vi, ui, wi;  对每个人判断，通过改变3项比赛的路程，是否能让该人获胜（严格获胜）。

思路：题目实际上是给出了n个式子方程，Ti  = Ai * x + Bi * y + Ci * z ， 0 < i < n

          要判断第i个人能否获胜，即判断不等式组   Tj - Ti > 0,      0 < j < n && j != i    有解

        即 （Aj - Ai）* x + （Bj - Bi） * y + ( Cj - Ci ) * z > 0,   0 < j < n && j != i 有解

         由于 z > 0, 所以 可以两边同时除以 z， 将 x / z, y / z 分别看成 x和 y ， 这样就化三维为二维，可用半平面交判断是否存在解了，

         对每个人构造一次，求一次半平面交即可。

我用的ZZY的 I＆Ｓ算法的模版，做的过程中要将Ａ＊x + B * y + C > 0表示的半平面，转化成由两点组成的向量表示， IQ问题，纠结挺久

首先，所有的半平面保证符号一致（我取的> ), 然后根据 A， B， C的正负构造向量，下面随便画了个小图，帮助理解

我取的每个向量的左边为半平面， 图中 d1， d2就是代码中的，表示的A， B的正负

[cpp]  view plain  copy

1. #include <iostream>  
2. #include <cstdio>  
3. #include <algorithm>  
4. #include <cmath>  
5. using namespace std;  
6.   
7. const double eps = 1e-10;  
8. const double maxl = 1e10;  
9. const int maxn = 105;  
10.   
11. int dq[maxn], top, bot, pn, order[maxn], n, A[maxn], B[maxn], C[maxn];  
12. struct Point {  
13.     double x, y;  
14. } p[maxn];  
15.   
16. struct Line {  
17.     Point a, b;  
18.     double angle;  
19. } tmpL[maxn];  
20.   
21. int dblcmp(double k) {  
22.     if (fabs(k) < eps) return 0;  
23.     return k > 0 ? 1 : -1;  
24. }  
25.   
26. double multi(Point p0, Point p1, Point p2) {  
27.     return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x);  
28. }  
29.   
30. bool cmp(int u, int v) {  
31.     int d = dblcmp(tmpL[u].angle-tmpL[v].angle);  
32.     if (!d) return dblcmp(multi(tmpL[u].a, tmpL[v].a, tmpL[v].b)) > 0;  
33.     return d < 0;  
34. }  
35.   
36. void getIntersect(Line l1, Line l2, Point& p) {  
37.     double dot1,dot2;  
38.     dot1 = multi(l2.a, l1.b, l1.a);  
39.     dot2 = multi(l1.b, l2.b, l1.a);  
40.     p.x = (l2.a.x * dot2 + l2.b.x * dot1) / (dot2 + dot1);  
41.     p.y = (l2.a.y * dot2 + l2.b.y * dot1) / (dot2 + dot1);  
42. }  
43.   
44.   
45. bool judge(Line l0, Line l1, Line l2) {  
46.     Point p;  
47.     getIntersect(l1, l2, p);  
48.     return dblcmp(multi(p, l0.a, l0.b)) <= 0;  /*此处有=，也就是交点p在 l0 上时， dq中最上面的半平面也去掉，因为题目要求严格获胜，也就是最终求出的 
49.                                                                         半平面交为一个点，也认为是无法获胜的 
50.                                                                      */  
51.  }  
52.   
53. void addLine(double x1, double y1, double x2, double y2, Line& l) {  
54.     l.a.x = x1; l.a.y = y1;  
55.     l.b.x = x2; l.b.y = y2;  
56.     l.angle = atan2(y2-y1, x2-x1);  
57. }  
58.   
59. bool halfPlaneIntersection(Line l[], int n) {  
60.     int i, j;  
61.   
62.     for (i = 0; i < n; i++) order[i] = i;  
63.     sort(order, order+n, cmp);  
64.     for (i = 1, j = 0; i < n; i++)  
65.         if (dblcmp(l[order[i]].angle-l[order[j]].angle) > 0)  
66.             order[++j] = order[i];  
67.     n = j + 1;  
68.     dq[0] = order[0];  
69.     dq[1] = order[1];  
70.     bot = 0;  
71.     top = 1;  
72.     for (i = 2; i < n; i++) {  
73.         while (bot < top && judge(l[order[i]], l[dq[top-1]], l[dq[top]])) top--;  
74.         while (bot < top && judge(l[order[i]], l[dq[bot+1]], l[dq[bot]])) bot++;  
75.         dq[++top] = order[i];  
76.     }  
77.     while (bot < top && judge(l[dq[bot]], l[dq[top-1]], l[dq[top]])) top--;  
78.     while (bot < top && judge(l[dq[top]], l[dq[bot+1]], l[dq[bot]])) bot++;  
79.     if (bot + 1 >= top) return false; //当dq中少于等于两条边时，说明半平面无交集  
80.     return true;  
81. }  
82.   
83. void solve() {  
84.     int i, j, k;  
85.     double x1, y1, x2, y2, a, b, c;  
86.      
87.    //给半平面加一个框，这样可以使解x，y都大于0，也可以避免所有半平面交起来后为不为凸多边形，而是一个敞开的区域  
88.   //如果题目输入的不是一个多边形，而是本题这种输入若干不等式组的情况，这样的限定就是必须的，不然有bug，例如，两条线是平行的(但是极角不同)，  
89.   //极角排序后又挨在一起, 那么就可能求它们的交点，就容易出错  
90.    addLine(0, 0, maxl, 0, tmpL[0]);  
91.     addLine(maxl, 0, maxl, maxl, tmpL[1]);  
92.     addLine(maxl, maxl, 0, maxl, tmpL[2]);  
93.     addLine(0, maxl, 0, 0, tmpL[3]);  
94.     for (i = 0; i < n; i++) {  
95.          bool flag = false;  
96.          for (k = 4, j = 0; j < n; j++)  
97.             if (i != j) {  
98.                 a = 1.0 / A[j] - 1.0 / A[i];  
99.                 b = 1.0 / B[j] - 1.0 / B[i];  
100.                 c = 1.0 / C[j] - 1.0 / C[i];  
101.                 int d1 = dblcmp(a);  
102.                 int d2 = dblcmp(b);  
103.                 int d3 = dblcmp(c);  
104.                 /*本人IQ较低，以下这段纠结一个小时。。。 
105.                   下面是根据a*x+b*y+c>0取向量p1p2， 
106.                   其中p1（x1，y1），p2（x2，y2） 
107.                   就是将直线转化为以两点的表示，取向量p1p2左半为半平面 
108.                 */  
109.                 if (!d1) {  
110.                     if (!d2) {  
111.                         if (d3 <= 0) {  
112.                             flag = true; break;  
113.                         }  
114.                         continue;  
115.                     }  
116.                     x1 = 0;  
117.                     x2 = d2;//d2的值为1或-1  
118.                     y1 = y2 = - c / b;  
119.                 }  
120.                 else {  
121.                     if (!d2) {  
122.                         x1 = x2 = - c / a;  
123.                         y1 = 0;  
124.                         y2 = -d1;  
125.                     }  
126.                     else {  
127.                         x1 = 0; y1 = - c / b;  
128.                         x2 = d2;  
129.                         y2 = -(c + a * x2) / b;  
130.                     }  
131.                 }  
132.                 addLine(x1, y1, x2, y2, tmpL[k]);  
133.                 k++;  
134.             }  
135.          if (flag || !halfPlaneIntersection(tmpL, k)) printf ("No\n");  
136.          else printf ("Yes\n");  
137.     }  
138. }  
139. int main()  
140. {  
141.     scanf ("%d", &n);  
142.     for (int i = 0; i < n; i++)  
143.         scanf ("%d%d%d", &A[i], &B[i], &C[i]);  
144.     solve();  
145.   
146.     return 0;  
147. }  

kuangbin写法：

/* ***********************************************
Author        :kuangbin
Created Time  :2013/8/18 19:47:45
File Name     :F:\2013ACM练习\专题学习\计算几何\半平面交\POJ1755_2.cpp
************************************************ */

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
const double eps = 1e-18;
int sgn(double x)
{
    if(fabs(x) < eps)return 0;
    if(x < 0)return -1;
    else return 1;
}
struct Point
{
    double x,y;
    Point(){}
    Point(double _x,double _y)
    {
        x = _x; y = _y;
    }
    Point operator -(const Point &b)const
    {
        return Point(x - b.x, y - b.y);
    }
    double operator ^(const Point &b)const
    {
        return x*b.y - y*b.x;
    }
    double operator *(const Point &b)const
    {
        return x*b.x + y*b.y;
    }
};
//计算多边形面积
double CalcArea(Point p[],int n)
{
    double res = 0;
    for(int i = 0;i < n;i++)
        res += (p[i]^p[(i+1)%n]);
    return fabs(res/2);
}
//通过两点，确定直线方程
void Get_equation(Point p1,Point p2,double &a,double &b,double &c)
{
    a = p2.y - p1.y;
    b = p1.x - p2.x;
    c = p2.x*p1.y - p1.x*p2.y;
}
//求交点
Point Intersection(Point p1,Point p2,double a,double b,double c)
{
    double u = fabs(a*p1.x + b*p1.y + c);
    double v = fabs(a*p2.x + b*p2.y + c);
    Point t;
    t.x = (p1.x*v + p2.x*u)/(u+v);
    t.y = (p1.y*v + p2.y*u)/(u+v);
    return t;
}
Point tp[110];
void Cut(double a,double b,double c,Point p[],int &cnt)
{
    int tmp = 0;
    for(int i = 1;i <= cnt;i++)
    {
        //当前点在左侧，逆时针的点
        if(a*p[i].x + b*p[i].y + c < eps)tp[++tmp] = p[i];
        else
        {
            if(a*p[i-1].x + b*p[i-1].y + c < -eps)
                tp[++tmp] = Intersection(p[i-1],p[i],a,b,c);
            if(a*p[i+1].x + b*p[i+1].y + c < -eps)
                tp[++tmp] = Intersection(p[i],p[i+1],a,b,c);
        }
    }
    for(int i = 1;i <= tmp;i++)
        p[i] = tp[i];
    p[0] = p[tmp];
    p[tmp+1] = p[1];
    cnt = tmp;
}
double V[110],U[110],W[110];
int n;
const double INF = 100000000000.0;
Point p[110];
bool solve(int id)
{
    p[1] = Point(0,0);
    p[2] = Point(INF,0);
    p[3] = Point(INF,INF);
    p[4] = Point(0,INF);
    p[0] = p[4];
    p[5] = p[1];
    int cnt = 4;
    for(int i = 0;i < n;i++)
        if(i != id)
        {
            double a = (V[i] - V[id])/(V[i]*V[id]);
            double b = (U[i] - U[id])/(U[i]*U[id]);
            double c = (W[i] - W[id])/(W[i]*W[id]);
            if(sgn(a) == 0 && sgn(b) == 0)
            {
                if(sgn(c) >= 0)return false;
                else continue;
            }
            Cut(a,b,c,p,cnt);
        }
    if(sgn(CalcArea(p,cnt)) == 0)return false;
    else return true;
}
int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    while(scanf("%d",&n) == 1)
    {
        for(int i = 0;i < n;i++)
            scanf("%lf%lf%lf",&V[i],&U[i],&W[i]);
        for(int i = 0;i < n;i++)
        {
            if(solve(i))printf("Yes\n");
            else printf("No\n");
        }
    }
    return 0;
}