\( \newcommand{\ord}[1]{\mathcal{O}\left(#1\right)} \newcommand{\abs}[1]{\lvert #1 \rvert} \newcommand{\floor}[1]{\lfloor #1 \rfloor} \newcommand{\ceil}[1]{\lceil #1 \rceil} \newcommand{\opord}{\operatorname{\mathcal{O}}} \newcommand{\argmax}{\operatorname{arg\,max}} \newcommand{\str}[1]{\texttt{"#1"}} \)

2016年12月7日 星期三

[ 0-1 Fractional Programming Dinkelbach algorithm ] 0-1分數規劃

【定義】

給你兩個數組,$benefit[i]$表示選$i$的利潤、$cost[i]$表示選$i$的花費。
如果選$i$,定義$x[i]=1$否則$x[i]=0$,通常題目會給你選擇的限制條件,像是必須是一顆生成樹之類的,這裡就把它忽略掉
我們的目標是要求$\frac{\sum{benefit[i] \times x[i]}}{\sum{cost[i] \times x[i]}}$的最大值。

【分析】

先不管最大化問題,令$\frac{\sum{benefit[i] \times x[i]}}{\sum{cost[i] \times x[i]}}=L$
可以把式子改寫成$(\sum{benefit[i] \times x[i]})-L \times (\sum{cost[i] \times x[i]})=0$
分離參數得到$\sum{(benefit[i]-L \times  cost[i]) \times x[i]}=0$
只要知道$L$,$(benefit[i]-L \times  cost[i])$是直接可以求出來的,先記他為$d(i,L)$吧
那可以設$f(L)=\sum{d(i,L) \times x[i]}$

我們的目標變成是要最大化$L$
那若存在一種$x$的選法使得$f(L)>0$能夠證明什麼呢?
$f(L)>0 \to \frac{\sum{benefit[i] \times x[i]}}{\sum{cost[i] \times x[i]}}>L$那表示存在另一個$L$會比現在的$L$還要好
那 $f(L)<0$?很明顯這種情況完全沒有意義,因為找不到這種$L$

最佳的$L$會讓所有$x$的選法都$f(L)\leq0$,只要找到這樣的$L$就說明他是最佳解了
也可以用類似的想法求最小化$L$,這裡就不贅述了。

【解法】

根據剛剛找出來的性質$f(L)$是單調性的,我們可以二分搜$L$,然後驗證是否存在一組解使得$f(L)>0$,這很簡單就不附code了。
另一個有效率的算法是Dinkelbach,他跟二分搜不一樣的地方是他要去求解,在驗證解跟求解的複雜度一樣的時候Dinkelbach是會比較快的,但有時候求解的複雜度會比驗證解的複雜度還要高,因此要看情況使用
以下附上虛擬碼:

常見的題型有:最優比率環、最優比率生成樹、最大密度子圖(有flow解)等

2016年12月6日 星期二

[ dominator tree ] 有向圖的支配樹

dominator tree,中文名稱叫支配樹
給一張有向圖$G$,我們從其中一個點$r$出發到$y$的所有路徑中,一定會經過點$x$,就稱$x$是$y$的必經點;我們把距離$y$最近的必經點稱為最近必經點,記為$idom(y)$。
支配樹上的有向邊$(u,v)$滿足$idom(v)=u$,那對於一個點$y$,$y$的所有必經點就會是支配樹上$r$到$y$路徑上的所有點。
這個東西可以求有向圖的割點、橋(在每個邊加入虛擬點)等等。
注意code裡的$idom$跟我定義的不一樣
我是看這份講義的偽代碼寫出來的:

2016年11月3日 星期四

[ Pi algorithms - John Machin's formula ] 梅欽公式計算圓周率

我們知道
$$\frac{\pi}{4}= \arctan 1$$
而$\arctan(x)$的泰勒展開式長這樣
$$\arctan(x)=x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^9}{9}- \cdots \; (x \leq 1)$$
故可以求出莱布尼茨公式如下:
$$\frac{\pi}{4}= 1 \,-\, \frac{1}{3} \,+\, \frac{1}{5} \,-\, \frac{1}{7} \,+\, \frac{1}{9} \,-\, \cdots$$
但是這用來計算$\pi$到小數點後指定位數是非常困難的,因此有了以下的梅欽公式:
$$\frac{\pi}{4} = 4 \arctan \frac{1}{5} - \arctan \frac{1}{239}$$
再用泰勒展開,可以得到:
$$\pi=(\frac{16}{5} - \frac{16} {3*5^3} + \frac{16}{5 * 5^5} - \frac{16} {7 * 5^7} + \cdots) \\ - (\frac{4}{239} - \frac{4}{3 * 239^3} + \frac{4}{5 * 239^5} - \frac{4}{7 * 239^7} + \cdots)$$
可以將這個公式整理為:
$$\pi=\frac{\frac{16}{5} - \frac{4}{239}}{1} -\frac{\frac{16}{5^3} - \frac{4}{239^3}}{3} + \frac{\frac{16}{5^5} - \frac{4}{239^5}}{5} - \cdots$$
如果我們要計算圓周率至10的負$L$次方,由於$\frac{16}{5^{2*n-1}} - \frac{4}{239^{2*n-1}}$中$\frac{16}{5^{2*n-1}}$比$\frac{4}{239^{2*n-1}}$來的大,具有決定性,所以表示至少必須計算至第n項:
$$\frac{16}{(2*n-1)*5^{2*n-1}}=10^{-L}$$
將上面的等式取log並經過化簡,我們可以求得(誤差什麼得先不管它):
$$n =\frac{L}{2log_{10} 5} = \frac{L}{1.39794}$$
這樣就可以計算$\pi$到小數點後第$L$位了

這裡有其他可以快速算出圓周率的梅欽類公式
http://jeux-et-mathematiques.davalan.org/calc/formulespi.pdf

2016年10月28日 星期五

[ Earley parser ] Earley語法解析器

有一天蚯蚓跟我說CYK算法很慢,推薦了這個東西,所以我就把他寫出來了。
他主要是先建構出類似於自動機的圖,有很多個狀態在跳來跳去,在我的實作裡有將其轉成語法樹的方法。
這裡需要自己實作的有lexer,還有最後的語法一定要有一個gamma rule,gamma rule就是有一個$gamma \Longrightarrow S$,$S$是初始符號。
首先是語法規則的部分:

這是分析器的本體,lexer的部分要自己去完成,lexer的結果會是parse的第二個參數,計算出來一些資料會存在table裡面
以下為模板:

最後是構造語法樹的部分,因為div是一顆左兄弟右兒子的樹,所以要將她轉換成正常的語法樹,還要去記錄曖昧文法的部分:

2016年10月11日 星期二

[ Cocke-Younger-Kasami algorithm ] CYK算法

它是一個用來判定任意給定的字符串是否屬於某一個上下文無關文法的算法,還可以順便構造出二元語法樹和判斷是否有曖昧文法(Ambiguous Grammar)。對於一個任意給定的上下文無關文法,都可以使用CYK算法來計算上述問題,所以它是一個幾乎萬能的算法,但首先要將該文法轉換成喬姆斯基範式(CNF)。
這個算法主要是利用動態規劃的思想,給一個有$n$個字符的字符串s,和$R$條語法規則,則他的計算複雜度為$\ord{n^3R}$,空間複雜度為$\ord{n^2R}$。類似於optimal binary search tree的算法,對於每個s[l,r]我們去計算s[l,k]和s[k+1,r]能符合的所有規則。

我們用2016 NCPC的題目當例子
這裡給一個簡單的範例,給你一個語法規則如下 :
$
A \Longrightarrow AAAAAAA \qquad 20 \\
A \Longrightarrow AA \qquad 15 \\
A \Longrightarrow a \qquad 5
$
如果看不懂的,他大約的意思如下:
一個合法字串$A$可以是七個合法字串$A$所組成,也可以是兩個合法字串$A$所組成,也可以是小寫字母$a$所組成。
每個規則都給他一個權重,表示經過這個規則的轉換會有這些花費,我們會給一個目標字串,要找出把整個字串轉換成$A$的最小花費,例如給定字串$aaaaaaaa$,他的最小花費就是75。
這裡舉另一個例子:
$
A \Longrightarrow BA \qquad 10 \\
A \Longrightarrow bcd \qquad 5 \\
B \Longrightarrow c \qquad 4
$
給定字串$ccbcd$,他的最小花費就是33。
當然也會發生無法轉換的情況,這個時候叫要輸出$NIR$ 。

首先必須要把規則轉成喬姆斯基範式,這裡我以數字當作規則的名稱,當y=-1時表示這個cnf是s->x的規則,否則是s->xy的規則。
接著就可以跑我們的演算法啦:
當然它也可以用來判斷是不是有曖昧語法(ambiguous),只要做一些小修改就行了

2016年10月3日 星期一

[ cantor expansion ] 康托展开

給一個$0 \sim n-1$共$n$個數字的全排列,求他是排名第幾的全排列;給一個名次,求全排列。這東西可以有效的減少hash_table的大小。
以下附上code:

使用前記得要呼叫init();

2016年8月5日 星期五

[ tree centroid ] 樹重心

一棵無向樹\(T=(V,E)\),定義:
\(w_v(u)\)為點\(u\)的權重,\(u \in V\)
\(w_e(e)\)為邊\(e\)的權重,\(e \in E\)
\(d(u,v)\)為\(u\)到\(v\)路徑的權重和,\(u,v \in V\)

重心的定義是:
以這個點為根,那麼所有的子樹(不算整個樹自身)的大小都不超過整個樹大小的一半
即為最大的子樹最小的點

廣義的樹的重心:
無向樹\(T=(V,E)\)滿足
 \(w_v(u)≧0, \; \forall u \in V \)
 \(w_e(e)>0, \; \forall e \in E \)

\(c(u)=\sum_{v \in V}d(u,v)*w_v(v)\)
則樹重心 \(tree \; centroid=u \; | \; min(\{c(u):u \in V\})\)


 \(w_v(u)=1, \; \forall u \in V \)
 \(w_e(e)=1, \; \forall e \in E \)
的情況下,就是一般的樹重心

性質:
  1. 把兩個樹通過一條邊相連得到一個新的樹,那麼新的樹的重心在連接原來兩個樹的重心的路徑上。
  2. 把一個樹添加或刪除一個葉子,那麼它的重心最多只移動一條邊的距離。
用途:
樹分治、動態求樹重心等

可以利用DFS找出最大的子樹最小的點,即為樹重心
樹重心求法(用vector存無向樹):

[ bipartite graph multiple matching ] 二分圖多重匹配增廣路算法

這種問題的題目通常是這樣:

給一張圖G有n1個點和n2個點,n1個點之間沒有邊,n2個點之間也沒有邊,但是n1和n2個點之間有m條邊(簡單的來說就是n1個點和n2個點的二分圖啦),沒有重邊;其中n2個點每個點\(u\)都有一個可接受匹配數 \(c_u\)。
n1的點只能跟一個點匹配,但n2的點在不超過可接受匹配數的情況下,可以跟多個點匹配,求這張圖的最大匹配

通常可以利用拆點的方法或是flow來做,但是這樣空間跟效率都不是很好,這裡提供一個有效率的方法:

複雜度分析:
設n2個點每個點\(u\)都有一個值\(E_u\)表示和\(u\)相鄰的邊數,則總複雜度為
$$ \ord{n1*(\sum{c_u*E_u}+n1+n2)} $$

2016年5月19日 星期四

[ Kuhn-Munkres Algorithm ] 二分圖最大權完美匹配KM算法

關於此算法的介紹及教學請看這份投影片,這篇文章注重於討論$\ord{N^4}$的slack優化及$\ord{N^4}$到$\ord{N^3}$的過程。

在網路上曾經看到這段話:

"實際上KM算法的複雜度是可以做到$\ord{N^3}$的。我們給每個y頂點一個“鬆弛量”函數slack,
每次開始找增廣路時初始化為無窮大。在尋找增廣路的過程中,檢查邊(i,j)時,如果它不在相等子圖中,則讓slack[j]變成原值與lx[i]+ly[j]-w[i,j]的較小值。這樣,在修改頂標時,取所有不在交錯樹中的Y頂點的slack值中的最小值作為d值即可。但還要注意一點:修 改頂標後,要把所有不在交錯樹中的Y頂點的slack值都減去d。"

這段話其實是錯的,請看底下的code:

$\ord{N^4}$常數優化版本
我們仔細看一下迴圈的層數。
第一層是一個for迴圈,執行N次;第二層有一個無線迴圈,但她最多執行N次;裡面有一個DFS,最多執行$\ord{N^2}$次,還有一些迴圈但是不影響複雜度所以不討論。
總複雜度是$\ord{N}*\ord{N}*\ord{N^2}$為$\ord{N^4}$,因此網路上大多數$\ord{N^3}$的code其實是錯的。

接下來看一下我修改後真正的$\ord{N^3}$code:
可以看到我在dfs裡面加了一個參數adjust,如果他是true,就跟原本的dfs沒兩樣,可以在找到增廣路後順便將增廣路擴充;如果他是false,整個dfs就變成只能判斷有沒有增廣路了。
主函數dfs的部分移到while的外面,而修改lx,ly和slack_y後的地方增加了一個迴圈來判斷這次的修改有沒有產生可行的增廣路,對於新增加的「等邊」,如果他有機會是增廣路的話,會對此「等邊」連像的點y進行dfs「判斷」有沒有增廣路(把adjust設成0)。如果有增廣路,就會跳出while,在進行一次dfs來擴充增廣路。
總複雜度是
$\ord{N}[第一層迴圈] * ( \ord{N^2}[dfs的時間] + \ord{N}[while的次數]*\ord{N} ) = \ord{N^3}$
這才是真正的$\ord{N^3}$算法,前者只是常數優化而已

最後是$\ord{N^3}$算法常數較小的版本:

2016年4月21日 星期四

[ Mersenne Twister ] 梅森旋轉算法

梅森旋轉算法是一種可以快速產生高級偽隨機數的算法,修正了很多以前的隨機速算法的缺陷(Ex:線性同於算法)。

這是專門用在蒙地卡羅測試用的,不要拿它來做持久化隨機二叉數的亂數產生器,會太慢

最為廣泛使用Mersenne Twister的一種變體是MT19937,可以產生32位整數序列,從C++11開始,C++也可以使用這種算法。在Boost C++,Glib和NAG數值庫中,作為插件提供。但是一般在競賽時可能會有不支援的情況發生,所以將它做成模板當作參考。

MT19937_64則是可以產生64位整數序列

以下提供模板:

2016年4月6日 星期三

[ chinese remainder theorem ] 中國剩餘定理

我是看維基百科實現的,所以就直接貼上模板
注意:
  • crt函數的m[i]是模數
  • crt函數的a[i]是原本數模m[i]後的值
  • 必須要保證m兩兩互質
  • 容易溢位
模板:

2016年4月2日 星期六

[ Number Theoretic Transform ] 快速數論變換)(NTT)

這也不知道怎麼說明,應該是我數學太爛的關係,直接看連結吧
有以下幾點要注意:
  1. P要是一個q*2^k+1的質數,G是P的原根
  2. N(數列長度)要是2的冪次
  3. (P-1)%N=0
  4. 做了逆轉換後的結果為原本數列mod P的結果,所以P要夠大
  5. 如果要求P是任意數的話,可以用中國剩餘定理處理
這裡有一些質數和它的原根(按質數順序由大到小排):
  • 2615053605667*(2^18)+1,3
  • 15*(2^27)+1,31
  • 479*(2^21)+1,3
  • 7*17*(2^23)+1,3
  • 3*3*211*(2^19)+1,5
  • 25*(2^22)+1,3
因為可能會有A*B%M的情況發生,所以需要用
long long 乘 long long 模 long long不溢位算法
只需要將code裡面所有的A*B%M的部分改掉即可

以下提供模板(bit_reverse參考自Morris):

[ fast fourier transform - Cooley Tukey ] 快速傅立葉轉換(FFT)

因為不知道怎麼說明所以提供一些網址:
注意:
  • N(數列長度)要是2的冪次
以下提供模板(bit_reverse參考自Morris):

2016年2月25日 星期四

[ dynamic kd tree ] 動態kd樹模板

參考資料:
kd樹的資料在網路上其實滿多的,但很多code效率不高或是功能不全,所以這邊稍微講解一下kd樹基本的一些操作
  1. 構造
    現在有一個point序列S[0~n-1],每個point都包含一個序列d[0~kd-1]表示維度,S[i].d[j]表示第i個點的第j維。
    節點的資料型別定義為node,裡面有l,r兩個指標分別指向左右兩顆子樹,還有一個point pid表示所存的point。
    構造的過程是一個遞迴結構,大致代碼如下:
    當前這個節點的k稱為分裂維度
  2. 查找
    kd樹支援兩種查找法:
    一、給定一個點p和一個數字k,求離p前k近的點有哪些
    二、給定一個範圍,求範圍內有哪些點
    兩種方法跟一般爆搜有點像,但是利用了kd樹可以做到有效的剪枝,有時候可以判斷被切分出來的範圍內有沒有機會存在前k近點或是在範圍內,直接看模板應該就懂了
  3. 刪除
    模板裡沒有附刪除的code所以會寫在這裡。
    首先如果找到要刪除的節點是葉節點直接刪除就好了;如果不是葉節點,假設現在這個點的分裂維度是k,就拿他右子樹第k維最小節點mi去替代他,接著遞迴刪除右子樹的mi;如果沒有右子樹,就拿他左子樹第k維最小節點mi去替代他,然後把左右子樹互換,接著遞迴刪除右子樹的mi。

    找到要刪除的點用查找中的第一種方法就好了,這裡p=要刪除的點,k=1,查找的時候順便維護最近節點其位置與其對p的距離,若最近距離!=0則刪除失敗。

    那對一個子樹找第k維最小節點呢?方法很簡單,也是遞迴定義的:
    首先如果當前節點o的分裂維度剛好是第k維,則若o有右子樹的話答案必在o的右子樹,否則答案為o,如果o的分裂維度!=k,則遞迴搜尋左右子樹,把得到的答案和o自己進行比較,求最小。
    接下來附上只有刪除操作的模板:

插入刪除通常配合替罪羊樹使用,插入很簡單模板有就不說了,接下來講一下複雜度:
設N=size(root)
  • 構造:
    很明顯就是\(\ord{N \; log \; N}\)就不說了
  • 查找:
    已經有人證明了,假設現在的維度有k維,則查找的最差複雜度是\(\ord{N^{1-1/k}}\)
    但是平均狀況下為\(\ord{log \; N}\)
  • 插入:
    複雜度為\(\ord{樹的高度}\),因為是套替罪羊樹,而且重構的時間是\(\ord{N \; log \; N}\),所以單次插入的均攤時間複雜度是\(\ord{log \; N \; * \; log \; N}\)
  • 刪除:
    有三個步驟需要分析,假設現在的維度有k維
    1. findmin:
      最多會遍歷\(\alpha^{(k-1)*(log_{\alpha}N)/k}=N^{1-1/k}\)個節點,所以是\(\ord{N^{1-1/k}}\)
      但實際操作量為\(N^{1-1/k}-n^{1-1/k}\),n是最小值的子樹大小
    2. nearest:
      就是查找操作,所以是\(\ord{N^{1-1/k}}\)
      但因為是找相同點,可以優化code,所以實際操作量為\(N^{1-1/k}-n^{1-1/k}\),n是相同點的子樹大小
    3. 刪除操作本身:
      看code明顯就是重複findmin直到把刪除的點變成葉節點為止,會感覺做了很多操作,但是我們發現把所有操作加起來後會=\(\ord{N^{1-1/k}}\)
      設\(o_1,o_2,...,o_d\)=findmin(root)找到的點,findmin(\(o_1\))找到的點,\(...\),要刪除的葉節點,\(n_1,n_2,...,n_d\)=findmin(root)找到的點的子樹大小,findmin(\(o_1\))找到的點的子樹大小,findmin(\(o_2\))找到的點的子樹大小,\(...\),要刪除的葉節點的子樹大小,\(d\)為樹的高度
      複雜度為:
      $$ N^{1-1/k}-n1^{1-1/k}+n1^{1-1/k}-n2^{1-1/k}+n2^{1-1/k}...-nd^{1-1/k}=\ord{N^{1-1/k}} $$
模板裡的code找的是曼哈頓距離,如果要找歐基里德距離只需要修改point裡的dist跟kd tree裡的heuristic即可,以下提供修改方法:
這樣查找時回傳的就是歐基里德距離的平方

以下附只有插入操作的模板:

以下附支援插入刪除的模板:
模板的使用方法請參考:

日月卦長的解題紀錄 [ IOICAMP2016 ] 動態曼哈頓最短距離


2016年1月31日 星期日

[ Miller-Rabin Strong Probable-prime Base ] 質數測試算法保證long long範圍內正確

這本來是隨機演算法,但是已經有一些人,找出一些特定底數,保證判定結果在某些範圍內正確。

如果是unsigned long long範圍內東西記得要用long long乘long long mod long long 的模板
這裡的範例已經加在code裡了,如果不需要可以把它拿掉(zerojudge a007不拿掉會變慢)

code:

2016年1月27日 星期三

[ Biconnected Component ] 雙連通分量(BCC)

一張無向圖上,不會產生割點的連通分量,稱作「雙連通分量」。一張無向圖的雙連通分量們,通常會互相重疊,重疊的部分都是原圖的割點。
可以利用堆疊+tarjan算法快速地找出那些點是割點和雙連通分量
以下提供模板:

2016年1月3日 星期日

[ Strongly Connected Component - Tarjan & Kosaraju ] 強連通分量的兩種做法(dfs)

一張有向圖,找到所有的 Strongly Connected Component ;亦可進一步收縮所有的 SCC 、收縮所有的環,讓原圖變成 DAG 。

方法1-Tarjan法
跟BCC差不多,直接提供模板:

方法2-Kosaraju
對圖做一次DFS,求他的時間戳,再把圖反向,逆著時間戳進行第二次DFS,所遍歷的點就會是同一個SCC,記得不要走重複的點
以下提供模板:

用法:

[ Bridge & Bridge-connected Component ] 橋和橋連通分量(BCC)

只要從一張無向圖上移除了,就會讓這張圖分離成更多部分,呈現不連通的狀態。
注意!必須要考慮重邊的問題:
  1. 對於有向圖的強連通分量,重邊是沒有影響的,因為強連通只要求任意兩點可以互相連通
  2. 對於無向圖的點雙連通分量,重邊也是沒有影響的,因為點雙連通要求是任意兩點之間至少存在兩條點不重複的路徑,對邊的重複並沒有要求
  3. 對於無向圖的邊雙連通分量,重邊就有影響了,因為邊雙連通要求任意兩點之間至少存在兩條邊不重複的路徑
這裡提供模板:

一張無向圖上,不會產生橋的連通分量,稱作橋連通分量,或是橋雙連通分量
這裡提供模板,code跟橋差不多,只是需要用堆疊去記錄當前節點而已:


[ Articulation Vertex ] 割點

關節點(割點)是讓一張無向圖維持連通,不可或缺的點。只要從一張無向圖上移除了關節點(以及與之相連的邊),就會讓這張圖分離成更多部分,呈現不連通的狀態。
存在一種線性算法找出一張圖所有的割點,一個dfs就能解決了

以下提供模板: