半平面交可以解决一些线性规划问题和纯几何问题，比较常用于求多边形核等。
半平面交有两种写法，一种是比较直观的$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 ",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 ",ans); 67 return 0; 68 }
BZOJ4445:
1 #include<cmath> 2 #include<cstdio> 3 #include<algorithm> 4 #define A double k2=(b.b.y-b.a.y)/(b.b.x-b.a.x),b2=b.a.y-k2*b.a.x 5 #define B double k1=(a.b.y-a.a.y)/(a.b.x-a.a.x),b1=a.a.y-k1*a.a.x 6 #define rep(i,l,r) for (int i=l; i<=r; i++) 7 using namespace std; 8 9 const int N=100010; 10 const double eps=1e-8; 11 struct P{ 12 double x,y; 13 P(double xx=0,double yy=0):x(xx),y(yy){} 14 }p[N],a[N]; 15 struct L{P a,b; double slop;}l[N],q[N]; 16 17 int n,k,cnt,tot; 18 double ans1,ans2; 19 20 P operator -(P a,P b){ return P(a.x-b.x,a.y-b.y); } 21 double operator *(P a,P b){ return a.x*b.y-a.y*b.x; } 22 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; } 23 24 P inter(L a,L b){ 25 if (a.a.x==a.b.x){ A; return (P){a.a.x,k2*a.a.x+b2}; } 26 if (b.a.x==b.b.x){ B; return (P){b.a.x,k1*b.a.x+b1}; } 27 A; B; double x=(b2-b1)/(k1-k2),y=k1*x+b1; 28 return (P){x,y}; 29 } 30 31 bool jud(L a,L b,L t){ P p=inter(a,b); return (t.b-t.a)*(p-t.a)<-eps; } 32 33 void work(){ 34 sort(l+1,l+cnt+1); 35 int L=1,R=0; tot=0; 36 rep(i,1,cnt){ 37 if (abs(l[i].slop-l[i-1].slop)>eps) tot++; 38 l[tot]=l[i]; 39 } 40 cnt=tot; tot=0; 41 q[++R]=l[1]; q[++R]=l[2]; 42 rep(i,3,cnt){ 43 while (L<R && jud(q[R-1],q[R],l[i])) R--; 44 while (L<R && jud(q[L+1],q[L],l[i])) L++; 45 q[++R]=l[i]; 46 } 47 while (L<R && jud(q[R-1],q[R],q[L])) R--; 48 while (L<R && jud(q[L+1],q[L],q[R])) L++; 49 q[R+1]=q[L]; 50 rep(i,L,R) a[++tot]=inter(q[i],q[i+1]); 51 } 52 53 void get1(){ 54 ans1=0; 55 if (tot<3) return; 56 a[++tot]=a[1]; 57 rep(i,1,tot) ans1+=a[i]*a[i+1]; 58 ans1=fabs(ans1)/2; 59 } 60 61 void get2(){ 62 ans2=0; p[n+1]=p[1]; 63 rep(i,1,n+1) ans2+=p[i]*p[i+1]; 64 ans2=fabs(ans2)/2; 65 } 66 67 void ins(int i,int j){ 68 double a=p[1].y-p[i].y-p[2].y+p[j].y; 69 double b=p[2].x-p[j].x-p[1].x+p[i].x; 70 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; 71 if (abs(a)<eps && abs(b)<eps) return; 72 if (abs(a)<eps){ 73 if (b>0) l[++cnt].a=P(1,-c/b),l[cnt].b=P(0,-c/b); 74 else l[++cnt].a=P(0,-c/b),l[cnt].b=P(1,-c/b); 75 } 76 if (abs(b)<eps){ 77 if (a>0) l[++cnt].a=P(-c/a,0),l[cnt].b=P(-c/a,1); 78 else l[++cnt].a=P(-c/a,1),l[cnt].b=P(-c/a,0); 79 } 80 if (abs(a)>eps && abs(b)>eps){ 81 if (b>0) l[++cnt].a=P(1,-a/b-c/b),l[cnt].b=P(0,-c/b); 82 else l[++cnt].a=P(0,-c/b),l[cnt].b=P(1,-a/b-c/b); 83 } 84 } 85 86 int main(){ 87 freopen("bzoj4445.in","r",stdin); 88 freopen("bzoj4445.out","w",stdout); 89 scanf("%d",&n); 90 rep(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y); 91 rep(i,2,n) ins(i,i%n+1); l[++cnt].a=p[1],l[cnt].b=p[2]; 92 rep(i,1,cnt) l[i].slop=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x); 93 work(); get1(); get2(); printf("%.4lf ",ans1/ans2); 94 return 0; 95 }