- 問題模型
有四只鴨子，隨機分布在圓形的池塘中，請問四只鴨子同時處於同一個半圓的概率有多大？
- 統計模擬
使用計算機進行大量的隨機實驗，統計推斷出未知變量的概率。我們將進行100萬次獨立采樣試驗，每次試驗中讓四只鴨子隨機分布在圓形池塘中，統計四只鴨子處於同一半圓的試驗次數。那么使用該次數除以100萬，就可以大概得到題目所求的概率。
- 建立模型
假設以原點為中心，半徑為1的單位圓即為問題里的池塘，四只鴨子為出現在單位圓內部隨機位置的四個點，那么這四個點位於同一半圓內的概率即為題目所求的概率值。
- 怎么判斷四個點位於同一個半圓？
四個點可能一下子有點抽象，那么不妨將問題弱化：考慮從四個點變為三個點的情況。怎么判斷三個點位於同一個半圓？
如上圖所示，在圓中有三個點，我們可以畫出三個點組成的三角形的三條邊，不難發現，圓的圓心在三角形的外邊，這時我們一定可以找到圓的一條直徑，將三角形分在圓的某個半圓中。
那如果圓心在三個點組成的三角形內部呢，能找到一條直徑將三個點分在同一個半圓嗎？稍微在頭腦中模擬一下，不難想到，答案是否定的。如下圖所示，無論怎樣選取直徑，都不能將三個點分到同一個半圓中。
於是我們可以得到一個結論，只要圓心位於三角形外，就一定能將三角形分在同一個半圓中。
- 那么推廣到四個點的情況，結論會是怎樣的？
首先，如果圓中隨機分布的四個點共半圓，那么這四個點中的任意三個也一定共半圓，這個很容易理解。那么結論反過來成立嗎？也就是說，圓中有四個點，任意三個點都共半圓，那么這四個點一定共半圓嗎？
我們可以使用反證法。
假設存在圓中有四個點，其中任意三個點都共半圓，但是四個點不共半圓的情況。不失一般性，我們假設四個點為a，b，c，d。如下圖所示，我們畫出a，b，c點和經過這三個點的圓直徑。
三條直徑將圓分為了6個區域，分別為S1~S6。現在我們考慮將第四個點放在什么地方，可以滿足任意三個點共半圓，但是四個點不共半圓的情況。
如果d點放在區域S1，S2或者S3，顯然，abcd共藍色直徑左邊的半圓，矛盾；
如果d點放在區域S6，顯然，abcd共紫色半徑下方的半圓，矛盾；
如果d點放在區域S4，此時abcd不共半圓，但是此時abd點也不共半圓，與假設任意三個點共半圓，矛盾；
如果d點放在區域S5，此時abcd不共半圓，但是此時acd點也不共半圓，與假設任意三個點共半圓，矛盾。
綜上，d點無論放在哪個區域，都不滿足假設的條件，因此假設不成立！
因此我們得到結論，圓內隨機分布的四個點，如果其中任意三個點都共半圓，那么這四個點也一定共半圓。
於是我們得到了判斷四只鴨子位於池塘的同一半圓區域的判別條件。
- 判斷點在三角形外部的方法
判斷一個點在三角形外部還是內部有很多種方法，我常用的是兩種，一種是內角和法：如果點在三角形內部，那么該點與三角形的三個頂點形成的三個角的和為360度，反之，如果點在三角形外部，那么三個角之和一定小於360度。
第二種是同向法。對於三角形Δabc，如果點p在三角形內部，那么對於三角形的每條邊構成的向量——例如向量ab——來說，其另外一個頂點c與該邊的起始點a構成的向量ca，這兩個向量的叉積方向與向量pa和向量ab的叉積方向是一致的，從數值上的表現來看，要么同為正值，要么同為負值。對於向量bc，向量ca的情況，結果也是一致的。而對於在三角形外部的點q來說，至少有一條邊的情況不滿足叉積方向一致。
如上圖中的例子，先看三角形內部的點p。對於邊ab來說，向量ab叉乘向量ca，根據右手定則，叉積的方向指向屏幕里面，向量ab叉乘向量pa，根據右手定則，叉積的方向指向屏幕里面，兩個叉積方向一致；對於邊bc來說，向量bc叉乘向量ab，根據右手定則，叉積的方向指向屏幕里面，向量bc叉乘向量pb，根據右手定則，方向指向屏幕里面，兩個叉積方向一致；對於邊ca來說，向量ca叉乘向量bc，根據右手定則，方向指向屏幕里面，ca叉乘pc，方向同樣是指向屏幕里面。因此可以判斷點p在三角形內部。
而對於q點，我們可以直接看邊ca，向量ca與向量bc的叉乘方向指向屏幕里面，但是向量ca與向量qc叉乘，根據右手定則，叉積的方向指向屏幕外面，兩者方向不一致，可以判斷q點在三角形外部。
我們可以從數值上來驗證這個結論。對於二維向量a(ax,ay)和b(bx,by)，其叉積的數值為ax*by-bx*ay。對於筆者所畫的上圖，如果以左下角為坐標原點，正右方向為x軸方向，正上方向為y軸方向，單個像素為單位長度的話，各點的坐標如下：
a(75, 145), b(266, 76), c(241, 302), p(204, 164), q(95, 252)
先看邊ab，
ab X ca = (191, -69) X (-166, -157) = 191 * (-157) - (-69) * (-166) = -41441
ab X pa = (191, -69) X (-129, -19) = 191 * (-19) - (-129) * (-69) = -12530
ab X qa = (191, -69) X (-20, -107) = 191 * (-107) - (-20) * (-69) = -21887
三個叉乘的方向一致，再看邊bc
bc X ab = (-25, 226) X (191, -69) = -25 * (-69) - 226 * 191 = -44891
bc X pb = (-25, 226) X (62, -88) = -25 * (-88) - 226 * 62 = -11812
bc X qb = (-25, 226) X (171, -176) = -25 * (-176) - 226 * 171 = -34246
三個叉乘方向一致，最后看邊ca,
ca X bc = (-166, -157) X (-25, 226) = -166 * 226 - (-25) * (-157) = -41441
ca X pc = (-166, -157) X (37, 138) = -166 * 138 - 37 * (-157) = -17099
ca X qc = (-166, -157) X (146, 50) = -166 * 50 - 146 * (-157) = 14622
向量ca與向量bc的叉積方向和向量ca與向量pc的叉積方向一致，但是和向量ca與向量qc的叉積方向不一致。所以得出結論，點p在三角形abc內部，而點q在外部。
詳細的論述請參考鏈接http://www.blackpawn.com/texts/pointinpoly/default.html
綜上，我們已經得到了求解題目的所有步驟。
- 代碼
使用同向法判斷點在三角形外部代碼
1 bool beyondTriangle(Point2f& a, Point2f& b, Point2f& c, Point2f& p) 2 { 3 Point2f ab = Point2f(b.x - a.x, b.y - a.y); 4 Point2f ca = Point2f(a.x - c.x, a.y - c.y); 5 Point2f pa = Point2f(a.x - p.x, a.y - p.y); 6 double cross_ca_ab = ca.x * ab.y - ab.x * ca.y; 7 double cross_pa_ab = pa.x * ab.y - ab.x * pa.y; 8 double orientation = cross_ca_ab * cross_pa_ab; 9 if (orientation < 0) return true; 10 Point2f bc = Point2f(c.x - b.x, c.y - b.y); 11 Point2f pb = Point2f(b.x - p.x, b.y - p.y); 12 double cross_ab_bc = ab.x * bc.y - bc.x * ab.y; 13 double cross_pb_bc = pb.x * bc.y - bc.x * pb.y; 14 orientation = cross_ab_bc * cross_pb_bc; 15 if (orientation < 0) return true; 16 Point2f pc = Point2f(c.x - p.x, c.y - p.y); 17 double cross_bc_ca = bc.x * ca.y - ca.x * bc.y; 18 double cross_pc_ca = pc.x * ca.y - ca.x * pc.y; 19 orientation = cross_bc_ca * cross_pc_ca; 20 if (orientation < 0) return true; 21 return false; 22 }
求解的主要代碼如下：
double calcFourDuckInCirclePoolPolarCoordinates() { int PIECE = 36000; int REGION = 100000; int ITER = 1000000; double PI = acos(-1); uniform_int_distribution<unsigned> theta_d(0, PIECE); uniform_int_distribution<unsigned> r_d(0, REGION); default_random_engine theta_e((unsigned int)time(NULL)); default_random_engine r_e(0); vector<Point2f> pos; pos.resize(4); int candies = 0; Point2f center(0.f, 0.f); for (int i = 0; i < ITER; ++i) { for (int duck = 0; duck < 4; ++duck) { auto theta_ = theta_d(theta_e); auto r_ = r_d(r_e); double theta = ((double)theta_ / 100.0 - 180) / 180 * PI; double r = (double)r_*1.0 / REGION; double x = r * cos(theta); double y = r * sin(theta); pos[duck] = Point2f(x, y); } if (!beyondTriangle(pos[0], pos[1], pos[2], center)) continue; if (!beyondTriangle(pos[0], pos[1], pos[3], center)) continue; if (!beyondTriangle(pos[0], pos[2], pos[3], center)) continue; if (!beyondTriangle(pos[1], pos[2], pos[3], center)) continue; ++candies; } double probability = candies*1.0 / ITER; cout << "probability: " << fixed << setprecision(5) << candies*1.0 / ITER << endl; return probability; }
- 結果
多次的采樣實驗結果顯示，通過統計模擬求解四只鴨子位於池塘同一個半圓的概率的答案為0.5！
