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

又來到了十分弱勢的計算幾何上，手中對付幾何題的武器又多了一把

本人對牛書上那個N^2演算法不感冒，便直接學習這個進階點的演算法（集訓隊朱澤園論文）

略讀一遍，英文太多十分之礙眼，加之論文講的又不太清楚，所以一度陷入迷惘

磨蹭了良久，無意中翻見一標程，竟然就是一個簡單的雙端隊列，豁然開朗，遂速切出此模板

好了，結束廢話，進入正題

演算法思想非常簡單，把直線按極角排序之後，在順序掃描的過程中求出可行域，整個過程有點類似於凸包那個掃描法

實現的時候，除了一個存當前可行域上的直線的雙端隊列之外，我還弄了個點隊列，記相鄰兩直線的交點

每次加入i號結點時，看上一個交點有沒有在i號直線的可行域內，如果沒有，就退掉上一條直線，一直這樣做（大家可以自己心裡想一下這個過程）

而如果只退隊尾，就可能判不出第一條直線是不是有用的了，所以還要類似地退隊首

需要注意的是，共線的情況是要分題目特殊處理的，有的取<e，有的卻要<-e，看允不允許共線；還有掃描完了之後，最後還得用隊首直線去退一下隊尾

還有一個小地方，我發現網上大多數選手寫極角排序都是先算極角再做，這樣不僅精度而且不太漂亮（我認為。。個人意見，僅供參考），我這裡用的是象限第一關鍵字，叉積再判來排的序，這個“象限”還可以簡化一點，兩個就夠了，因為在pi內，叉積是能判出左右來的

代碼（模板）：

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

 

很多時候，網上幾百行的恐怖代碼都是紙老虎，自己編才發現。。並不難！

學完了知識，搞幾道題，就簡略地總結一下吧：

 

pku2451

純裸。。。

 

pku3525

求凸多邊形內最大圓

二分將凸多邊形“內縮”mid，半平面交判斷是否存在可行域

這裡貼出內縮的那段代碼吧，其他的差不多

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

 

還是計算幾何害死人，雖說每個程式都沿用模板，但基本上每道題都要調蠻久，功夫還是不夠深啊，做幾何題也算是一種曆練吧！

好像有個很nb的單純型演算法也可以做，有空再去研究