半平面交可以解决一些线性规划问题和纯几何问题，比较常用于求多边形核等。
半平面交有两种写法，一种是比较直观的$O(n^2)$写法，一种是$O(n \log n)$写法。
$O(n^2)$的写法的思路是：设当前已经得到了一个半平面交的有序点集p（若半平面交封闭则必然是凸包），对于当前正在处理的直线，$O(n)$查找它能否通过切割将凸包的面积缩小，同时算出它与凸包边界的交点，从而更新凸包。
$O(n \log n)$的思想于此类似，但是使用双端队列优化，主要思路是：
step1. 将所有半平面按极角排序，对于极角相同的，选择性的保留一个。 O(nlogn)
step2. 使用一个双端队列(deque)，加入最开始2个半平面。
step3. 每次考虑一个新的半平面：
a.while deque顶端的两个半平面的交点在当前半平面外:删除deque顶端的半平面
b.while deque底部的两个半平面的交点在当前半平面外:删除deque底部的半平面
c.将新半平面加入deque顶端
step4.删除两端多余的半平面，具体方法是：
a.while deque顶端的两个半平面的交点在底部半平面外:删除deque顶端的半平面
b.while deque底部的两个半平面的交点在顶端半平面外:删除deque底部的半平面
step5:计算出deque顶端和底部的交点即可。
首先第一步将所有直线按照斜率排序，同时每条直线都确定一个方向使得半平面一定在它的“左侧”，这样对于斜率相同的直线一定只有最“靠左“的直线会起作用。
然后先加入前两个半平面，接着依次加入半平面。最后删除多余的不符合条件的半平面。最后要记得计算顶端和底部的交点使凸包封闭。
对于两点式直线求交的问题，可以直接列式解方程。
$y_1=k_1x+b_1$，$k_1=\frac{l_1.y_2-l_1.y_1}{l_1.x_2-l_1.x_1}$（$l_1.y_2$表示第一条直线的第二个点的纵坐标）， $y_2$同理
将上式代入$k_1x+b_1=k_2x+b_2$化简即可。
用两点式时对于平行与x轴或y轴的直线也一样处理，解直线不等式时，对于直线$ax+by+c \leqslant 0$,若$a=0$则加入$(0,-c/b)$和$(1,-c/b)$，若$b=0$则加入$(-c/a,0)$和$(-c/a,1)$，若均不为0则加入$(0,-c/b)$和$(1,-a/b-c/b)$。
$O(n^2)$的写法：
POJ3335求多边形核：
1 #include<cstdio> 2 #include<cmath> 3 #include<algorithm> 4 #define rep(i,l,r) for (i=l; i<=r; i++) 5 using namespace std; 6 7 const int N=210,inf=0x3f3f3f3f; 8 const double eps=1e-8; 9 inline double sgn(double x){ return (fabs(x)<eps) ? 0 : ( x>0 ? 1 : -1); } 10 11 int i,T,ntot,n,tot; 12 struct P{ double x,y; }p[N],a[N],s[N]; 13 struct S{ P s,t; }; 14 15 double cross(P a,P b,P c) { return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x); } 16 bool outside(S seg,P p){ return cross(seg.s,seg.t,p)>eps; } 17 bool inside(S seg,P p){ return cross(seg.s,seg.t,p)<-eps; } 18 19 void work(P p1,P p2,P p3,P p4){ 20 double t1=cross(p1,p2,p3),t2=cross(p1,p2,p4); 21 if (sgn(t1)*sgn(t2)<0){ 22 double x=(fabs(t2)*p3.x+fabs(t1)*p4.x)/(fabs(t1)+fabs(t2)); 23 double y=(fabs(t2)*p3.y+fabs(t1)*p4.y)/(fabs(t1)+fabs(t2)); 24 a[++tot]=(P){x,y}; 25 } 26 } 27 28 void cut(S seg){ 29 int i; tot=0; 30 rep(i,1,ntot){ 31 if (!outside(seg,p[i])) a[++tot]=p[i]; 32 else work(seg.s,seg.t,p[i-1],p[i]),work(seg.s,seg.t,p[i],p[i+1]); 33 } 34 ntot=tot; swap(a,p); p[0]=p[tot]; p[tot+1]=p[1]; 35 } 36 37 void solve(){ 38 rep(i,1,n) scanf("%lf%lf",&s[i].x,&s[i].y); s[n+1]=s[1]; 39 p[1]=(P){-inf,-inf}; p[2]=(P){-inf,inf}; p[3]=(P){inf,inf}; p[4]=(P){inf,-inf}; 40 ntot=4; p[0]=p[4]; p[5]=p[1]; 41 rep(i,1,n) cut((S){s[i],s[i+1]}); 42 if (tot==0) puts("NO"); else puts("YES"); 43 } 44 45 int main(){ 46 for (scanf("%d",&T); T--; ) scanf("%d",&n),solve(); 47 return 0; 48 }
POJ2540(求直线不等式的解集)
1 #include<cstdio> 2 #include<cmath> 3 #include<algorithm> 4 #define rep(i,l,r) for (i=l; i<=r; i++) 5 using namespace std; 6 7 const double eps=1e-8; 8 char str[10]; 9 int np,nq; 10 struct P { double x,y; }pre,cur,p[110],q[110]; 11 double det(P a,P b,P c){ return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y); } 12 13 void get(P p1,P p2,double &a,double &b,double &c){ 14 if (fabs(p1.x-p2.x)<eps) a=0,b=1; 15 else if (fabs(p1.y-p2.y)<eps) a=1,b=0; 16 else b=p2.y-p1.y,a=p2.x-p1.x; 17 c=-a*(p2.x+p1.x)/2-b*(p2.y+p1.y)/2; 18 } 19 20 void work(P p1,P p2,double a,double b,double c){ 21 double u=a*p1.x+b*p1.y+c,v=a*p2.x+b*p2.y+c; 22 if (u*v<-eps){ 23 double x=(fabs(u)*p2.x+fabs(v)*p1.x)/(fabs(u)+fabs(v)); 24 double y=(fabs(u)*p2.y+fabs(v)*p1.y)/(fabs(u)+fabs(v)); 25 q[++nq]=(P){x,y}; 26 } 27 } 28 29 void cut(double a,double b,double c){ 30 int i; nq=0; 31 rep(i,1,np){ 32 if (a*p[i].x+b*p[i].y+c>-eps) q[++nq]=p[i]; 33 else work(p[i-1],p[i],a,b,c),work(p[i],p[i+1],a,b,c); 34 } 35 swap(q,p); swap(nq,np); p[np+1]=p[1]; p[0]=p[np]; 36 } 37 38 double area(int n){ 39 int i; double ar=0; 40 rep(i,2,n-1) ar+=det(p[1],p[i],p[i+1]); 41 return fabs(ar)/2.0; 42 } 43 44 int main(){ 45 p[1]=(P){0,0}; p[2]=(P){0,10}; p[3]=(P){10,10}; p[4]=(P){10,0}; 46 np=4; p[0]=p[4]; p[5]=p[1]; 47 double ar=1.0; 48 while (scanf("%lf%lf%s",&cur.x,&cur.y,str)!=EOF){ 49 double a,b,c; get(pre,cur,a,b,c); 50 if (str[0]==‘C‘) 51 { if (a*pre.x+b*pre.y+c<-eps) a=-a,b=-b,c=-c; } 52 else 53 if (str[0]==‘H‘) 54 { if (a*cur.x+b*cur.y+c<-eps) a=-a,b=-b,c=-c;} 55 else ar=0; 56 if (fabs(ar)<eps) { puts("0.00"); continue; } 57 cut(a,b,c); ar=area(np); 58 printf("%.2f\n",ar); pre=cur; 59 } 60 return 0; 61 }
$O(n \log n)$写法：
BZOJ2618:
1 #include<cmath> 2 #include<cstdio> 3 #include<algorithm> 4 #define rep(i,l,r) for (int i=l; i<=r; i++) 5 using namespace std; 6 7 const int N=1010; 8 const double eps=1e-8; 9 struct P{ 10 double x,y; 11 P(double xx=0,double yy=0):x(xx),y(yy){} 12 }p[N],a[N]; 13 struct L{P a,b; double slop;}l[1005],q[1005]; 14 15 int n,k,cnt,tot; 16 double ans; 17 18 P operator -(P a,P b){ return P(a.x-b.x,a.y-b.y); } 19 double operator *(P a,P b){ return a.x*b.y-a.y*b.x; } 20 bool operator <(L a,L b){ return (a.slop!=b.slop) ? (a.slop-b.slop)<-eps : (a.b-a.a)*(b.b-a.a)>eps; } 21 22 P inter(L a,L b){ 23 double A1=a.b.y-a.a.y,B1=a.a.x-a.b.x,C1=(a.b.x-a.a.x)*a.a.y-(a.b.y-a.a.y)*a.a.x; 24 double A2=b.b.y-b.a.y,B2=b.a.x-b.b.x,C2=(b.b.x-b.a.x)*b.a.y-(b.b.y-b.a.y)*b.a.x; 25 return P((C2*B1-C1*B2)/(A1*B2-A2*B1),(C1*A2-C2*A1)/(A1*B2-A2*B1)); 26 } 27 28 bool jud(L a,L b,L t){ P p=inter(a,b); return (t.b-t.a)*(p-t.a)<0; } 29 30 void work(){ 31 sort(l+1,l+cnt+1); 32 int L=1,R=0; tot=0; 33 rep(i,1,cnt){ 34 if (l[i].slop!=l[i-1].slop) tot++; 35 l[tot]=l[i]; 36 } 37 cnt=tot; tot=0; 38 q[++R]=l[1]; q[++R]=l[2]; 39 rep(i,3,cnt){ 40 while (L<R && jud(q[R-1],q[R],l[i])) R--; 41 while (L<R && jud(q[L+1],q[L],l[i])) L++; 42 q[++R]=l[i]; 43 } 44 while (L<R && jud(q[R-1],q[R],q[L])) R--; 45 while (L<R && jud(q[L+1],q[L],q[R])) L++; 46 q[R+1]=q[L]; 47 rep(i,L,R) a[++tot]=inter(q[i],q[i+1]); 48 } 49 50 void getans(){ 51 if (tot<3) return; 52 a[++tot]=a[1]; 53 rep(i,1,tot) ans+=a[i]*a[i+1]; 54 ans=fabs(ans)/2; 55 } 56 57 int main(){ 58 scanf("%d",&n); 59 rep(i,1,n){ 60 scanf("%d",&k); 61 rep(j,1,k) scanf("%lf%lf",&p[j].x,&p[j].y); 62 p[k+1]=p[1]; 63 rep(j,1,k) l[++cnt].a=p[j],l[cnt].b=p[j+1]; 64 } 65 rep(i,1,cnt) l[i].slop=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x); 66 work(); getans(); printf("%.3lf\n",ans); 67 return 0; 68 }
BZOJ4445:
1 #include<cmath> 2 #include<cstdio> 3 #include<algorithm> 4 #define rep(i,l,r) for (int i=l; i<=r; i++) 5 using namespace std; 6 7 const int N=100010; 8 const double eps=1e-8; 9 struct P{ 10 double x,y; 11 P(double xx=0,double yy=0):x(xx),y(yy){} 12 }p[N],a[N]; 13 struct L{P a,b; double slop;}l[N],q[N]; 14 15 int n,k,cnt,tot; 16 double ans1,ans2; 17 18 P operator -(P a,P b){ return P(a.x-b.x,a.y-b.y); } 19 double operator *(P a,P b){ return a.x*b.y-a.y*b.x; } 20 bool operator <(L a,L b){ return (a.slop!=b.slop) ? (a.slop-b.slop)<-eps : (a.b-a.a)*(b.b-a.a)>eps; } 21 22 P inter(L a,L b){ 23 double A1=a.b.y-a.a.y,B1=a.a.x-a.b.x,C1=(a.b.x-a.a.x)*a.a.y-(a.b.y-a.a.y)*a.a.x; 24 double A2=b.b.y-b.a.y,B2=b.a.x-b.b.x,C2=(b.b.x-b.a.x)*b.a.y-(b.b.y-b.a.y)*b.a.x; 25 return P((C2*B1-C1*B2)/(A1*B2-A2*B1),(C1*A2-C2*A1)/(A1*B2-A2*B1)); 26 } 27 28 bool jud(L a,L b,L t){ P p=inter(a,b); return (t.b-t.a)*(p-t.a)<-eps; } 29 void access(int x){ for (int y=0; x; y=x,x=f[x]) splay(x),ch[x][1]=y,upd(x); } 30 void work(){ 31 sort(l+1,l+cnt+1); 32 int L=1,R=0; tot=0; 33 rep(i,1,cnt){ 34 if (abs(l[i].slop-l[i-1].slop)>eps) tot++; 35 l[tot]=l[i]; 36 } 37 cnt=tot; tot=0; 38 q[++R]=l[1]; q[++R]=l[2]; 39 rep(i,3,cnt){ 40 while (L<R && jud(q[R-1],q[R],l[i])) R--; 41 while (L<R && jud(q[L+1],q[L],l[i])) L++; 42 q[++R]=l[i]; 43 } 44 while (L<R && jud(q[R-1],q[R],q[L])) R--; 45 while (L<R && jud(q[L+1],q[L],q[R])) L++; 46 q[R+1]=q[L]; 47 rep(i,L,R) a[++tot]=inter(q[i],q[i+1]); 48 } 49 50 void get1(){ 51 ans1=0; 52 if (tot<3) return; 53 a[++tot]=a[1]; 54 rep(i,1,tot) ans1+=a[i]*a[i+1]; 55 ans1=fabs(ans1)/2; 56 } 57 58 void get2(){ 59 ans2=0; p[n+1]=p[1]; 60 rep(i,1,n+1) ans2+=p[i]*p[i+1]; 61 ans2=fabs(ans2)/2; 62 } 63 64 void ins(int i,int j){ 65 double a=p[1].y-p[i].y-p[2].y+p[j].y; 66 double b=p[2].x-p[j].x-p[1].x+p[i].x; 67 double c=p[1].x*p[2].y-p[2].x*p[1].y-p[i].x*p[j].y+p[j].x*p[i].y; 68 if (abs(a)<eps && abs(b)<eps) return; 69 if (abs(a)<eps){ 70 if (b>0) l[++cnt].a=P(1,-c/b),l[cnt].b=P(0,-c/b); 71 else l[++cnt].a=P(0,-c/b),l[cnt].b=P(1,-c/b); 72 } 73 if (abs(b)<eps){ 74 if (a>0) l[++cnt].a=P(-c/a,0),l[cnt].b=P(-c/a,1); 75 else l[++cnt].a=P(-c/a,1),l[cnt].b=P(-c/a,0); 76 } 77 if (abs(a)>eps && abs(b)>eps){ 78 if (b>0) l[++cnt].a=P(1,-a/b-c/b),l[cnt].b=P(0,-c/b); 79 else l[++cnt].a=P(0,-c/b),l[cnt].b=P(1,-a/b-c/b); 80 } 81 } 82 83 int main(){ 84 freopen("bzoj4445.in","r",stdin); 85 freopen("bzoj4445.out","w",stdout); 86 scanf("%d",&n); 87 rep(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y); 88 rep(i,2,n) ins(i,i%n+1); l[++cnt].a=p[1],l[cnt].b=p[2]; 89 rep(i,1,cnt) l[i].slop=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x); 90 work(); get1(); get2(); printf("%.4lf\n",ans1/ans2); 91 return 0; 92 }