數學と計算幾何
hansonyu123
2015 年 11 月 4 日（水曜日）
1 有關數論
1.1 最大公因數和歐幾里德算法
如果d同時整除a, b，就稱d是a, b的公因數 最大的公因數叫作最大公因數
聽說最大公因數叫作gcd(a, b),性質條列如下 1. gcd(a, b)=gcd(b, a) =gcd(a,−b)
2. gcd(a,0)=|a|
3. gcd(a, b)=gcd(a, b−a)
4. (Bezout’s Theorem) 對於所有整數a, b，存在兩個整數x, y 使得ax+by =gcd(a, b) 5. 若gcd(n, b) = 1且n|ab，則n|a
6. gcd(a, b)=gcd(a, b mod a)(3.的直接推論)
利用第2條性質和第6條性質可以輕鬆地算出最大公因數 Algorithm 1: Euclid Algorithm 1 f u n c t i o n gcd ( a , b )
2 i f b==0
3 r e t u r n a
4 else
5 r e t u r n gcd ( b , a mod b )
6 end i f
7 end f u n c t i o n
I
複雜度是O(log(min(a, b)))
注意：上面的虛擬碼是針對非負整數的 a, b所寫的。如果你需要計算負整數的最大公因 數，請務必取絕對值。請務必取絕對值。這實在是太重要了，因為筆者曾被這個bug卡了 很久(orz
第 4 條性質告訴我們存在那樣的 x, y，事實上我們可以將歐幾里德算法擴展一下得到 答案。所以他就叫作擴展歐幾里德囉…
基本的想法是得到(b, a mod b)的答案之後設法算出 (a, b)的答案。實際上去做就是，假 設我們已經有bx′ + (a mod b)y′ =d，那麼就會有ay′ +b(x′−(a−(a mod b))y′/b) = d，
化簡一下就是ay′ +b(x′− ⌊ab⌋y′) =d，所以就可以輕鬆寫出擴展歐幾里德了 Algorithm 2: Extended Euclid Algorithm
1 f u n c t i o n gcd ( a , b )
2 i f b==0
3 r e t u r n ( a , 1 , 0 )
4 else
5 ( d , x , y ) = gcd ( b , a mod b ) 6 r e t u r n ( d , y , x−(a / b ) y )
7 end i f
8 end f u n c t i o n
複雜度仍然是O(log(min(a, b)))
同樣地，這演算法也是為非負整數而生。丟負整數進去之前請記得取絕對值。
另外，因為同時回傳三個值有點麻煩，所以迴圈版本也滿常見的。
1.2 同餘
m|a−b 若且唯若a≡b mod m
同餘有一些很基本的性質，基本到我不想寫出來(誤) 列出幾個可能會用到又不太基本的性質吧
1. (費馬小定理)如果p是質數，那麼ap ≡a mod p
2. (尤拉定理)如果a, n互質，那麼 aϕ(n) ≡1 mod n，這裡 ϕ(n)是小於等於 n且和n 互 質的數的個數。
3. (尤拉判別法) 如果 p 是質數且d 不是 p 的倍數。那麼存在 x2 ≡ d mod p若且唯若
dp−21 ≡1 mod p
值得一提的是，這三個性質配著後面的快速冪一起吃都格外地美味。
同餘裡面有加法，有減法，有乘法，可惜的就是除法不是那麼的完全。如果想要「造」出
II
一個c−1，我們會希望他滿足c−1c≡ 1 mod m。可惜的是根據 Bezout’s Theorem，這要
求gcd(c, m) = 1。所以如果c, m互質，這樣的c−1 一定存在，稱之為模逆元。
模逆元可以直接由擴展歐幾里德求到，複雜度是O(log(min(c, m))).不過如果 m是質數，
由費馬小定理知道cp−2 其實就是c的模逆元。配著快速冪的話複雜度是O(logp)，雖然複 雜度差不多，不過快速冪更好寫(而且通常是需要寫)，所以通常質數的模逆元都是用快速 冪求的。
1.3 快速冪
快速冪剛剛出現了好幾遍，不介紹快速冪好像對不太起自己。
快速冪的目的是快速地求出冪次 (廢話)。想法就是a, a2, a4, . . . 可以透過自己乘自己很快 地求出來，然後對於任何一個 an，都可以在 a, a2, a4, . . . 選不超過 O(log(n))個數字乘起 來得到。實際寫起來大概長得像這樣
Algorithm 3: Fast Power 1 f u n c t i o n power ( a , n )
2 ans = 1
3 while n > 0
4 i f n mod 2 = 1
5 ans = ans * a
6 a = a * a
7 n = n / 2
8 end i f
9 end while
10 r e t u r n ans 11 end f u n c t i o n
複雜度是O(logn).
通常快速冪都是要來求模某個數字的冪次，如果是這樣的話記得在每個地方模一個m 就 好了。另一種想法是an= ((a2)n2)×an mod2，寫起來大同小異。
1.4 質數
建質數表，判斷質數之類的東西已經考到快爛了…怒直接講最常用的XD
要找出2∼N 之間的所有質數，只要使用改進過的篩法就好了。第一個觀察是：只需要把 質數拿出來，將他的倍數們篩掉就好，不需要將合數拿出來篩。第二個觀察是：如果我現 在要把i的倍數全都篩掉，只需要從 i×i開始就好，因為比他小的i的倍數一定會有其它 更小的質因數，所以老早就被篩掉了。
III
實作起來長得像這樣
Algorithm 4: Sieve of Eratosthenes 1 f u n c t i o n e r a t o s t h e n e s ( n )
2 f o r i = 2 t o n
3 prime [ i ] = 1 ( prime [ i ]代 表 i 是 否 是 質 數. 這 行 只 是 初 始 化 )
4 end f o r
5 f o r i =2 , i * i <=n , i ++
6 i f prime [ i ]
7 f o r j = i * i , j <n , j = j + i
8 prime [ j ] = 0
9 end f o r
10 end i f
11 end f o r
12 end f u n c t i o n
複雜度是O(nlog logn)
至於要判斷一個數N 是否是質數，你可以從 2到√
N 判斷是否整除N。 複雜度是O(√
N)
也可以採用Miller-Rabin演算法。其精髓是費馬小定理以及一個質數才有的性質：x2 ≡1 mod p若且唯若x ≡ ±1 mod p.隨便給一個a，要看看他是否通過了這兩個性質的考驗。
如果沒通過，那麼他就是合數。如果通過了，代表他有可能是質數。所以這演算法不保證 正確性，不過隨機用很多種a 測試的話其正確性是頗高的 (99.9趴以上)。令外，據說取 a= 2,7,61，如果三個都通過考驗了，那麼他就一定是質數(在n在int範圍的時候，夠用 了)。
Algorithm 5: Miller-Rabin 1 f u n c t i o n i s p r i m e ( n , a )
2 i f n mod 2 == 0
3 r e t u r n n == 2
4 end i f
5 u = n − 1
6 t = 0
7 while u mod 2 == 1
8 u = u / 2
9 t ++
10 end while
IV
11 a = power ( a , u ) (這 裡 的power是 模n的 快 速 冪 , 複 雜 度O( l o g n ) ) 12 i f a == 1 o r a == n−1
13 r e t u r n t r u e
14 end i f
15 f o r i = 0 t o t−1
16 a = a * a mod n
17 i f a == 1
18 r e t u r n f a l s e
19 else i f a == n−1
20 r e t u r n t r u e
21 end i f
22 end f o r
23 r e t u r n f a l s e 24 end f u n c t i o n
整體的複雜度是O(logn)，相當快速。根據測試，萬惡的ZJ判斷質數可以0.6s過(200000 筆測資)。
質因數分解一個很大的數可以直接從2到√
n找因數，複雜度O(√ n)
也可以採用Pollard’s ρ Algorithm,過程過於複雜，有興趣的請自行搜尋。(演算法筆記有) 而同時將 2至n的數質因數分解只要在篩法時同時實行 DP就好。複雜度不變(但是將一 個數的質因數分解式寫出需要額外代價，複雜度會變為質因數個數，全部展開的話會到 O(nlogn))
1.5 稍微冷僻的東東
這裡再舉兩個和數論有關的東西，然而他們比較不常出現。不過不常出現不代表不重要就 是了。
第一個是盧卡斯定理。這個定理是說，對於所有的質數p和非負整數m, m0, n, n0，都有 Cpn+npm+m00 ≡CnmCnm00 mod p
這裡約定若m < n,則Cnm = 0.如此一來就可以O(p2)的預處理，O(log min(m, n))的時間 內算出Cnm 除以p的餘數。
Algorithm 6: Lucas’ Theorem 1 ( 利 用 巴 斯 卡 恆 等 式 預 處 理C [ 0 ] [ 0 ] 到C [ p−1 ] [ p−1]) 2 f u n c t i o n comb (m, n )
3 i f n == 0
V
4 r e t u r n 1
5 else i f m == 0
6 r e t u r n 0
7 ens i f
8 r e t u r n comb (m/ p , n / p ) * C [m mod p ] [ n mod p ] mod p 9 end f u n c t i o n
第二個是少見複雜度是O(√
n)的演算法。
一個事實是給定一個 n，⌊nx⌋的取值只會有大約2√
n 種。那要怎麼枚舉出這些值呢？最笨 的方法是x從1到n 直接掃一遍。不過這樣太費工了，事實上有更快的方法：
1 f u n c t i o n f ( n )
2 a = 1
3 ( do something )
4 while a < n
5 a = n / ( n / ( a + 1 ) )
6 ( do something )
7 end while
8 end f u n c t i o n
如此一來就會將所有⌊nx⌋枚舉過一遍了。複雜度O(√ n)
1.6 習題
1. (中國剩餘定理) 請設計一套演算法，能計算滿足x ≡ai mod bi(i= 0,1, . . . , n−1)的
x的通式。
2. 給你兩個數d, n.請在O(√
n)的時間內判斷是否存在 x使得x2 ≡d mod n 3. 請試著實作1.4中的所有演算法。
4. 令 f(n)代表 n 的相異質因數個數。給你一個n ≤ 107，請求出x+f(⌊nx⌋)的最小值。
多筆測資，共有104 筆測資。
5. <Codeforces 547A Mike and Frog>
有隻熊，他養了個青蛙和花。青蛙和花一開始身長/高h1, h2.每過一天,青蛙的身長會 變為x1h1+y1 mod m,花會變為x2h2+y2 mod m.請問至少過幾天之後，青蛙和花 的身長/高分別會是a1, a2。
(貼心小提醒：有可能永遠到不了喔ww)
6. <ZJ a289>
求a模m 的模逆元
VI
2 有關矩陣
2.1 基礎的矩陣
相信大家都會矩陣的加法以及乘法的定義，這裡就清描淡寫就好了。
1. 如果A, B 都是m×n的矩陣，那麼(A+B)ij =Aij +Bij 2. 如果A是m×l的矩陣,B 是l×n的矩陣，那麼(AB)ij =∑l
k=0AikBkj
3. 矩陣加法有交換律，結合律
4. 矩陣乘法沒有交換律，有結合律。
矩陣加法的演算法…呃O(n2)不好嗎
矩陣乘法的樸素作法是 O(n3)，然而有更快的 O(nlog27) 的演算法，叫作Strassen algo-
rithm.雖然想法只是一個簡單的Divide and Conqure, 但是細節過於複雜，所以請自己搜
尋∼
矩 陣 的 冪 次 只 要 結 合 矩 陣 乘 法 跟 前 面 所 講 的 快 速 冪 就 好 了， 複 雜 度 O(n3logk) 或 O(nlog27logk)
2.2 高斯消去法
接著來講一下高斯消去法。高斯消去法的意思是說用以下的三種操作使得一個矩陣除了左 上至右下的對角線是1之外，其它都是0
1. 第i列和第j 列交換 2. 第j 列同乘一個常數
3. 第i列加上α倍的第j 列(α 是一個常數)
想法是：對於某一行，挑某一列在那行不是 0，然後把這列拿起來將其它列的那行消成
0(使用第 3 種操作).必要的時候可以交換兩列。至於實作的方法因人而異，而且有點繁
雜，這裡就不詳細寫了。複雜度是O(n2m)
然而如果是整數運算的話，直接使用上面的想法是無法將左上至右下的對角線都消成 1 的。這個時候有兩種替代方案。第一種是折衷，允許對角線上的數字可以不是1(也就是說 可以是2,3之類的，但是不能是0)。另一種方法是輾轉相除法。後者的複雜度會多乘一個 來自輾轉相除的logc.
高斯消去法求聯立方程的解是高中數學範圍，就不多提了。
VII
2.3 計算反矩陣
如果A是一個方陣，在他行列式不為0的時候都存在一個反矩陣A−1。高斯消去法也可以 幫助我們找出反矩陣。想法：實作的時候在A 右邊放一個單位矩陣，一起施行高斯消去 法。假設過程中左邊的矩陣變為 X,右邊的矩陣變為Y，可以證明X =Y A永遠都會保持 住。當高斯消去法結束時，X = I，所以 Y =A−1. 複雜度是O(n3)(如果採用輾轉相除的 話，複雜度是O(n3logc))
2.4 計算行列式
高斯消去法也可以幫助我們算行列式。每次將兩列交換時行列式變號，進行第2個操作時 行列式也乘上相同的常數。而高斯消去法終止時，行列式就是對角線元素的乘積。
然而如果要算行列式，只要讓整個矩陣是上三角矩陣就好了，也就是說高斯消去法可以不 用做得那麼澈底，可以稍微加快效率。複雜度同樣是O(n3)
2.5 快速求線性遞迴
給定 k1, k2, . . . , kn,還有初始項 a0, a1, . . . , an−1. 之後的每一項都滿足ai =∑n
j=1kjai−j(i ≥ n).相信DP課有教過怎麼樣求出a0, . . . , aN 的值。不過矩陣可以幫助我們更快地求出單一 個aN 的值。
以下以最常見的費波那契數列為例。a0 =a1 = 1,ai =ai−1+ai−2.一個神奇的發現是 [
ai+1 ai
] [1 1 1 0 ]
= [
ai+2 ai+1 ]
所以 [
a1 a0
] [1 1 1 0
]N
= [
aN+1 aN ]
因此只要利用快速冪算出那個N 次方就好了。複雜度O(logN).
至於將這個算法推廣就當作作業囉lol
2.6 習題
1. 實作高斯消去法和輾轉相除的高斯消去法 2. 實作反矩陣的演算法
3. 實作計算行列式的演算法 (如果你覺得數字太大的話，模個1000000007) 4. 快速求n 階線性遞迴數列的第N 項
VIII
5. <ZJ d639>
有 一 個 數 列 a1 = a2 = a3 = 1,ai = ai−1 +ai−2 + ai−3 ∀i ≥ 4, 請 求 出 an.n ≤ 2147483647,單筆測資。請輸出模10007後的結果。
6. <ZJ a410>
解二元一次聯立方程式
3 有關組合
3.1 組合計數
以下提供幾個比較常用到的東西。大部分都是高中數學範圍喲 1. 加法原理和乘法原理
2. 排容原理
3. m 個東西中取n 個東西排列的方法數Pnm = (mm!−n)!
4. m 個東西中取n 個東西的方法數Cnm = n!(mm!−n)!
比較值得一提的是P 和 C 的大小都誇張到了一個極點，所以通常不是要模一個數字就是 要用浮點數運算(Ex:算機率).與之相對的，排容原理的複雜度大到了一個極點，所以通常 會在測資範圍很小的時候出沒。千萬不要小看排容原理的威力，只要是用到他的暴搜題，
實作難度都滿高的...
3.2 組合恆等式
以下列了幾個，不過其實也沒有很常用。為了以防萬一，知道一下比較好。
1. Cnm =Cnm−−11Cnm−1 2. ∑m−n
i=0 Cnn+i =Cn+1m+1
3. (二項式定理)(1 +x)n =∑n
i=0Cinxi
第二條貌似比較冷僻，不過如果你沒看過的話第一次見到應該會被嚇到。
3.3 遞迴與組合計數
當你要計算一個東西的個數的時候，一個好方法是令ai 代表n=i的答案，想辦法寫出遞 迴式。
事實上你平常在做 DP 的時候也要做類似的事：寫出轉移的式子。不過DP 通常都會有
IX
max跟min之類的，遞迴式就輕鬆多了。
通常，一個題目能被遞迴輕鬆殺掉若且唯若他是線性遞迴，所以配著前面教過的矩陣快速 冪用吧。
3.4 卡特蘭數
卡特蘭數莫名地熱門，真是奇怪…
有一個 n×n 的網格紙，被一個左下到右上的對角線切。有一個小螞蟻從(0,0)出發，要
到(n, n).請問小螞蟻不跨越對角線的走法有幾個？(螞蟻只能向右或向上走)
我們假設這樣的答案叫作Cn.一開始你會發現
Cn =
n−1
∑
i=0
CiCn−1−i
結合C0 = 1 你可以用O(n2)的時間算出Cn.然而這樣還不夠快。事實上可以證明
Cn=Cn2n−Cn2n−1 = 1 n+ 1Cn2n 這樣就可以在O(n)的時間內算出卡特蘭數了。真是可喜可賀。
為什麼卡特蘭數有名字呢？因為他很重要 (廢話)。重要的點在於以下的題目答案都是卡特 蘭數：
1. n 個節點的二元樹個數 2. n 對括號合法匹配的方法數 3. 凸n+ 2邊形三角化的方法數
4. 著名stack練習題Rails(UVa 514)中答案是”Yes”的排列數
3.5 忽略足夠小的數
通常是在算機率的時候會用到。因為通常輸出浮點數的時候要求的精確度有限，所以可以 將絕對不會影響到精確度的部分捨去以降低複雜度。這要求非常高的估計技巧，所以如果 你不會的話…自己調個常數看看會不會既不TLE又不WA 囉…
啊對了我從來都沒有成功辦到這種事…如果有興趣的話拜CBD為師吧(?)
3.6 習題
1. <2009 TOI初選第一題>
從原點開始，每次只能往右，往上或往左走一單位，且不能走重複的邊。請問長度為 n 的路徑總數是多少？n ≤50
X
2. 承上題，請輸出模 109+ 7 後的答案，n≤108 3. <TIOJ 1607>
計算寬度為n的山梭線總數模109+ 7。n≤106,105 筆測資。
4. (承上題)把答案改模997.n≤2147483647,106 筆測資。
5. <2014 TOI二模pC>
有 n 個人排成一列，編號從 1到n，接下來死神會依照 1, 2,. . . , n 的順序去找每個 人，一個人被死神找到的時候有 12 的機率會死掉，並且死神在找完編號n 的人之後 從編號最小的還沒死的人繼續按照編號由小到大找，重複這個過程直到只剩下一個 人為止。求第 m個人存活下來的機率。1 ≤ m ≤n ≤100000.絕對誤差< 10−9 都算 AC(炸
6. <Codeforces 559D Randomizer>
給你一個凸多邊形(頂點座標依逆時鐘方向給)。從凸多邊形的頂點中隨機挑個子集使 得圍出的面積不是零(也就是不退化), 求所挑出的多邊形嚴格地包住的格點數的期望 值。
4 有關計算幾何
4.1 浮點數誤差
有興趣的人可以試著用C++編譯執行這段code
Algorithm 8: BOOM!!!
1 # i n c l u d e < b i t s / s t d c ++. h>
2 u s i n g namespace s t d ; 3
4 i n t main ( ) 5 {
6 f o r ( double i = 0 ; i ! = 1 ; i += 0 . 0 0 0 1 ) i f ( i >1) p u t s ( ”BOOM” ) ; 7 r e t u r n 0 ;
8 }
他會不斷地吐BOOM給你。為什麼呢？不是迴圈執行了 10000次之後 i就會變成1，就 會跳離迴圈了嗎？
要記注，浮點數是有誤差的。當他加了 10000次0.0001之後，並不會變成 1，而會跟1 差一點點。所以當我們要判斷兩個浮點數是否相等的時候，事實上是問說這兩個數夠近 嗎？所謂夠不夠近的「標準」，我們通常叫他 eps(也就是微小誤差的意思).如果兩個數的
XI
差值小於eps，我們就說他們是相等的。
eps的取值同樣來自高超的估計技巧，不會的話就自己亂試常數吧。
Algorithm 9: Correct code using eps 1 # i n c l u d e < b i t s / s t d c ++. h>
2 u s i n g namespace s t d ; 3
4 i n t main ( ) 5 {
6 i n t c n t = 0 ;
7 double eps = 0 . 0 0 0 0 0 1 ;
8 f o r ( double i = 0 ; abs ( i −1.0) > eps ; i += 0 . 0 0 0 1 ) c n t ++;
9 cout << cn t << e n d l ; 10 r e t u r n 0 ;
11 }
現在他會乖乖地吐出10000給你了。
4.2 內積
兩個向量v⃗1(x1, y1), ⃗v2(x2, y2)的內積定義為v⃗1·v⃗2 =x1x2+y1y2.事實上他也是|v⃗1||v⃗2|cosθ(這 裡θ 是兩向量的夾角)
比較常見的用法是判斷兩向量是否垂直，如果垂直的話內積會是0。
4.3 外積
天啊外積實在是太重要啦!!兩個向量v⃗1(x1, y1), ⃗v2(x2, y2)的外積定義為v⃗1×v⃗2 =x1y2−x1y2. 事實上他也是|v⃗1||v⃗2|sinθ(這裡θ 是兩向量的夾角) 利用這件事實和arccos不難算出兩向 量間的夾角。
聰明的你應該發現外積有負交換律了。這是因為θ 的意思是由v⃗1「逆時鐘」轉到v⃗2 的夾 角。所以外積常常用來判斷兩向量之間的夾角是順時鐘還是逆時鐘，在凸包算法中扮演著 極重要的角色。除此之外，因為長度為 a和b 夾角為θ的三角形面積是 12absinθ，跟外積 公式長得好像喔!所以外積也常來計算面積。
注意：其實外積是個向量，但是計算幾何中似乎都把他當作一個純量對待(他的正負號代 表了一切)
XII
4.4 判斷喵喵相交
判斷相交算是基礎又惹人厭的一個東西。如果你同時要判斷線段相交，半線和線段相交，
半線和半線相交，. . . 諸如此類．勸你還是直接寫求交點的函數判斷吧…
想法是(x1, y1)(x2, y2)可以定出一條直線(x1 +t(x2−x1), y1+t(y2 −y1))(t 是參數，這叫 做直線的參數式)。直線代表t可以是任意數，射線代表t≥0(或者t≤1)，半線代表t≤0
或t ≥ 1，線段代表0 ≤t ≤ 1。好處是可以同時解決一堆case而不用重寫函式，而且如
果座標都是整數的話其實可以只用整數運算。
坊間還有流傳一種判斷線段相交的方法，如果你只需要判斷線段相交的話是不錯的選擇。
想法就是：如果要AB 和 CD 相交的話，A, B 要在 CD 的異側，C, D 要在 AB 的異側。
所以我們只需要判斷兩個點是否在直線的兩側。
這時外積就出動了！不難發現A, B 在CD 的異側若且唯若AC, AD的夾角和 BC, BD的 夾角一個是逆時鐘的，一個是順時鐘的。也就是說AC, AD的外積和BC, BD的外積一正 一負。如此就完成了判斷。
4.5 計算面積
剛剛講過外積可以計算面積。事實上，如果P1, . . . , Pn 是一個凸n邊形，那麼隨便取一個 原點O,這個凸n 邊形的面積會是
1 2|
∑n
i=1
OP⃗ i×OP⃗i+1|
可以用有向面積推得這個結果。你也可以把O 取在凸n 邊形內驗證看看。
4.6 習題
1. 實作判斷各種線相交的函數 2. 實作判斷線段相交的演算法
3. <Codeforces 559A Gerald’s Hexagon>
給你一個內角全都是 120◦ 的六邊形還有他的六個邊長 (都是正整數)，請算出這個六 邊形可以被幾個邊長為 1的正三角形填滿。
(註：想想不用面積怎麼做？用面積怎麼做？)
4. <HOJ404 現充爆炸裝置>
平面上有一些紅點和一些藍點。任三點不共線。兩個紅點連線構成紅線段，兩個藍點 連線構成藍線段。請算出紅線段與藍線段相交的次數。(也就是說如果有兩個紅線段跟 一個藍線段共點的話，還是算兩次)請嘗試在O(n2logn)的時間內完成
註：這題極度血腥請小心 XD
XIII
