關於半平面交的探討
參考文獻
你穀日報真的是太好啦:
半平面交中雙端佇列的一種錯誤寫法的講解:https://www.luogu.com.cn/blog/wjyyy/geometry1
本文會有許多99%和本文相似的地方,因為就是從這個日報才終於看懂了半平面交,對於半平面交的理解才開始加深。
除小部分程式碼其餘皆採用此部落格的程式碼,圖片也大部分是基本上只要是好看的都是,好像只要看水印就行了。
我以前的證明都是什麼垃圾啊,嚴謹又不嚴謹,就是一坨*
當然,本文的證明還是更加的偏感性,因為計算幾何理性證明是在太難了QAQ,為什麼我不會理性的計算幾何啊(╯‵□′)╯︵┻━┻。
題目
【題意】 在一個有限大(-10 0000<=x,y<=10 0000)的平面座標系上有n個半平面(注意有限的),每個半平面給出一條有向線段(x1,y1)——>(x2,y2)。 每個半平面的有效區域都是左側(包括此直線)。求這n個半平面的交的面積。 【輸入檔案】 第一行一個整數n(1 <= n <= 2 0000) 下來n個半平面。每個半平面一行四個整數x1,y1,x2,y2。(-100000 <= x1,y1,x2,y2 <= 100000) 【輸出檔案】 一行,輸出n個半平面的交的面積(保留一位小數),如果有效面積不存在則輸出"0.0"。 【樣例1輸入】 3 1 1 9 9 5 10 4 10 0 3 0 2 【樣例1輸出】 50.0 【樣例2輸入】 4 0 0 10 0 10 0 10 10 10 10 0 10 11 11 11 0 【樣例2輸出】 0.0
半平面交中的直線更像是有向直線,你可以理解為把向量擴充套件到了直線。
而有向直線的極角則為向量的幅角。
考慮\(ax+by+c≥0\),首先有\(ax+by+c=0\),如果\(a>0\),那麼\(x'=x+k(k≥0)\)都滿足此要求,此時解集是在這條直線的右邊,如果\(a<0\)是在左邊,\(a=0\)則要看\(b\),此時就是上下之分了,因此,有一條直線的半平面更多的是這個直線的一個二元一次不等式的體現,而其解集就是其的半平面,這裡預設半平面包括這條直線。
但是,這裡卻說是有向直線的左邊,有向直線和直線不同,多了方向,因此其左右邊更多的清晰,即把這個平面旋轉,旋轉到直線與轉前的\(y\)
因此,\(x+y≥0\)和\(-x-y≥0\)是同一條直線,但是方向卻不同,至於如何把二元一次不等式轉換成有向直線,下面的練習會講,這裡先不講。
半平面交求的是多條直線的半平面的重合部分,也就是多個二元一次方程的解集。
這裡講的無交情況總共有兩種,一種是解集為空,一種是解集無限。(當然,大部分的題目都會讓你在一個有限的正方形或者有限的範圍跑半平面交,下文會大概的提一提這種情況如何處理,但是呢,其不是主要無交情況,程式碼中主要也不會考慮到這種情況,請大家閱讀時發現無限交沒有提到也不要吃驚,至於如何限定在一個範圍?如果限定在一個正方形,就把正方形的四條邊以有向直線的方式加入到半平面交即可,其餘類似)
當然,因為有向線段的左邊更加直觀,所以一般在求解的過程中都是儲存代表有向直線的向量,而且這能更好的利用叉積。
求法
約定
以下來自https://www.luogu.com.cn/blog/105254/dui-ban-ping-mian-jiao-suan-fa-zheng-que-xing-xie-shi-di-tan-suo。
程式碼約定:
typedef pair<int, int> pad;/*pad awa*/
const double pi =acos(-1), eps =1e-6;
/*點或者說向量*/
struct vect{
double x, y;
vect(){}
vect(double xx, double yy):x(xx), y(yy){}
vect operator + (vect v){ return vect(x+v.x, y+v.y); }
vect operator - (vect v){ return vect(x-v.x, y-v.y); }
vect operator * (double mu){ return vect(mu*x, mu*y); }
double operator / (vect v){ return x*v.y-y*v.x; }/*叉積*/
};
/*直線*/
struct line{
vect u, v;
double angle;
line(){}
line(vect uu, vect vv):u(uu), v(vv){ angle =atan2(vv.y-uu.y, vv.x-uu.x); }
};
line hull[];/*儲存 凸包 / 環(定義見下) / 答案*/
line ls[];/*代表加入直線序列*/
int l, r;/*左右指標,左閉右開*/
/*(a > b)*/
inline bool gtr(double a, double b){ return (a-b > eps); }
/*(a == b)*/
inline bool eq(double a, double b){ return (a-b < eps && a-b > -eps); }
/*點是否在有向直線右側,是返回 true,否(在左側或直線上)返回 false*/
inline bool onright(line f, vect w){ return (gtr((w-f.u)/(f.v-f.u), 0)); }
/*求兩條直線交點(無交情況未定義)*/
vect getIntersection(line f, line g){
double w =((g.u-f.u)/(f.u-f.v))/((f.u-f.v)/(g.u-g.v));
return g.u+(g.u-g.v)*w;
}
/*用於排序的比較函式*/
int cmp(line A, line B){
if(eq(A.angle, B.angle)) return onright(B, A.u);/*有向直線最左的會排在最後面,並被保留*/
else return (gtr(B.angle, A.angle));
}
說明:
- 文中的半平面交預設是有向直線左側交,包括邊界。
- 實現中所用的計較都是直接使用
atan2
的返回值,並未做其他處理(因此值域為 \((−π,π]\))、 - 這裡的增量法實現都是按直線極角從小到大加入的。
- 話說下面
以及上面提到的 "增量法"、"標準雙端佇列做法",一種典型示範就是劉汝佳書中有關半平面交章節的示例程式碼。或者關於這種做法的實現細節及說明,可以先看看之前的 計算幾何日報。當然,本文的程式碼來自此篇日報
如何求一個點是否在其半平面
額,其實你會發現一個事情,對於有向直線\(a->b\),存在點\(c\)在其半平面,那麼\(\vec{ab}*\vec{ac}≥0\),這是判斷是否在半平面的充要條件,當然,上面約定函式中看是否在直線右邊也是同樣的原理。
簡略求法
現在,我們考慮把所有的直線按照極角排序,什麼?你問我極角怎麼求?已知有向直線\(AB\),其極角為\(atan2(y_B-y_A,x_B-x_A)\)。
其實基本上很多也就半平面交和凸包好吧幾何題都可以用極角排序的方式確定遍歷直線的順序,設排序後的陣列為\(a\)陣列。
不妨先講個粗略的求法,最後一步步彌補缺漏,最後完成程式碼。
因為面積有限,不難先猜想其是個多個多邊形,後面發現不能是非凸多邊形(自己畫一個就知道了),於是又發現其必須只能是一個(自己畫以下),所以其是一個凸包(這也就可以解釋為什麼前面要用極角排序了,為了方便下文求凸包)。
我們先講個類似凸包的做法,看看大家有沒有異議,用棧(棧用\(A\)代表)儲存與半平面交部分有直接接觸的有向直線(也就是構成半平面交輪廓的直線),對於前兩條直線,直接儲存(紅色半平面和藍色半平面)。
現在已經有了\(top\)條直線,對於第\(i\)(這裡\(i=3\))條有向直線(綠色半平面)
如果\(A_{top}\)和\(A_{top-1}\)的交點在\(a_i\)的半平面,那麼說明直接插入即可,如果剛好在直線上也是一樣,如果不在,那麼說明什麼,如圖,你會發現藍色半平面的有向(這裡根本看不出有向好吧,看半平面吧)直線根本和半平面交部分沒有任何接觸,此時我們稱其為無用直線,直接彈出,然後接著判斷。
然後這樣不斷的判斷下去,最終直線和其相鄰間的交點便會長的像凸包一樣。
這就是粗略的做法,但是其是極其不完善的,甚至普通情況都難以處理,下面開始講解不同情況以及其處理手法。
對於棧中元素更加形象的儲存和講解方式
哇,日報中講的模型是真的直觀易懂。
這裡先說以下,如果\(A_{i}\) ~ \(A_{j}\)都交於一個點的話,那麼只保留\(A_{i}\),\(A_{j}\),方便下文討論,在你看懂之後,再將其加入回來再思考一波正確性,你會發現其完全不會影響我們的討論,可以讓討論存在的爭議更少。
我們先將存下的直線想象成一條由一段段線段組成的 "鏈"。
或者具體地來講,是:
鏈由一組有向線段 \(v_1, v_2, v_3, ..., v_n\) 組成,並滿足 \(v_1\) 的結尾點是 \(v_2\) 的起始點, \(v_2\) 的結尾點是 \(v_3\) 的起始點,...,$ v_{n-1}$ 的結尾點是 \(v_n\) 的起始點,同時這個序列的極角遞增(這裡也可以反著定義遞減,下面的說明類似),且滿足最後一個元素的極角與第一個元素的極角差不超過 \(2π\)(弧度制),這裡並不明確頭尾元素的長度。
而這裡,除了\(v_{1}\)和\(v_{n}\)以外的點都是相鄰直線的交點,而其的有向線段就是有向直線與半平面交接觸的部分。
而\(v_1,v_n\)嗎,還記得最後一句話嗎?這裡並不明確頭尾元素的長度,這個是什麼意思?因為我們只是處理了每條直線和\(top\)以及\(top-1\)的關係,而沒有處理頭尾元素的關係,所以我們認為首尾元素與半平面交的部分接觸的長度是無限的,因此,我們可以認為,\(v_{1}\)和\(v_{2}\)所代表的有向直線是\(A_{1}\),\(v_{n-1}\)和\(v_{n}\)所代表的有向直線是\(A_{top}\),而\(v_{1}\)和\(v_{n-1}\)則是滿足上面兩個要求的任意點。
我們成功的把半平面交的模型轉換成容易理解的長得像凸包的模型。
環的出現及處理
在跑完半平面後,若是有交情況,就要考慮環的處理了。
而其說的環是什麼呢?
額,我猜日報所說的環分三種:
-
原本模型中的有向線段(不包括\(v_{1}->v_{2}\)和\(v_{n-1}->v_{n}\))的出現了交點,存在環,即下圖的情況四。
-
原本的有向線段不存在交點,但棧底或者棧頂的直線與有向線段(不包括\(v_{1}->v_{2}\)和\(v_{n-1}->v_{n}\))存在交點,即為情況\(2,3\)
-
以上兩種情況都不存在,且棧底和棧頂直線存在交點。
總之你會發現,就是棧頂和棧底這兩種情況瘋狂的出問題。
其實我們會發現,情況\(1\)就是我們想要的結果,但是\(2,3,4\)都出現了無用的有向線段(也就是無用的有向直線,你是可以證明這些無用的有向線段所代表的直線也是無用的)。
這個時候我們發現,棧底的元素也是可能存在問題的,不妨將其換為雙端佇列(只不過插入時能看成棧)。
首先可以看出對於 [情況2], [情況3]
,只要判斷線段的交點在不在 隊尾、隊頭 元素的右側即可。
/*情況2*/
while(l < r-1 && onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
/*情況3*/
while(l < r-1 && onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;
但是對於情況4呢?
其實你會發現:
onright(hull[r-1], getIntersection(hull[l], hull[l+1]))
onright(hull[l], getIntersection(hull[r-1], hull[r-2]))
//即上文的兩個條件
情況4肯定滿足其中一個。
你肯定會問,憑啥?
我們不妨把整個平面旋轉,把\(hull[l]\)變成\(x\)軸,方向為\(x\)負半軸的方向,這樣保證了\(hull[l]\)和\(hull[l+1]\)的交點是在\(x\)軸上,且\(hull[l]\)的半平面為\(y≤0\),更加的方便討論。
不妨通過構造反例來實現,首先\(hull[r]\)與\(hull[r-1]\)的交點必須在\(y≤0\)的部分。
不難發現,\(hull[r]\)的極角絕對大於\(0\)好像這個沒用,且無論怎麼拖動被圈住的點,都一定滿足\(onright(hull[r-1], getIntersection(hull[l], hull[l+1]))\)
當然,你會發現嚴謹的證明我完全不會,但是感性的證明還是算比較直觀的吧。
因此,你只要改造一下程式碼就能處理所有的情況。
while(l < r-1){
if(onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
else if(onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;
else break;/*已經沒有更新了*/
}
至此,如果有交,則已經處理完畢了。
無交情況判斷
先預設無平行線段(即幅角之差的絕對值\(=180°\)和幅角相等的情況),後面講。
無限交
無限交的情況是什麼?
首先先講一定有解的情況,還記得凸包的充要條件嗎?
沒錯,就是滿足那東西,還記得凸包的內角一定小於\(180°\)嗎?所以相鄰兩條直線的幅角差的絕對值要小於\(180°\),頭尾的元素也是如此。
如果出現了\(≥\)的情況,說明交是無限交,那麼如何判定此情況呢?
上圖不滿足極角排序
管他,旋轉一下不就滿足了
此時你會發現,其會在下述判定中直接把佇列變成兩個元素。
while(l < r-1){
if(onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
else if(onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;
else break;/*已經沒有更新了*/
}
因此,只要最後判斷佇列是否中小於等於兩個元素即可(那有沒有可能小於等於兩個元素還有有限交的爬,都說了只是去掉無效元素,有效元素小於等於2個怎麼可能有限交)。
但是如果滿足了佇列中相鄰兩條直線的幅角差的絕對值要小於\(180°\)就一定有交嗎?
廢話,充要條件啊(╯‵□′)╯︵┻━┻。
當然,有時候還會出現頭尾的情況,因此還要判斷頭尾的幅角差,不過要注意還要模\(360°\)(弧度制模\(2\pi\)),不過一般不用判斷,因為無限交的情況一般不會出現,因為會在一定的範圍內(那為什麼還要討論有限交的情況,因為下文中的無交就是轉換成這種情況,只不過不會轉換乘頭尾的情況罷了)。
無交的情況
首先,如果存在無交的情況,則一定存在一條\(a_{i}\),在其插入的時候,其半平面與當時的半平面交的交集為空(充要條件)。
即存在一條直線能夠直接把佇列踢到\(l=r\)(必要條件)。
好,接下來較為理性的證明,存在一條\(a_{i}\),在其插入的時候,其半平面與當時的半平面交的交集為空的充要條件是,\(a_{i}\)與\(hull[l]\)的幅角差的絕對值\(≥180°\)。
對於\(hull[l]\),傳統異能,旋轉到\(x\)軸,方向負半軸,那麼對於幅角\(<0°\)的直線而言,其和\(hull[l]\)的半平面的交集都一定包含點\((INF,0)\),所以其半平面交絕對不為空。
這樣子的話,插入就一定存在兩條相鄰直線幅角差的絕對值大於等於\(180°\),轉換為上面的情況,直接判斷即可。
平行的情況
為什麼要特別的拿出來講一下,
因為你還記得函式中有個求交點的嗎?交點是處理不了平行的。
對於幅角相等的兩條直線,如果不共線,一個必然淘汰另外一個,然後就不用求交點了,但是共線呢?
因此,半平面交還有一個騷操作,在半平面交極角排序時,預設半平面比較大的優先。
int cmp(line A, line B){
if(eq(A.angle, B.angle)) return onright(B, A.u);/*有向直線最左的會排在最後面,並被保留*/
else return (gtr(B.angle, A.angle));
}
然後在for迴圈插入中加入此語句。
while(i < totl-1 && eq(ls[i].angle, ls[i+1].angle)) ++i;/*去重*/
對於幅角之差的絕對值\(=180°\)的情況。
額。
如果有一條直線在這兩條直線之間,那麼佇列中不相鄰,也就不用考慮他們兩個交點情況。
但是對於沒有直線的情況呢?那麼這個時候就會出現兩種情況。
- 無限交,這個時候只要在入隊時暴力特判\(hull[r]\)是否是平行的即可。
- 無交,這個時候,就是下圖的情況了。
這個時候,不難發現的事情是,佇列會一直踢出只剩隊頭,所以只需要在入隊時判斷\(hull[r]\)是否是平行的即可。
當然,一般不存在無限交的情況,因為一般會限定在一個正方形內。
解集為點或者線段
還記得之前我說的如果交點在直線\(a_{i}\)上,不彈出棧嗎?就是為了應付這種情況的,如果在直線上彈出,那麼就有可能導致這些情況被誤判成無解。
但是在上述的證明中,又有可能會被理解為多條有向線段交於同一個點(實際上有向線段長度為\(0\)),因此交大家忽略掉\(A_{l+1}\) ~ \(A_{r-1}\)的直線,但是如果你懂了證明,你會發現這些直線根本不會影響證明,刪除與否根本無關,只不過刪除了更加好的表述而已。
因此,你會發現,其實如果解集是點和線段可以視作無解的話,那麼可以直接認為在交點在直線上可以直接彈出。(其實等價於認為半平面不包括直線。)
程式碼
時間複雜度:\(O(nlogn)\)。
/*這裡求得的交包括邊界,主要取決於 onright() 函式*/
inline pad getHPI(line ls[], int totl, line hull[]){
sort(ls, ls+totl, cmp);
int l =0, r =0;
for(int i =0; i < totl; ++i){
while(i < totl-1 && eq(ls[i].angle, ls[i+1].angle)) ++i;/*去重*/
while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
if(r > 0/*首次迴圈可能為空*/ && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);/*判平行的無交*/
hull[r++] =ls[i];
}
while(r-l > 1){
if(onright(hull[r-1], getIntersection(hull[l], hull[l+1]))) ++l;
else if(onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;
else break;/*已經沒有更新了*/
}
if(r-l < 3) return pad(0, 0);
else return pad(l, r);/*返回答案所在區間*/
}
另外一種雙端佇列的寫法
講解
這實際上更常用,但是上面一份程式碼的講解更加直觀。(好像上面的做法叫標準增量法,這個叫標準雙端佇列做法)
如果我們在插入的時候不是利用類似棧的形式,而是充分發揮佇列的優勢呢?
什麼意思呢?為什麼會產生第\(4\)種情況?如果我們在他還是第\(2\)種情況時就馬上變成第一種情況,那麼他最後不就只剩下第\(1,3\)種情況了嗎?因此由此產生了在判斷完隊尾後判斷隊首的情況,最後只需要判斷第\(3\)種情況即可。
/*這裡求得的交包括邊界,主要取決於 onright() 函式*/
inline pad getHPI(line ls[], int totl, line hull[]){
sort(ls, ls+totl, cmp);
int l =0, r =0;
for(int i =0; i < totl; ++i){
while(i < totl-1 && eq(ls[i].angle, ls[i+1].angle)) ++i;/*去重*/
while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
while(r-l > 1 && onright(ls[i], getIntersection(hull[l], hull[l+1]))) ++l;
if(r-l > 0/*首次迴圈可能為空*/ && eq(ls[i].angle-hull[r-1].angle, pi)) return pad(0, 0);/*判平行的無交*/
hull[r++] =ls[i];
}
while(r-l > 2 && onright(hull[l], getIntersection(hull[r-1], hull[r-2]))) --r;/*僅彈出隊尾即可*/
if(r-l < 3) return pad(0, 0);
else return pad(l, r);/*返回答案所在區間*/
}
程式碼會出的問題
考慮調換下列程式碼的先後順序:
while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
while(r-l > 1 && onright(ls[i], getIntersection(hull[l], hull[l+1]))) ++l;
變成:
while(r-l > 1 && onright(ls[i], getIntersection(hull[l], hull[l+1]))) ++l;
while(r-l > 1 && onright(ls[i], getIntersection(hull[r-1], hull[r-2]))) --r;
這樣竟然會錯?
其實一般情況下是不會錯的,但是有一種情況會錯,就是直接把佇列彈出到只剩\(l=r\)的情況,這種情況,下面的寫法會導致錯誤,首先,很明顯的一個事情是,這樣子可能會導致無交的情況直接變成頭尾幅度差大於等於\(180°\)的無限交的情況,但是常規的程式碼一般不會判斷這種情況,固會導致錯誤。
無限交情況省略我沒細想,大致差不多吧。
此處採用wjyyy的證明:
更加一般的有限交情況,其也會在把佇列彈出到只剩\(l=r\)的情況出錯。
一般情況下,我們在佇列 (佇列順序為 \(\left\{\vec{u},\vec{v}\right\}\)) 後面加一條邊(向量 \(\vec{w}\)),會產生一個交點 \(N\),縮小 \(\vec{v}\)後面的範圍。
但是畢竟每次操作都是一般的,因此可能會有把 \(M\) 點“擠出去”的情況。
如果此時出現了向量 \(\vec{a}\) ,使得 \(M\) 在 \(\vec{a}\) 的右側,那麼 \(M\) 就要出隊了。此時如果從隊首列舉++l
,顯然是擴大了範圍。實際上 \(M\)點是由 \(\vec{u}\) 和 \(\vec{v}\) 共同構成的,因此需要考慮影響到現有程序的是 \(\vec{u}\) 還是 \(\vec{v}\)。而因為我們在極角排序後,向量是逆時針順序,所以 \(\vec{v}\) 的影響要更大一些。
就如上圖,如果 \(M\) 確認在 \(\vec{a}\) 的右側,那麼此時 \(\vec{v}\) 的影響一定不會對半平面交的答案作出任何貢獻。
而我們排除隊首的原因是當前向量的限制比隊首向量要大,這個條件的前提是佇列裡有不止兩個線段(向量),不然就會出現上面的情況。
所以一定要先排除隊尾再排除隊首。(實際上排除了隊首的情況基本上就是情況\(2,4\)了)
線上做法的idea
一下引用日報《對半平面交演算法正確性解釋的探索》的講解:
現在考慮給已知半平面交新增一條有向直線的情況。
顯然新交不包含新加入直線右側的所有點(在直線上的點保不保留沒大區別)。
同時,還發現這些點還滿足:
- 它們一定是在凸包邊界上 "連續" 的一段
- 我們將凸包所有邊按極角排序,第一個極角大於和第一個極角小於插入直線的極角的兩條邊的交點一定是需要彈出的點(即使在最壞情況下(只彈出一個點))。
(如下圖;對更多點情況可平移直線並由凸包性質發現)
實現時,我們先按 性質二 嘗試找到一個需要彈出的點。如果這個點不在加入直線的右側,則這條直線就不會對交產生影響;否則只需逆/順時針依次判斷即可(這些點是相鄰的),並對於一個方向在第一個不滿足要求(在加入直線右側)的點停下。
(至於點在當前直線上...其實這時候不特判去除也不會影響貢獻。但在下一部分的做法中,不判這種情況可能會在某些特殊情況下出錯)
同時,為了保證精度,我們仍儲存加入直線的資訊(通常用兩個點表示直線...加入後並不重求點為交點),實際的凸包點在判斷時求出。而直線刪除的條件是它的兩個交點都被彈出(一個交點被彈出的話保留,由於相鄰的直線(邊)被更新了,那個交點下次求出來也是不同的,就相當於彈出了);同時注意對可能將所有直線彈出的加入操作(即無交)需要特判。(或者乾脆存點...這段可以無視)
為了能支援二分查詢而快速找到上述兩條直線,儲存的直線需按極角排序,考慮到還可能要在序列中間插入點,我們可能需要用平衡樹來維護。
(由於在競賽沒什麼實用價值程式碼就暫時咕咕了 QAQ (話說半平面交也真的常考嗎...))
這種線上做法也可以很方便地改成離線。不過單純照搬的話其實有很多部分是可以簡化一下的。(其實這裡應該早就有人看出這種做法和標準增量法很相似了...)
題目程式碼
多少年前的程式碼了
採用標準兩端佇列做法。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 21000
#define MAXN 100000
using namespace std;
double eps=1e-8;
struct dian
{
double x,y;
dian(double xx=0.0,double yy=0.0){x=xx;y=yy;}
}pi[N];//點
inline dian operator-(dian x,dian y){return dian(x.x-y.x,x.y-y.y);}
inline dian operator+(dian x,dian y){return dian(x.x+y.x,x.y+y.y);}
inline dian operator*(dian x,double y){return dian(x.x*y,x.y*y);}
struct line
{
dian x,y;double agr;
void set(double x1,double y1,double x2,double y2){x.x=x1;x.y=y1;y.x=x2;y.y=y2;agr=atan2(y2-y1,x2-x1);/*求角度*/}
}ss[N],li[N];int head,tail,n,tp;
inline double mu(dian x1,dian x2,dian y){return (x1.x-y.x)*(x2.y-y.y)-(x2.x-y.x)*(x1.y-y.y);}//叉積
inline double muu(dian x1,dian x2){return x1.x*x2.y-x2.x*x1.y;}//叉積
inline dian jd(line x,line y)
{
double w =muu((y.y-y.x),(x.x-y.x))/muu((x.y-x.x),(y.y-y.x));
return x.x+(x.y-x.x)*w;
}
inline bool safe(dian x,line y){return mu(y.x,x,y.y)<=eps;}//在不在此半平面
inline bool cmp(line x,line y)//角度排序
{
if(x.agr<y.agr)return true;
else if(fabs(x.agr-y.agr)<=eps && safe(x.x,y)==true/*如果同一個角度那麼就判斷哪個直線的半平面比較小(在有限範圍內)*/)return true;
return false;
}
void HPI()
{
sort(ss+1,ss+n+1,cmp);
tp=1;for(int i=2;i<=n;i++)if(ss[i].agr-ss[i-1].agr>eps)ss[++tp]=ss[i];//去重
li[head=1]=ss[1];li[tail=2]=ss[2];n=tp;
for(int i=3;i<=n;i++)
{
while(head<tail && safe(jd(li[tail],li[tail-1]),ss[i])==false)tail--;//踢隊尾
while(head<tail && safe(jd(li[head],li[head+1]),ss[i])==false)head++;//踢隊頭
li[++tail]=ss[i];//加入
}
while(head<tail && safe(jd(li[tail],li[tail-1]),li[head])==false)tail--;//再來一次
if(tail-head<2){printf("0.0\n");return ;}//無解
int pl=0;for(int i=head;i<tail;i++)pi[++pl]=jd(li[i],li[i+1]);pi[++pl]=jd(li[head],li[tail]);//求出凸多邊形
double ans=0;
for(int i=1;i<=pl;i++)ans+=mu(pi[i-1],pi[i],pi[1]);//算面積
if(ans<0)ans=-ans;
printf("%.1lf\n",ans/2.0);
}
int main()
{
scanf("%d",&n);
ss[1].set(MAXN,MAXN,-MAXN,MAXN);ss[2].set(-MAXN,MAXN,-MAXN,-MAXN);
ss[3].set(-MAXN,-MAXN,MAXN,-MAXN);ss[4].set(MAXN,-MAXN,MAXN,MAXN);n+=4;//新增四條直線,因為題目添加了,所以不用考慮無限的情況
for(int i=5;i<=n;i++)
{
double x1,y1,x2,y2;scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
ss[i].set(x1,y1,x2,y2);
}
HPI();
return 0;
}
一道練手題
根據古典概形得知,概率就是\(P\)點的可行位置的面積\(÷\)操場面積。
對於凸多邊形的每條線段,推導其和\(0-1\)邊面積之間的不等式,如果是解集是個半平面,則使用半平面交(事實證明單個不等式的解集確實是半平面)。
設是\(i-(i+1)\)的邊和\(0,1\)邊進行討論,利用叉積算面積列出\(P\)的不等式:
\((x_1-x_0)(y_P-y_0)-(y_1-y_0)(x_P-x_0)≤(x_{i+1}-x_i)(y_{P}-y_i)-(y_{i+1}-y_i)(x_P-x_{i})\)
\(-(x_1-x_0)y_0+(x_1-x_0)y_P-(y_1-y_0)x_P+(y_1-y_0)x_0≤(x_{i+1}-x_i)y_{P}-(x_{i+1}-x_i)y_i-(y_{i+1}-y_i)x_P+x_{i}(y_{i+1}-y_i)\)
\((y_{i+1}-y_i-y_1+y_0)x_{P}+(x_{1}-x_0+x_{i}-x_{i+1})y_{P}+(x_0y_1-x_1y_0+x_{i+1}y_i-y_{i+1}x_{i})≤0\)
注:因為是逆時針給出點,所以叉積為正,所以直接這樣列不等式。
PS:許多人說是\(<\)號,但是我JO得說他是最小的,但是沒說最小是唯一的啊,所以個人覺得是\(≤\)這個得看對最小的定義是什麼,當然,不管\(<\)還是\(≤\),也就差直線這麼一個面積,但是直線的面積可以是為\(0\),因此怎麼看都不會錯。
不難發現,是形如:\(ax+by+c≤0\)的形式,因此是半平面,用半平面交求解集面積(當然,別忘了半平面是限定在一個凸包當中的)。
因此,我們就得到了其解析式,但是換成有向直線呢?
隨便啦,自己想象就知道啦。(看程式碼也不難理解啦)
時間複雜度:\(O(nlogn)\)。
PS:這裡點的下標從\(1\)開始,而不是從\(0\)開始。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 110000
#define NN 210000
using namespace std;
double eps=1e-6;
struct dian
{
double x,y;
dian(double xx=0,double yy=0){x=xx;y=yy;}
};
inline dian operator+(dian x,dian y){return dian(x.x+y.x,x.y+y.y);}
inline dian operator-(dian x,dian y){return dian(x.x-y.x,x.y-y.y);}
inline dian operator*(dian x,double y){return dian(x.x*y,x.y*y);}
inline double operator/(dian x,dian y){return (x.x*y.y-x.y*y.x);}//叉積
struct line
{
dian x,y;double arg;
void set(){arg=atan2(y.y-x.y,y.x-x.x);}//atan2算極角
}li[NN],st[NN]/*佇列*/;int len,head,tail;
inline dian operator*(line x,line y)//求兩線交點
{
double tt=((y.y-y.x)/(x.x-y.x))/((y.y-y.x)/(x.x-x.y));
return x.x+(x.y-x.x)*tt;
}
inline bool eq(double x,double y){return (x-y)<=eps && (x-y)>=-eps;}//判斷x==y的
inline bool safe(line x,dian y){return (x.y-x.x)/(y-x.x)>0;}//左邊為safe
inline bool cmp(line x,line y)
{
if(!eq(x.arg,y.arg))return x.arg<y.arg;
else return safe(x,y.x);
}
double BPMJ()
{
sort(li+1,li+len+1,cmp);
head=1;tail=0;
for(int i=1;i<=len;i++)
{
while(i<len && eq(li[i].arg,li[i+1].arg))i++;
while(head<tail && !safe(li[i],st[tail]*st[tail-1]))tail--;
while(head<tail && !safe(li[i],st[head]*st[head+1]))head++;
if(head<=tail && eq(li[i].arg-st[tail].arg,3.141592653))return 0;
st[++tail]=li[i];
}
while(head+1<tail && !safe(st[head],st[tail]*st[tail-1]))tail--;
if(tail-head<=1)return 0;
double ans=0;
dian ji=st[tail]*st[head];
for(int i=head+1;i<tail;i++)ans+=((st[i-1]*st[i])-ji)/((st[i]*st[i+1])-ji);
return ans/2;
}
dian sh[N];
int n;
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
scanf("%lf%lf",&sh[i].x,&sh[i].y);
}
len=n;sh[n+1]=sh[1];
for(int i=1;i<=n;i++)li[i].x=sh[i],li[i].y=sh[i+1],li[i].set();//限定在凸包範圍
double ans=0;
for(int i=2;i<n;i++)ans+=(sh[i]-sh[1])/(sh[i+1]-sh[1]);//算凸包面積
ans/=2;
//轉換部分
for(int i=2;i<=n;i++)
{
double a=(sh[i+1].y-sh[i].y-sh[2].y+sh[1].y),b=(sh[2].x-sh[1].x+sh[i].x-sh[i+1].x),c=(sh[1].x*sh[2].y-sh[2].x*sh[1].y+sh[i+1].x*sh[i].y-sh[i].x*sh[i+1].y);
if(eq(a,0) && eq(b,0))
{
if(c>=0)
{
printf("0\n");
return 0;
}
continue;
}
if(eq(b,0)==1)
{
len++;
if(a>=0)li[len].x=dian(-c/a,1),li[len].y=dian(-c/a,2);
else li[len].x=dian(-c/a,2),li[len].y=dian(-c/a,1);
li[len].set();
continue;
}
if(b<=0)
{
len++;
li[len].x=dian(1,-(a+c)/b);li[len].y=dian(2,-(2*a+c)/b);
}
else if(b>0)
{
len++;
li[len].x=dian(2,-(2*a+c)/b);li[len].y=dian(1,-(a+c)/b);
}
li[len].set();
}
ans=BPMJ()/ans;
printf("%.4lf\n",ans);
return 0;
}
小結
真噁心