【pku2451+3525+3384+1755】半平面交O(nlogn)

又来到了十分弱势的计算几何上，手中对付几何题的武器又多了一把

本人对牛书上那个N^2算法不感冒，便直接学习这个高级点的算法（集训队朱泽园论文）

略读一遍，英文太多十分之碍眼，加之论文讲的又不太清楚，所以一度陷入迷惘

磨蹭了良久，无意中翻见一标程，竟然就是一个简单的双端队列，豁然开朗，遂速切出此模板

好了，结束废话，进入正题

算法思想非常简单，把直线按极角排序之后，在顺序扫描的过程中求出可行域，整个过程有点类似于凸包那个扫描法

实现的时候，除了一个存当前可行域上的直线的双端队列之外，我还弄了个点队列，记相邻两直线的交点

每次加入i号结点时，看上一个交点有没有在i号直线的可行域内，如果没有，就退掉上一条直线，一直这样做（大家可以自己心里想一下这个过程）

而如果只退队尾，就可能判不出第一条直线是不是有用的了，所以还要类似地退队首

需要注意的是，共线的情况是要分题目特殊处理的，有的取<e，有的却要<-e，看允不允许共线；还有扫描完了之后，最后还得用队首直线去退一下队尾

还有一个小地方，我发现网上大多数选手写极角排序都是先算极角再做，这样不仅精度而且不太漂亮（我认为。。个人意见，仅供参考），我这里用的是象限第一关键字，叉积再判来排的序，这个“象限”还可以简化一点，两个就够了，因为在pi内，叉积是能判出左右来的

代码（模板）：

program syj; {s&i alogrithm} type real=double; const e=1e-6; maxn=20000+5; var n,i,st,ed:longint; x,y,xx,yy,xi,a,b,c,px,py,d:array[0..maxn]of real; l:array[1..maxn]of longint; ans:real; function cj(x1,y1,x2,y2,x3,y3:real):real; begin cj:=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1); end; function first(i,j:longint):boolean; begin first:=(xi[i]-xi[j]>e)or (abs(xi[i]-xi[j])<e)and(cj(x[i],y[i],xx[i],yy[i],xx[j]-x[j]+x[i],yy[j]-y[j]+y[i])>e) end; procedure swap(i,j:longint); var z:real; begin z:=x[i];x[i]:=x[j];x[j]:=z;z:=y[i];y[i]:=y[j];y[j]:=z; z:=xx[i];xx[i]:=xx[j];xx[j]:=z;z:=yy[i];yy[i]:=yy[j];yy[j]:=z; z:=xi[i];xi[i]:=xi[j];xi[j]:=z; end; procedure sort(l,r:longint); var i,j,k:longint; begin i:=l;j:=r;k:=(l+r) div 2; xi[0]:=xi[k];x[0]:=x[k];y[0]:=y[k];xx[0]:=xx[k];yy[0]:=yy[k]; repeat while first(i,0) do inc(i); while first(0,j) do dec(j); if i<=j then begin swap(i,j);inc(i);dec(j); end; until i>j; if l<j then sort(l,j); if i<r then sort(i,r); end; function f(i:longint;x,y:real):boolean; begin f:=a[i]*x+b[i]*y+c[i]<-e; end; function mj(i,j:longint):real; begin mj:=px[i]*py[j]-px[j]*py[i]; end; procedure jiaodian(a1,b1,c1,a2,b2,c2:real;var xx,yy:real); var s,x,y:real; begin s:=a2*b1-a1*b2;x:=c1*b2-c2*b1;y:=a1*c2-a2*c1; xx:=x/s;yy:=y/s; end; begin assign(input,'input.txt');reset(input); assign(output,'output.txt');rewrite(output); readln(n); for i:=1 to n do begin readln(x[i],y[i],xx[i],yy[i]);xi[i]:=ord((xx[i]-x[i]>0)or(xx[i]-x[i]=0)and(yy[i]-y[i]<0)); end; x[n+1]:=0;y[n+1]:=0;xx[n+1]:=10000;y[n+1]:=0;xi[n+1]:=1; x[n+2]:=10000;y[n+2]:=0;xx[n+2]:=10000;yy[n+2]:=10000;xi[n+2]:=0; x[n+3]:=10000;y[n+3]:=10000;xx[n+3]:=0;yy[n+3]:=10000;xi[n+3]:=0; x[n+4]:=0;y[n+4]:=10000;xx[n+4]:=0;yy[n+4]:=0;xi[n+4]:=1; n:=n+4; sort(1,n); for i:=1 to n do begin a[i]:=y[i]-yy[i];b[i]:=xx[i]-x[i];c[i]:=x[i]*yy[i]-xx[i]*y[i]; if abs(b[i])>e then d[i]:=c[i]/abs(b[i]) else d[i]:=c[i]/abs(a[i]); end; st:=1;ed:=1;l[1]:=1; for i:=2 to n do begin while (st<ed) and f(i,px[ed-1],py[ed-1]) do dec(ed); while (st<ed) and f(i,px[st],py[st]) do inc(st); inc(ed);l[ed]:=i; if abs(a[i]*b[l[ed-1]]-a[l[ed-1]]*b[i])<e then begin dec(ed);if d[i]<d[l[ed]] then l[ed]:=i; end; if ed>st then jiaodian(a[l[ed]],b[l[ed]],c[l[ed]],a[l[ed-1]],b[l[ed-1]],c[l[ed-1]],px[ed-1],py[ed-1]); end; while (st<ed) and f(l[st],px[ed-1],py[ed-1]) do dec(ed); ans:=0; if ed-st>1 then begin jiaodian(a[l[ed]],b[l[ed]],c[l[ed]],a[l[st]],b[l[st]],c[l[st]],px[ed],py[ed]);px[ed+1]:=px[st];py[ed+1]:=py[st]; for i:=st to ed do ans:=ans+mj(i,i+1); end; writeln(ans/2:0:1); close(input); close(output); end.

 

很多时候，网上几百行的恐怖代码都是纸老虎，自己编才发现。。并不难！

学完了知识，搞几道题，就简略地总结一下吧：

 

pku2451

纯裸。。。

 

pku3525

求凸多边形内最大圆

二分将凸多边形“内缩”mid，半平面交判断是否存在可行域

这里贴出内缩的那段代码吧，其他的差不多

function calc(i:longint):real; begin if x[i]<>x[i+1] then calc:=mid*dis(i,i+1)/abs(x[i]-x[i+1]) else calc:=mid; end;

a[i]:=y[i]-y[i+1];b[i]:=x[i+1]-x[i];c[i]:=x[i]*y[i+1]-x[i+1]*y[i]; if abs(b[i])>e then c[i]:=c[i]+o[ord(b[i]<0)]*b[i]*calc(i) else c[i]:=c[i]+o[ord(a[i]<0)]*a[i]*calc(i); if abs(b[i])>e then d[i]:=c[i]/abs(b[i]) else d[i]:=c[i]/abs(a[i]);

啊，这里其实还是编丑了，我用相似三角形算的。。。

 

pku3384

题意是用两个给定半径的圆覆盖一个多边形，问最多能覆盖多边形的面积

这题就每条直线内缩R，求个最远点对就可以了

代码差不多，就不贴了

p.s discuss里说spj很诡异，但可能我比较走运吧，ms没出什么精度啊，输出顺序之类的乱七八糟错误！HAPPY

 

pku1755

这是我所做的几道题里面最巧妙的一道了，参考了一下discuss里的某神牛做法

首先把让i号人获胜转化为n-1个不等式组是很显然的，但三个参量一般的线性规划做不了，怎么办？

记a,b,c为三段的长度，s=a+b+c，先同除s，又(a+b+c)/s=1，所以可以消掉c，就可以用半平面交的算法了！

需要注意的是，记得要添加3条直线,a>0,b>0,a+b<1

无聊，贴个代码

program syj; {s&i alogrithm} type real=extended; const e=1e-12; maxn=200+5; var n,i,j,k,m,st,ed:longint; x,y,xx,yy,xi,a,b,c,px,py,d:array[0..maxn]of real; l,s1,s2,s3:array[1..maxn]of longint; ans:real; function cj(x1,y1,x2,y2,x3,y3:real):real; begin cj:=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1); end; function first(i,j:longint):boolean; begin first:=(xi[i]-xi[j]>e)or (abs(xi[i]-xi[j])<e)and(cj(x[i],y[i],xx[i],yy[i],xx[j]-x[j]+x[i],yy[j]-y[j]+y[i])>e) end; procedure swap(i,j:longint); var z:real; begin z:=x[i];x[i]:=x[j];x[j]:=z;z:=y[i];y[i]:=y[j];y[j]:=z; z:=xx[i];xx[i]:=xx[j];xx[j]:=z;z:=yy[i];yy[i]:=yy[j];yy[j]:=z; z:=xi[i];xi[i]:=xi[j];xi[j]:=z; end; procedure swap(var x,y:real); var z:real; begin z:=x;x:=y;y:=z; end; procedure sort(l,r:longint); var i,j,k:longint; begin i:=l;j:=r;k:=(l+r) div 2; xi[0]:=xi[k];x[0]:=x[k];y[0]:=y[k];xx[0]:=xx[k];yy[0]:=yy[k]; repeat while first(i,0) do inc(i); while first(0,j) do dec(j); if i<=j then begin swap(i,j);inc(i);dec(j); end; until i>j; if l<j then sort(l,j); if i<r then sort(i,r); end; function f(i:longint;x,y:real):boolean; begin f:=a[i]*x+b[i]*y+c[i]<e; end; function mj(i,j:longint):real; begin mj:=px[i]*py[j]-px[j]*py[i]; end; procedure jiaodian(a1,b1,c1,a2,b2,c2:real;var xx,yy:real); var s,x,y:real; begin s:=a2*b1-a1*b2;x:=c1*b2-c2*b1;y:=a1*c2-a2*c1; xx:=x/s;yy:=y/s; end; function f(i:longint;x:real):real; begin if abs(b[i])>e then f:=-(a[i]*x+c[i])/b[i] else f:=a[i]; end; begin assign(input,'input.txt');reset(input); assign(output,'output.txt');rewrite(output); readln(m); for i:=1 to m do readln(s1[i],s2[i],s3[i]); if m=1 then begin writeln('Yes'); close(input);close(output); halt; end; for j:=1 to m do begin n:=0; for k:=1 to m do if j<>k then begin inc(n);a[n]:=1/s1[k]-1/s3[k]-1/s1[j]+1/s3[j]; b[n]:=1/s2[k]-1/s3[k]-1/s2[j]+1/s3[j];c[n]:=1/s3[k]-1/s3[j]; if abs(b[n])>e then begin x[n]:=0;y[n]:=-c[n]/b[n]; xx[n]:=1;yy[n]:=-(a[n]+c[n])/b[n]; if b[n]<0 then begin swap(x[n],xx[n]);swap(y[n],yy[n]); end; d[n]:=c[n]/abs(b[n]); end else if abs(a[n])<e then begin dec(n); if c[n]<e then begin n:=0;break; end end else begin x[n]:=-c[n]/a[n];y[n]:=0; xx[n]:=-c[n]/a[n];yy[n]:=1; if a[n]>0 then swap(y[n],yy[n]); d[n]:=c[n]/abs(a[n]); end; end; if n=0 then begin writeln('No');continue; end; x[n+1]:=0;y[n+1]:=0;xx[n+1]:=1;yy[n+1]:=0; x[n+2]:=0;y[n+2]:=1;xx[n+2]:=0;yy[n+2]:=0; x[n+3]:=1;y[n+3]:=0;xx[n+3]:=0;yy[n+3]:=1; n:=n+3; for i:=1 to n do xi[i]:=ord((xx[i]-x[i]>e)or(abs(xx[i]-x[i])<e)and(yy[i]-y[i]<-e)); sort(1,n); for i:=1 to n do begin a[i]:=y[i]-yy[i];b[i]:=xx[i]-x[i];c[i]:=x[i]*yy[i]-xx[i]*y[i]; if abs(b[i])>e then d[i]:=c[i]/abs(b[i]) else d[i]:=c[i]/abs(a[i]); end; st:=1;ed:=1;l[1]:=1; for i:=2 to n do begin while (st<ed) and f(i,px[ed-1],py[ed-1]) do dec(ed); while (st<ed) and f(i,px[st],py[st]) do inc(st); inc(ed);l[ed]:=i; if abs(a[i]*b[l[ed-1]]-a[l[ed-1]]*b[i])<e then begin dec(ed);if d[i]<d[l[ed]] then l[ed]:=i; end; if ed>st then jiaodian(a[l[ed]],b[l[ed]],c[l[ed]],a[l[ed-1]],b[l[ed-1]],c[l[ed-1]],px[ed-1],py[ed-1]); end; while (st<ed) and f(l[st],px[ed-1],py[ed-1]) do dec(ed); if ed-st<=1 then writeln('No') else begin ans:=0; jiaodian(a[l[ed]],b[l[ed]],c[l[ed]],a[l[st]],b[l[st]],c[l[st]],px[ed],py[ed]);px[ed+1]:=px[st];py[ed+1]:=py[st]; for i:=st to ed do ans:=ans+mj(i,i+1); if ans>1e-16 then writeln('Yes') else writeln('No'); end; end; close(input); close(output); end.

 

还是计算几何害死人，虽说每个程序都沿用模板，但基本上每道题都要调蛮久，功夫还是不够深啊，做几何题也算是一种历练吧！

好像有个很nb的单纯型算法也可以做，有空再去研究