[ round-off error eps ] 浮點數誤差模板

有些計算幾何的過程中會產生捨位誤差，這時候就需要決定一個很小的數，去避開這個捨位誤差問題，我們稱那個很小的數為EPS

• A < B ....... A - B < - EPS
• A > B ....... A - B > EPS
• A == B ....... abs(A - B) <= EPS
• A <= B ....... A - B <= EPS
• A >= B ....... A - B >= - EPS

因為我的計算幾何模板在實作時是用template實作的，所以如果需要處理浮點數誤差時，需要傳入一個新的資料結構，所以寫了這份模板

	#ifndef SUNMOON_DOUBLE
	#define SUNMOON_DOUBLE
	#include <cmath>
	#include <iostream>
	const double EPS = 1e-9;
	struct Double {
	double d;
	Double(double d = 0) : d(d) {}
	Double operator-() const { return -d; }
	Double operator+(const Double &b) const { return d + b.d; }
	Double operator-(const Double &b) const { return d - b.d; }
	Double operator*(const Double &b) const { return d * b.d; }
	Double operator/(const Double &b) const { return d / b.d; }
	Double operator+=(const Double &b) { return d += b.d; }
	Double operator-=(const Double &b) { return d -= b.d; }
	Double operator*=(const Double &b) { return d *= b.d; }
	Double operator/=(const Double &b) { return d /= b.d; }
	bool operator<(const Double &b) const { return d - b.d < -EPS; }
	bool operator>(const Double &b) const { return d - b.d > EPS; }
	bool operator==(const Double &b) const { return fabs(d - b.d) <= EPS; }
	bool operator!=(const Double &b) const { return fabs(d - b.d) > EPS; }
	bool operator<=(const Double &b) const { return d - b.d <= EPS; }
	bool operator>=(const Double &b) const { return d - b.d >= -EPS; }
	friend std::ostream &operator<<(std::ostream &os, const Double &db) {
	return os << db.d;
	}
	friend std::istream &operator>>(std::istream &is, Double &db) {
	return is >> db.d;
	}
	};
	#endif

view raw eps.cpp hosted with ❤ by GitHub

[ long long multiply long long mod long long ] long long 乘 long long 模 long long不溢位算法

這標題有點長，不過也很清楚的顯示了這篇的重點，所以就直接附模板吧
如果a,b<m的話a%=m,b%=m;可以拿掉沒關係

code(比較慢，但適用範圍更大，算法取自維基百科:同餘):

	using ull = unsigned long long;
	ull _mod_mul(ull a,ull b,ull m){
	a%=m,b%=m;
	ull ans=0;
	for(;b;b>>=1){
	if(b&1){
	ans+=a;
	if(ans>=m)ans-=m;
	}
	a<<=1;
	if(a>=m)a-=m;
	}
	if(ans>=m)ans-=m;
	else if(ans<0)ans+=m;
	return ans;
	}

view raw mod_mul1.cpp hosted with ❤ by GitHub

[ Maze generation-Randomized Kruskal's algorithm ] 迷宮生成-隨機Kruskal算法

說得通俗點，就是"隨機拆牆"。
一開始假設所有牆都在，也就是圖的每個點都不連通。如果當前不滿足“任意兩個格子間都有唯一路徑”，那麼就隨機選擇一堵牆，要求牆兩端的格子不連通，即它們對應的兩個頂點在兩個不同的連通分量。把這堵牆拆掉，即在後牆兩端的格子對應的頂點間建立一條邊。重複作下去，直到所有的點都在同一個連通分量裡面。
會發現這就是隨機合併點的Kruskal演算法，需要用到並查集

範例生成圖:
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃
＃　　　＃　　　　　　　＃　　　＃　　　　　　　＃
＃＃＃　＃　＃　＃　＃＃＃　＃＃＃＃＃＃＃　＃＃＃
＃　　　　　＃　＃　　　＃　　　＃　＃　　　　　＃
＃＃＃＃＃＃＃＃＃　＃＃＃　＃　＃　＃＃＃＃＃　＃
＃　　　　　　　＃　＃　　　＃　　　＃　　　　　＃
＃＃＃　＃＃＃　＃　＃　＃　＃＃＃　＃＃＃　＃＃＃
＃　　　＃　　　＃　　　＃　　　＃　　　　　　　＃
＃＃＃　＃　＃＃＃＃＃　＃＃＃　＃＃＃＃＃　＃＃＃
＃　　　＃　　　　　　　＃　＃　　　　　＃　＃　＃
＃＃＃＃＃＃＃　＃　＃＃＃　＃　＃＃＃＃＃＃＃　＃
＃　　　＃　＃　＃　　　＃　　　　　　　＃　＃　＃
＃　＃＃＃　＃＃＃　＃＃＃　＃　＃＃＃＃＃　＃　＃
＃　＃　　　　　＃　＃　　　＃　＃　＃　＃　＃　＃
＃　＃　＃＃＃　＃　＃＃＃　＃＃＃　＃　＃　＃　＃
＃　＃　＃　　　　　＃　　　＃　　　　　　　　　＃
＃　＃　＃＃＃　＃＃＃＃＃＃＃＃＃　＃　＃＃＃＃＃
＃　　　＃　　　＃　＃　＃　　　　　＃　　　＃　＃
＃＃＃＃＃＃＃　＃　＃　＃　＃　＃＃＃　＃＃＃　＃
＃　　　　　＃　　　　　　　＃　　　＃　　　　　＃
＃　＃　＃＃＃＃＃　＃＃＃＃＃＃＃＃＃　＃　＃＃＃
＃　＃　　　　　　　　　　　　　　　＃　＃　　　＃
＃＃＃　＃　＃　＃＃＃　＃　＃＃＃＃＃＃＃　＃＃＃
＃　　　＃　＃　　　＃　＃　　　＃　　　　　　　＃
＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

code(n,m必須是奇數):

	#include<bits/stdc++.h>
	using namespace std;
	struct P{
	P(int q,int w,int e,int r):a(q),b(w),x(e),y(r){}
	int a,b;
	int x,y;
	};
	char s[105][105];
	int id[105][105],cnt;
	int st[10005];
	vector<P>p;
	inline int find(int x){
	return st[x]==x?x:st[x]=find(st[x]);
	}
	inline void print(int n,int m){
	for(int i=0;i<n;++i){
	for(int j=0;j<m;++j){
	if(s[i][j]=='#')printf("＃");
	else printf("　");
	}
	puts("");
	}
	}
	int main(){
	int n=25,m=25;
	for(int i=0;i<n;++i){
	for(int j=0;j<m;++j){
	s[i][j]='#';
	}
	}
	for(int i=1;i<n;i+=2){
	for(int j=1;j<m;j+=2){
	id[i][j]=++cnt;
	s[i][j]=' ';
	}
	}
	for(int i=1;i<n-1;++i){
	for(int j=1;j<m-1;++j){
	if(s[i][j]=='#'){
	if(s[i][j-1]!='#'&&s[i][j+1]!='#'){
	p.push_back(P(id[i][j-1],id[i][j+1],i,j));
	}else if(s[i-1][j]!='#'&&s[i+1][j]!='#'){
	p.push_back(P(id[i-1][j],id[i+1][j],i,j));
	}
	}
	}
	}
	for(int i=1;i<=cnt;++i)st[i]=i;
	srand(time(0));
	random_shuffle(p.begin(),p.end());
	for(int i=0;i<(int)p.size();++i){
	if(find(p[i].a)!=find(p[i].b)){
	st[find(p[i].a)]=find(p[i].b);
	s[p[i].x][p[i].y]=' ';
	}
	}
	print(n,m);
	return 0;
	}

view raw random kruskal.cpp hosted with ❤ by GitHub

[ min-max heap ] 最小-最大堆

最小-最大堆(Min-Max Heaps)是一個完整二元樹。此二元樹是交替的階層方式呈現，分別為最小階層 ( min level ) 和最大階層 ( max level ) ，這裡實作樹根為最小鍵值。

最小-最大堆是一種堆積支援以下幾種操作:

1. 求最大值 - max()
2. 求最小值 - min()
3. 刪除最大值 - pop_max()
4. 刪除最小值 - pop_min()
5. 元素入堆 - push()

詳細演算法可以參考原始論文
Min-Max Heaps and Generalized Priority Queues

如果不懂敘述就看code吧
以下提供模板:

	#ifndef SUNMOON_MIN_MAX_HEAP
	#define SUNMOON_MIN_MAX_HEAP
	#include<vector>
	#include<algorithm>
	template<typename T,typename _Seq=std::vector<T>,typename _Compare=std::less<T> >
	class minmax_heap{
	_Seq s;
	_Compare cmp;
	bool is_min_level(size_t id){
	return !(std::__lg(++id)%2);
	}
	void TrickleDown(size_t id){
	if(is_min_level(id)) TrickleDownMin(id);
	else TrickleDownMax(id);
	}
	size_t get_min_descendant(size_t id){
	size_t res=(++id)*2-1;
	for(size_t i=2;i<=4;i*=2){
	for(size_t j=0;j<i;++j){
	int descendant=id*i+j-1;
	if(descendant>=s.size()) return res;
	if(cmp(s[descendant],s[res])){
	res=descendant;
	}
	}
	}
	return res;
	}
	void TrickleDownMin(size_t id){
	while((id+1)*2<=s.size()){
	size_t m=get_min_descendant(id);
	if(cmp(s[m],s[id])){
	std::swap(s[m],s[id]);
	if(m<=(id+1)*2) return;
	if(cmp(s[(m+1)/2-1],s[m])){
	std::swap(s[(m+1)/2-1],s[m]);
	}
	id=m;
	}else return;
	}
	}
	size_t get_max_descendant(size_t id){
	size_t res=(++id)*2-1;
	for(size_t i=2;i<=4;i*=2){
	for(size_t j=0;j<i;++j){
	int descendant=id*i+j-1;
	if(descendant>=s.size()) return res;
	if(cmp(s[res],s[descendant])){
	res=descendant;
	}
	}
	}
	return res;
	}
	void TrickleDownMax(size_t id){
	while((id+1)*2<=s.size()){
	size_t m=get_max_descendant(id);
	if(cmp(s[id],s[m])){
	std::swap(s[m],s[id]);
	if(m<=(id+1)*2) return;
	if(cmp(s[m],s[(m+1)/2-1])){
	std::swap(s[(m+1)/2-1],s[m]);
	}
	id=m;
	}else return;
	}
	}
	void BubbleUp(size_t id){
	if(!id) return;
	int parent=(id+1)/2-1;
	if(is_min_level(id)){
	if(cmp(s[parent],s[id])){
	std::swap(s[parent],s[id]);
	BubbleUpMax(parent);
	}else BubbleUpMin(id);
	}else{
	if(cmp(s[id],s[parent])){
	std::swap(s[parent],s[id]);
	BubbleUpMin(parent);
	}else BubbleUpMax(id);
	}
	}
	void BubbleUpMin(size_t id){
	while(id>2){
	int grandparent=(id+1)/4-1;
	if(cmp(s[id],s[grandparent])){
	std::swap(s[id],s[grandparent]);
	id=grandparent;
	}else break;
	}
	}
	void BubbleUpMax(size_t id){
	while(id>2){
	int grandparent=(id+1)/4-1;
	if(cmp(s[grandparent],s[id])){
	std::swap(s[id],s[grandparent]);
	id=grandparent;
	}else break;
	}
	}
	size_t find_max_id()const{
	size_t res=0,n=std::min(s.size(),size_t(3));
	while(--n) if(cmp(s[res],s[n])) res=n;
	return res;
	}
	void erase_id(size_t id){
	s[id]=s.back();
	s.pop_back();
	TrickleDown(id);
	}
	public:
	bool empty()const{
	return s.empty();
	}
	int size()const{
	return s.size();
	}
	void push(const T&data){
	s.push_back(data);
	BubbleUp(s.size()-1);
	}
	const T&max()const{
	return s[find_max_id()];
	}
	const T&min()const{
	return s[0];
	}
	void pop_max(){
	erase_id(find_max_id());
	}
	void pop_min(){
	erase_id(0);
	}
	};
	#endif

view raw minmax_heap.cpp hosted with ❤ by GitHub

[ link-cut tree ] 動態樹教學+模板

因為code太多容易顯示不出來，這裡有Gist連結
動態樹，一種支持斷邊、連邊的資料結構（要保證是樹或森林），可以維護路徑、但不能維護子樹。
類似於樹鏈剖分（建議先學完樹鏈剖分，再學LCT）。
只不過樹鏈剖分的輕重邊是根據子樹大小而決定的，這棵樹只能是靜態的。
而LCT中的“輕重邊”是根據最後一次操作決定的，最後一次處理了x這個點，那麼從x到根的路徑就是重鏈（每個父親只與一個兒子的連邊為重邊）。
用splay來維護每一條重鏈，以點的深度為關鍵字，即splay中的每個點左兒子深度都比他小，右兒子的深度都比他大（若只有一個點，那麼這個點單獨是一棵splay）。
如果還不知道splay tree的可以參考這裡

splay tree及其節點的一些函數:

	#include<vector>
	struct splay_tree{
	int ch[2],pa;/*子節點跟父母*/
	bool rev;/*反轉的懶惰標記*/
	splay_tree():pa(0),rev(0){ch[0]=ch[1]=0;}
	};
	std::vector<splay_tree> node;
	/*
	有的時候用vector會TLE，要注意
	這邊以node[0]作為null節點
	*/
	inline bool isroot(int x){/*判斷是否為這棵splay tree的根*/
	return node[node[x].pa].ch[0]!=x&&node[node[x].pa].ch[1]!=x;
	}
	inline void down(int x){/*懶惰標記下推*/
	if(node[x].rev){
	if(node[x].ch[0])node[node[x].ch[0]].rev^=1;
	if(node[x].ch[1])node[node[x].ch[1]].rev^=1;
	std::swap(node[x].ch[0],node[x].ch[1]);
	node[x].rev^=1;
	}
	}
	void push_down(int x){/*將所有祖先的懶惰標記下推*/
	if(!isroot(x))push_down(node[x].pa);
	down(x);
	}
	inline void up(int x){}/*將子節點的資訊向上更新*/
	inline void rotate(int x){/*旋轉，會自行判斷轉的方向*/
	int y=node[x].pa,z=node[y].pa,d=(node[y].ch[1]==x);
	node[x].pa=z;
	if(!isroot(y))node[z].ch[node[z].ch[1]==y]=x;
	node[y].ch[d]=node[x].ch[d^1];
	node[node[y].ch[d]].pa=y;
	node[y].pa=x,node[x].ch[d^1]=y;
	up(y);
	up(x);
	}
	inline void splay(int x){/*將節點x伸展到所在splay tree的根*/
	push_down(x);
	while(!isroot(x)){
	int y=node[x].pa;
	if(!isroot(y)){
	int z=node[y].pa;
	if((node[z].ch[0]==y)^(node[y].ch[0]==x))rotate(x);
	else rotate(y);
	}
	rotate(x);
	}
	}

view raw lct_splay_tree.cpp hosted with ❤ by GitHub

定義：

1. Preferred child：
   如果一個點是他父親結點的所有兒子中最後被訪問的，那麼他是他父親的Preferred child。
   每個點最多有一個Preferred child。如果x結點是最後被訪問的點，那麼他沒有Preferred child。
2. Preferred edge:
   每個點到自己的Preferred child的邊叫Preferred edge（相當於樹鏈剖分中的重邊）
3. Preferred path:
   由Preferred edge組成的鏈叫Preferred path（相當於樹鏈剖分的重鏈）

簡單來說就是一個節點在重鏈上的直接兒子定義為preferred child，連向該兒子的邊定義為preferred edge。

access(x):
我們用操作access(x)來表示訪問節點x為LCT的最基礎操作，使x到根結點的路徑成為重鏈。

如上圖，粗邊為重邊。

首先把x旋到他所在的splay的根部。
因為當前的x點是最後一個被訪問的點，所以他沒有Preferred child，因此要把右兒子與他斷開，右兒子成為新的一條鏈。
x的父親是y，把y旋到他所在splay的根結點，把此時y的右兒子與他斷開，把x連過來即可，相當於合併y和x所在的splay。
然後，再找y的父親，合併y與y的父親所在的splay，一直到根結束。
最終，x到根結點的路徑成為新的Preferred edge，而曾經的Preferred child都被拋棄。

	inline int access(int x){
	int last=0;
	while(x){
	splay(x);
	node[x].ch[1]=last;
	up(x);
	last=x;
	x=node[x].pa;
	}
	return last;/*回傳access後splay tree的根*/
	}

view raw lct_access.cpp hosted with ❤ by GitHub

有了access操作後，其他的操作就更容易實現了

make_root(x):
將x設為其所在的樹的根
首先access(x)，這樣x與根結點處在同一棵splay中。
設a=access(x)，而a是splay tree的根，因為要讓x成為整棵樹的根，所以x的深度要最小，我們只要在a上加一個區間反轉的標記rev就可以完成。
這裡提供兩種方式，最後都會把節點x轉到splay tree根的位置(這很重要)
1.

	inline void make_root(int x){
	access(x),splay(x);
	node[x].rev^=1;
	}

view raw lct_make_root1.cpp hosted with ❤ by GitHub

2.

	inline void make_root(int x){
	node[access(x)].rev^=1;
	splay(x);
	}

view raw lct_make_root2.cpp hosted with ❤ by GitHub

link(x,y):
連接x與y所在的兩棵樹。
先判斷是否在同一棵splay tree裡，如果是的話就不做任何事
首先Makeroot(x)，x就成為了他所在樹的根。
然後直接把x的父親賦值為y即可。

	inline void link(int x,int y){
	make_root(x);
	node[x].pa=y;
	}

view raw lct_link.cpp hosted with ❤ by GitHub

cut(x,y):
斷開x與y所相連的邊。
首先Makeroot(x)，然後Access(y)，此時x與y在同一棵splay中，再Splay(y)，此時x是y的左兒子，然後將其左兒子及左兒子的父母設成空節點即可。

	inline void cut(int x,int y){
	make_root(x);
	access(y);
	splay(y);
	node[y].ch[0]=0;
	node[x].pa=0;
	}

view raw lct_cut.cpp hosted with ❤ by GitHub

cut_parents(x):
跟cut類似，斷掉x和其父親節點之間的邊。首先access節點x。然後讓x旋轉到所在輔助樹的根節點上。斷掉左子樹。

	inline void cut_parents(int x){
	access(x);
	splay(x);
	node[node[x].ch[0]].pa=0;
	node[x].ch[0]=0;
	}

view raw lct_cut_parents.cpp hosted with ❤ by GitHub

find_root(x):
尋找點x在所在樹的根節點。這個操作很簡單，我們只要對x做一次access操作，這樣x和它的根就應該在一個splay中了。那麼此時的根就應該是這個splay最左邊的節點。
因為是splay tree，所以在搜尋完後仍要進行splay操作

	inline int find_root(int x){
	x=access(x);
	while(node[x].ch[0])x=node[x].ch[0];
	splay(x);
	return x;
	}

view raw lct_find_root.cpp hosted with ❤ by GitHub

對於很多問u到v的路徑問題，我們可以將u到v的路徑構造成一顆splay tree並求出其根節點維護的值即可。

	inline int query(int u,int v){
	/*
	傳回uv路徑splay tree的根結點
	這種寫法無法求LCA
	*/
	make_root(u);
	return access(v);
	}

view raw lct_query.cpp hosted with ❤ by GitHub

如果要不改變根節點的話可以求出LCA再進行處理，以下用u,v路徑上的點權總和為例:

	inline int query_lca(int u,int v){
	/*假設求鏈上點權的總和，sum是子樹的權重和，data是節點的權重*/
	access(u);
	int lca=access(v);
	splay(u);
	if(u==lca){
	return node[lca].data+node[node[lca].ch[1]].sum;
	}else{
	return node[lca].data+node[node[lca].ch[1]].sum+node[u].sum;
	}
	}

view raw lct_query_lca.cpp hosted with ❤ by GitHub

修改點權:

	inline void change(int x,int b){
	splay(x);
	node[x].data=b;
	up(x);
	}

view raw lct_change.cpp hosted with ❤ by GitHub

將一棵有邊權的樹構建成link cut tree:

	struct EDGE{
	int a,b,w;
	}e[10005];
	int n;
	std::vector<std::pair<int ,int > >G[10005];
	/*first表示子節點，second表示邊的編號*/
	int pa[10005],edge_node[10005];
	/*pa是父母節點，暫存用的，edge_node是每個編被存在哪個點裡面的陣列*/
	inline void bfs(int root){
	/*在建構的時候把每個點都設成一個splay tree，不會壞掉*/
	std::queue<int > q;
	for(int i=1;i<=n;++i)pa[i]=0;
	q.push(root);
	while(q.size()){
	int u=q.front();
	q.pop();
	for(int i=0;i<(int)G[u].size();++i){
	int v=G[u][i].first;
	if(v!=pa[u]){
	pa[v]=u;
	node[v].pa=u;
	node[v].data=e[G[u][i].second].w;
	edge_node[G[u][i].second]=v;
	up(v);
	q.push(v);
	}
	}
	}
	}

view raw lct_edge_build.cpp hosted with ❤ by GitHub

查詢這棵樹的邊權，因為不能改變root，所以要對access函數作一點修改，以下以u,v路徑上的最大值為例:

	inline void access(int x,bool is=0){/*is=0就是一般的access*/
	int last=0;
	while(x){
	splay(x);
	if(is&&!node[x].pa){
	printf("%d\n",max(node[last].ma,node[node[x].ch[1]].ma));
	}
	node[x].ch[1]=last;
	up(x);
	last=x;
	x=node[x].pa;
	}
	}

view raw lct_access2.cpp hosted with ❤ by GitHub

查詢只需要這樣寫就好了:

	inline void query_edge(int u,int v){
	access(u);
	access(v,1);
	}

view raw lct_query_edge.cpp hosted with ❤ by GitHub

[ Karatsuba Multiply ] 大數模板乘法Karatsuba法

複雜度大約是$\ord{N^{log3}}$，約為$\ord{N^{1.5}}$比一般直式乘法還快，附上連結:
Karatsuba算法

附模板，配合大數模板使用:

	//operator*的函數用以下這些函數取代
	bigN convert_base(int old_width,int new_width)const{
	vector<long long> p(max(old_width,new_width)+1,1);
	for(size_t i=1;i<p.size();++i)p[i]=p[i-1]*10;
	bigN ans;
	long long cur=0;
	int cur_id=0;
	for(size_t i=0;i<size();++i){
	cur+=at(i)*p[cur_id];
	cur_id+=old_width;
	while(cur_id>=new_width){
	ans.push_back(cur%p[new_width]);
	cur/=p[new_width];
	cur_id-=new_width;
	}
	}
	return ans.push_back(cur),ans.trim(),ans;
	}
	bigN karatsuba(const bigN &b)const{
	bigN res;res.resize(size()*2);
	if(size()<=32){
	for(size_t i=0;i<size();++i)
	for(size_t j=0;j<size();++j)
	res[i+j]+=at(i)*b[j];
	return res;
	}
	size_t k=size()/2;
	bigN a1(begin(),begin()+k);
	bigN a2(begin()+k,end());
	bigN b1(b.begin(),b.begin()+k);
	bigN b2(b.begin()+k,b.end());
	bigN a1b1=a1.karatsuba(b1);
	bigN a2b2=a2.karatsuba(b2);
	for(size_t i=0;i<k;++i)a2[i]+=a1[i];
	for(size_t i=0;i<k;++i)b2[i]+=b1[i];
	bigN r=a2.karatsuba(b2);
	for(size_t i=0;i<a1b1.size();++i)r[i]-=a1b1[i];
	for(size_t i=0;i<a2b2.size();++i)r[i]-=a2b2[i];
	for(size_t i=0;i<r.size();++i)res[i+k]+=r[i];
	for(size_t i=0;i<a1b1.size();++i)res[i]+=a1b1[i];
	for(size_t i=0;i<a2b2.size();++i)res[i+size()]+=a2b2[i];
	return res;
	}
	bigN operator*(const bigN &b)const{
	const static int mul_base=1000000,mul_width=log10(mul_base);
	bigN A=convert_base(width,mul_width);
	bigN B=b.convert_base(width,mul_width);
	int n=max(A.size(),B.size());
	while(n&(n-1))++n;
	A.resize(n),B.resize(n);
	bigN res=A.karatsuba(B);
	res.negative=negative!=b.negative;
	res.carry(mul_base);
	res=res.convert_base(mul_width,width);
	return res.trim(),res;
	}
	bigN operator*(long long b)const{//除法會用到O(N)
	bigN res=*this;
	if(b<0)res.negative=!negative,b=-b;
	for(size_t i=0,is=0;i<res.size()||is;++i){
	if(i==res.size())res.push_back(0);
	long long a=res[i]*b+is;
	is=a/base;
	res[i]=a%base;
	}
	return res.trim(),res;
	}

view raw karatsuba.cpp hosted with ❤ by GitHub

[ Augmenting Path Algorithm ] 二分圖匹配增廣路徑演算法

就是一直找增廣路徑，找到不能再找就是最大匹配惹
複雜度$\ord{n1*(m+n2+n1)}$，用dinic會比較快但是code比較長
模板:

	#define MAXN1 505
	#define MAXN2 505
	int n1,n2;/*n1個點連向n2個點*/
	int match[MAXN2];/*每個屬於n2的點匹配了哪個點*/
	vector<int > g[MAXN1];/*圖*/
	bool vis[MAXN2];/*是否走訪過*/
	bool dfs(int u){
	for(size_t i=0;i<g[u].size();++i){
	int v=g[u][i];
	if(vis[v])continue;
	vis[v]=1;
	if(match[v]==-1||dfs(match[v])){
	match[v]=u;
	return 1;
	}
	}
	return 0;
	}
	inline int max_match(){
	int ans=0;
	memset(match,-1,sizeof(int)*n2);
	for(int i=0;i<n1;++i){
	memset(vis,0,sizeof(bool)*n2);
	if(dfs(i))++ans;
	}
	return ans;
	}

view raw Augmenting Path.cpp hosted with ❤ by GitHub

[ Improved Shortest Augmenting Paths,ISAP ,Flow ] 網路流ISAP算法

使用gap優化(間隙優化)
使用當前弧優化

這裡提供遞迴與非遞迴兩種實現
基本上兩個版本效能差不多
遞迴版本:

	#include<string.h>
	#include<limits.h>
	#include<vector>
	#include<algorithm>
	template<typename T>
	struct ISAP{
	static const int MAXN=105;
	static const T INF=INT_MAX;
	int n;//點數
	int d[MAXN],gap[MAXN],cur[MAXN];
	//層次、gap[i]=層次為i的點之個數、當前弧優化
	struct edge{
	int v,pre;
	T cap,flow,r;
	edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
	};
	int g[MAXN];
	std::vector<edge> e;
	void init(int _n){
	memset(g,-1,sizeof(int)*((n=_n)+1));
	e.clear();
	}
	void add_edge(int u,int v,T cap,bool directed=false){
	e.push_back(edge(v,g[u],cap));
	g[u]=e.size()-1;
	e.push_back(edge(u,g[v],directed?0:cap));
	g[v]=e.size()-1;
	}
	T dfs(int u,int s,int t,T cur_flow=INF){
	if(u==t)return cur_flow;
	T tf=cur_flow,df;
	for(int &i=cur[u];~i;i=e[i].pre){
	if(e[i].r&&d[u]==d[e[i].v]+1){
	df=dfs(e[i].v,s,t,std::min(tf,e[i].r));
	e[i].flow+=df;
	e[i^1].flow-=df;
	e[i].r-=df;
	e[i^1].r+=df;
	if(!(tf-=df)||d[s]==n)return cur_flow-tf;
	}
	}
	int minh=n;
	for(int i=cur[u]=g[u];~i;i=e[i].pre){
	if(e[i].r&&d[e[i].v]<minh)minh=d[e[i].v];
	}
	if(!--gap[d[u]])d[s]=n;
	else ++gap[d[u]=++minh];
	return cur_flow-tf;
	}
	T isap(int s,int t,bool clean=true){
	memset(d,0,sizeof(int)*(n+1));
	memset(gap,0,sizeof(int)*(n+1));
	memcpy(cur,g,sizeof(int)*(n+1));
	if(clean){
	for(size_t i=0;i<e.size();++i){
	e[i].flow=0;
	e[i].r=e[i].cap;
	}
	}
	T max_flow=0;
	for(gap[0]=n;d[s]<n;)max_flow+=dfs(s,s,t);
	return max_flow;
	}
	std::vector<int> cut_e;//最小割邊集
	bool vis[MAXN];
	void dfs_cut(int u){
	vis[u]=1;//表示u屬於source的最小割集
	for(int i=g[u];~i;i=e[i].pre){
	if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
	}
	}
	T min_cut(int s,int t){
	T ans=isap(s,t);
	memset(vis,0,sizeof(bool)*(n+1));
	dfs_cut(s),cut_e.clear();
	for(int u=0;u<=n;++u){
	if(vis[u])for(int i=g[u];~i;i=e[i].pre){
	if(!vis[e[i].v])cut_e.push_back(i);
	}
	}
	return ans;
	}
	};

view raw isap_dfs.cpp hosted with ❤ by GitHub

非遞迴版本:

	#include<string.h>
	#include<limits.h>
	#include<vector>
	#include<algorithm>
	template<typename T>
	struct ISAP{
	static const int MAXN=105;
	static const T INF=INT_MAX;
	int n;//點數
	int d[MAXN],gap[MAXN],pre[MAXN],cur[MAXN];
	//層次、gap[i]=層次為i的點之個數、前一個節點、當前弧優化
	struct edge{
	int v,pre;
	T cap,flow,r;
	edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
	};
	int g[MAXN];
	std::vector<edge> e;
	void init(int _n){
	memset(g,-1,sizeof(int)*((n=_n)+1));
	e.clear();
	}
	void add_edge(int u,int v,T cap,bool directed=false){
	e.push_back(edge(v,g[u],cap));
	g[u]=e.size()-1;
	e.push_back(edge(u,g[v],directed?0:cap));
	g[v]=e.size()-1;
	}
	T isap(int s,int t,bool clean=true){
	T cur_flow,max_flow=0;
	int u=s;
	memset(d,0,sizeof(int)*(n+1));
	memset(gap,0,sizeof(int)*(n+1));
	memcpy(cur,g,sizeof(int)*(n+1));
	if(clean){
	for(size_t i=0;i<e.size();++i){
	e[i].flow=0;
	e[i].r=e[i].cap;
	}
	}
	for(gap[0]=n;d[s]<n;){
	if(u==t){
	cur_flow=INF;
	for(int i=s;i!=t;i=e[cur[i]].v){
	if(cur_flow>e[cur[i]].r)cur_flow=e[cur[i]].r;
	}
	for(int i=s;i!=t;i=e[cur[i]].v){
	e[cur[i]].flow+=cur_flow;
	e[cur[i]^1].flow-=cur_flow;
	e[cur[i]].r-=cur_flow;
	e[cur[i]^1].r+=cur_flow;
	}
	max_flow+=cur_flow;
	u=s;
	}
	int i;
	for(i=cur[u];~i;i=e[i].pre){
	if(e[i].r&&d[u]==d[e[i].v]+1)break;
	}
	if(~i){
	cur[u]=i;
	pre[e[i].v]=u;
	u=e[i].v;
	}else{
	if(!--gap[d[u]])break;
	for(d[u]=n,i=cur[u]=g[u];~i;i=e[i].pre){
	if(e[i].r&&d[e[i].v]<d[u])d[u]=d[e[i].v];
	}
	++gap[++d[u]];
	if(u!=s)u=pre[u];
	}
	}
	return max_flow;
	}
	std::vector<int> cut_e;//最小割邊集
	bool vis[MAXN];
	void dfs_cut(int u){
	vis[u]=1;//表示u屬於source的最小割集
	for(int i=g[u];~i;i=e[i].pre){
	if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
	}
	}
	T min_cut(int s,int t){
	T ans=isap(s,t);
	memset(vis,0,sizeof(bool)*(n+1));
	dfs_cut(s),cut_e.clear();
	for(int u=0;u<=n;++u){
	if(vis[u])for(int i=g[u];~i;i=e[i].pre){
	if(!vis[e[i].v])cut_e.push_back(i);
	}
	}
	return ans;
	}
	};

view raw isap.cpp hosted with ❤ by GitHub

[ Dinic's algorithm ,Flow ] 網路流Dinic模板

使用當前弧優化

多路增廣模板:

	#include<string.h>
	#include<limits.h>
	#include<vector>
	#include<queue>
	#include<algorithm>
	template<typename T>
	struct DINIC{
	static const int MAXN=105;
	static const T INF=INT_MAX;
	int n;/*點數*/
	int level[MAXN],cur[MAXN];/*層次、當前弧優化*/
	struct edge{
	int v,pre;
	T cap,flow,r;
	edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
	};
	int g[MAXN];
	std::vector<edge> e;
	void init(int _n){
	memset(g,-1,sizeof(int)*((n=_n)+1));
	e.clear();
	}
	void add_edge(int u,int v,T cap,bool directed=false){
	e.push_back(edge(v,g[u],cap));
	g[u]=e.size()-1;
	e.push_back(edge(u,g[v],directed?0:cap));
	g[v]=e.size()-1;
	}
	int bfs(int s,int t){
	memset(level,0,sizeof(int)*(n+1));
	memcpy(cur,g,sizeof(int)*(n+1));
	std::queue<int >q;
	q.push(s);
	level[s]=1;
	while(q.size()){
	int u=q.front();q.pop();
	for(int i=g[u];~i;i=e[i].pre){
	if(!level[e[i].v]&&e[i].r){
	level[e[i].v]=level[u]+1;
	q.push(e[i].v);
	if(e[i].v==t)return 1;
	}
	}
	}
	return 0;
	}
	T dfs(int u,int t,T cur_flow=INF){
	if(u==t||!cur_flow)return cur_flow;
	T df,tf=0;
	for(int &i=cur[u];~i;i=e[i].pre){
	if(level[e[i].v]==level[u]+1&&e[i].r){
	if(df=dfs(e[i].v,t,std::min(cur_flow,e[i].r))){
	e[i].flow+=df;
	e[i^1].flow-=df;
	e[i].r-=df;
	e[i^1].r+=df;
	tf+=df;
	if(!(cur_flow-=df))break;
	}
	}
	}
	if(!tf)level[u]=0;
	return tf;
	}
	T dinic(int s,int t,bool clean=true){
	if(clean){
	for(size_t i=0;i<e.size();++i){
	e[i].flow=0;
	e[i].r=e[i].cap;
	}
	}
	T ans=0;
	while(bfs(s,t))ans+=dfs(s,t);
	return ans;
	}
	std::vector<int> cut_e;//最小割邊集
	bool vis[MAXN];
	void dfs_cut(int u){
	vis[u]=1;//表示u屬於source的最小割集
	for(int i=g[u];~i;i=e[i].pre){
	if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
	}
	}
	T min_cut(int s,int t){
	T ans=dinic(s,t);
	memset(vis,0,sizeof(bool)*(n+1));
	dfs_cut(s),cut_e.clear();
	for(int u=0;u<=n;++u){
	if(vis[u])for(int i=g[u];~i;i=e[i].pre){
	if(!vis[e[i].v])cut_e.push_back(i);
	}
	}
	return ans;
	}
	};

view raw dinic.cpp hosted with ❤ by GitHub

單路增廣模板:

	#include<string.h>
	#include<limits.h>
	#include<vector>
	#include<queue>
	#include<algorithm>
	template<typename T>
	struct DINIC{
	static const int MAXN=105;
	static const T INF=INT_MAX;
	int n;/*點數*/
	int level[MAXN],cur[MAXN];/*層次、當前弧優化*/
	struct edge{
	int v,pre;
	T cap,flow,r;
	edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
	};
	int g[MAXN];
	std::vector<edge> e;
	void init(int _n){
	memset(g,-1,sizeof(int)*((n=_n)+1));
	e.clear();
	}
	void add_edge(int u,int v,T cap,bool directed=false){
	e.push_back(edge(v,g[u],cap));
	g[u]=e.size()-1;
	e.push_back(edge(u,g[v],directed?0:cap));
	g[v]=e.size()-1;
	}
	int bfs(int s,int t){
	memset(level,0,sizeof(int)*(n+1));
	memcpy(cur,g,sizeof(int)*(n+1));
	std::queue<int >q;
	q.push(s);
	level[s]=1;
	while(q.size()){
	int u=q.front();q.pop();
	for(int i=g[u];~i;i=e[i].pre){
	if(!level[e[i].v]&&e[i].r){
	level[e[i].v]=level[u]+1;
	q.push(e[i].v);
	if(e[i].v==t)return 1;
	}
	}
	}
	return 0;
	}
	T dfs(int u,int t,T cur_flow=INF){
	if(u==t)return cur_flow;
	T df;
	for(int &i=cur[u];~i;i=e[i].pre){
	if(level[e[i].v]==level[u]+1&&e[i].r){
	if(df=dfs(e[i].v,t,std::min(cur_flow,e[i].r))){
	e[i].flow+=df;
	e[i^1].flow-=df;
	e[i].r-=df;
	e[i^1].r+=df;
	return df;
	}
	}
	}
	return level[u]=0;
	}
	T dinic(int s,int t,bool clean=true){
	if(clean){
	for(size_t i=0;i<e.size();++i){
	e[i].flow=0;
	e[i].r=e[i].cap;
	}
	}
	T ans=0,mf=0;
	while(bfs(s,t))while(mf=dfs(s,t))ans+=mf;
	return ans;
	}
	std::vector<int> cut_e;//最小割邊集
	bool vis[MAXN];
	void dfs_cut(int u){
	vis[u]=1;//表示u屬於source的最小割集
	for(int i=g[u];~i;i=e[i].pre){
	if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
	}
	}
	T min_cut(int s,int t){
	T ans=dinic(s,t);
	memset(vis,0,sizeof(bool)*(n+1));
	dfs_cut(s),cut_e.clear();
	for(int u=0;u<=n;++u){
	if(vis[u])for(int i=g[u];~i;i=e[i].pre){
	if(!vis[e[i].v])cut_e.push_back(i);
	}
	}
	return ans;
	}
	};

view raw dinic2.cpp hosted with ❤ by GitHub

[ Edmonds–Karp algorithm ,Flow] 網路流Edmonds–Karp算法

這是最基本的多項式時間網路流演算法，因此在這裡我多了一些變化，首先是一般用鄰接矩陣，複雜度$\ord{\abs{V}^5}$的做法：

	#include<string.h>
	#include<limits.h>
	#include<queue>
	#include<algorithm>
	#define MAXN 100
	int n,m;/*點數、邊數*/
	int g[MAXN+5][MAXN+5],f[MAXN+5][MAXN+5],r[MAXN+5][MAXN+5];
	/*容量上限、流量、剩餘容量*/
	int v[MAXN+5];/*經過的點*/
	int pa[MAXN+5];/*紀錄路徑*/
	int d[MAXN+5];/*紀錄流量瓶頸*/
	inline int bfs(int s,int t){
	memset(v,0,sizeof(v));
	std::queue<int >q;
	q.push(s);
	v[s]=1;
	pa[s]=-1;
	d[s]=INT_MAX;
	while(q.size()){
	int x=q.front();q.pop();
	for(int i=1;i<=n;++i){
	if(v[i]||!r[x][i])continue;
	v[i]=1;
	pa[i]=x;
	d[i]=std::min(d[x],r[x][i]);
	q.push(i);
	if(i==t)return d[t];
	}
	}
	return 0;
	}
	inline int edmonds_karp(int s,int t){
	memset(f,0,sizeof(f));
	memcpy(r,g,sizeof(g));
	int ans=0,tmd;
	for(;(tmd=bfs(s,t));ans+=tmd){
	for(int i=t;~pa[i];i=pa[i]){
	f[pa[i]][i]+=tmd;
	f[i][pa[i]]-=tmd;
	r[i][pa[i]]=g[i][pa[i]]-f[i][pa[i]];
	r[pa[i]][i]=g[pa[i]][i]-f[pa[i]][i];
	}
	}
	return ans;
	}

view raw Edmonds Karp.cpp hosted with ❤ by GitHub

接著是使用鄰接串列，複雜度是$\ord{\abs{E}^2\abs{V}}$的作法：

	#include<cstring>
	#include<climits>
	#include<queue>
	#include<algorithm>
	template<typename T>
	struct Edmonds_Karp{
	static const int MAXN=105;
	static const T INF=INT_MAX;
	struct edge{
	int v,pre;
	T cap,r;
	edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),r(cap){}
	};
	int g[MAXN],n;
	bool vis[MAXN];
	int pa[MAXN];
	T d[MAXN];
	std::vector<edge> e;
	void init(int _n){
	memset(g,-1,sizeof(int)*((n=_n)+1));
	e.clear();
	}
	void add_edge(int u,int v,T cap,bool directed=false){
	e.push_back(edge(v,g[u],cap));
	g[u]=e.size()-1;
	e.push_back(edge(u,g[v],directed?0:cap));
	g[v]=e.size()-1;
	}
	T bfs(int s,int t){
	memset(vis,0,sizeof(bool)*(n+1));
	std::queue<int> q;
	q.push(s);
	vis[s]=1;
	pa[s]=-1;
	d[s]=INF;
	while(q.size()){
	int u=q.front();q.pop();
	for(int i=g[u];~i;i=e[i].pre){
	if(vis[e[i].v]||e[i].r==0)continue;
	vis[e[i].v]=1;
	pa[e[i].v]=i;
	d[e[i].v]=std::min(d[u],e[i].r);
	q.push(e[i].v);
	if(e[i].v==t)return d[t];
	}
	}
	return 0;
	}
	T edmonds_karp(int s,int t){
	T ans=0,tmd;
	for(;(tmd=bfs(s,t));ans+=tmd){
	for(int u=t;~pa[u];){
	int i=pa[u];
	e[i].r-=tmd;
	e[i^1].r+=tmd;
	u=e[i^1].v;
	}
	}
	return ans;
	}
	};

view raw Edmonds Karp3.cpp hosted with ❤ by GitHub

最近幫學長寫作業的時候，又發現了一個不需要多紀錄反邊的做法，但是每個點要同時記錄入邊和出邊，這樣可以用較少的記憶體完成演算法：

	#include<cstring>
	#include<climits>
	#include<queue>
	#include<algorithm>
	template<typename T>
	struct Edmonds_Karp{
	static const int MAXN=105;
	static const T INF=INT_MAX;
	struct Edge{
	int u,v;
	T cap,r;
	Edge(int u,int v,const T &cap):u(u),v(v),cap(cap),r(cap){}
	};
	std::vector<Edge> edge;
	std::vector<int> in[MAXN];
	std::vector<int> out[MAXN];
	bool vis[MAXN];
	int pa[MAXN];
	T d[MAXN];
	int n;
	void init(int _n){
	edge.clear();
	for(n=_n++;_n--;){
	in[_n].clear();
	out[_n].clear();
	}
	}
	void add_edge(int u,int v,const T &cap){
	in[u].emplace_back(edge.size());
	out[v].emplace_back(edge.size());
	edge.emplace_back(u,v,cap);
	}
	T bfs(int s,int t){
	memset(vis,0,sizeof(bool)*(n+1));
	std::queue<int> q;
	q.push(s);
	vis[s]=1;
	pa[s]=-1;
	d[s]=INF;
	while(q.size()){
	int u=q.front();q.pop();
	for(auto i:in[u]){
	const auto &e=edge[i];
	if(vis[e.v]||e.r==0)continue;
	vis[e.v]=1;
	pa[e.v]=i;
	d[e.v]=std::min(d[u],e.r);
	q.push(e.v);
	if(e.v==t)return d[t];
	}
	for(auto i:out[u]){
	const auto &e=edge[i];
	if(vis[e.u]||e.r==e.cap)continue;
	vis[e.u]=1;
	pa[e.u]=i;
	d[e.u]=std::min(d[u],e.cap-e.r);
	q.push(e.u);
	if(e.u==t)return d[t];
	}
	}
	return 0;
	}
	T edmonds_karp(int s,int t){
	T ans=0,tmd;
	for(;(tmd=bfs(s,t));ans+=tmd){
	for(int u=t;~pa[u];){
	auto &e=edge[pa[u]];
	if(e.v==u){
	e.r-=tmd;
	u=e.u;
	}else{
	e.r+=tmd;
	u=e.v;
	}
	}
	}
	return ans;
	}
	};

view raw Edmonds Karp2.cpp hosted with ❤ by GitHub

[ Travelling Salesman Problem, TSP ] 旅行推銷員問題

假設圖是完全圖的情況下，從原點出發經過所有點一次並回到原點的最短路徑問題
以下為模板:

	#include<limits.h>
	#include<algorithm>
	#define MAXN 16
	int n;/*點數*/
	int dp[MAXN][1<<MAXN];/*DP陣列*/
	int g[MAXN][MAXN];/*圖*/
	int dfs(int s,int d){
	if(dp[d][s])return dp[d][s];
	if((1<<d)==s)return g[0][d];
	dp[d][s]=INT_MAX;
	for(int i=0;i<n;++i){
	if(((1<<i)&s)&&i!=d){
	dp[d][s]=std::min(dp[d][s],dfs(s^(1<<d),i)+g[i][d]);
	}
	}
	return dp[d][s];
	}

view raw TSP.cpp hosted with ❤ by GitHub

[ print itself ] 印出自己的code

一個會輸出自己的code:

	#include<stdio.h>
	#define s "#include<stdio.h>%c#define s %c%s%c%cint main(){%c%cprintf(s,10,34,s,34,10,10,9,10,9,10);%c%creturn 0;%c}"
	int main(){
	printf(s,10,34,s,34,10,10,9,10,9,10);
	return 0;
	}

view raw print self.cpp hosted with ❤ by GitHub

可以看一下輸出:http://codepad.org/ZFGeRjkn

[ Rabin-Karp rolling hash ] Rabin-Karp 字串hash演算法

Michael O. Rabin和Richard M. Karp在1987年提出一個想法，即可以對模式串進行哈希運算並將其哈希值與文本中子串的哈希值進行比對。
設字串陣列為$S[n]$，字元編號為$0$到$n-1$
首先建立陣列$Hash[n+1]$，$Hash[i]$表示
$(S[0]*P^{i-1}+S[1]*P^{i-2}+...+S[i-1])\%PM,Hash[0]=0$,其中$P$跟$PM$是質數
假設要取編號$L$到編號$R$左閉右開區間的hash值，則其為$(Hash[R]-Hash[L]*P^{R-L})\%PM$
如果要用double hash則用不同質數計算兩種不同的hash值，用pair拼起來即可
以下提供模板:

	template <unsigned P, unsigned PM> class RollingHash {
	static vector<unsigned> Base;
	vector<unsigned> Hash;
	using LL = long long;
	public:
	RollingHash(const string &S) : Hash(S.size() + 1) {
	if (Base.empty())
	Base.resize(1, 1);
	for (size_t i = 1; i <= S.size(); ++i)
	Hash[i] = (LL(Hash[i - 1]) * P + S[i - 1]) % PM;
	while (Base.size() < Hash.size())
	Base.emplace_back((LL(Base.back()) * P) % PM);
	}
	size_t size() const { return Hash.size() - 1; }
	const vector<unsigned> &getHash() const { return Hash; }
	unsigned hash() const { return Hash.back(); }
	unsigned hash(unsigned L, unsigned R) const {
	return (Hash[R] + PM - (LL(Hash[L]) * Base[R - L]) % PM) % PM;
	}
	};
	template <unsigned P, unsigned PM> vector<unsigned> RollingHash<P, PM>::Base;

view raw rolling_hash.cpp hosted with ❤ by GitHub

[ Fraction ] 分數模板

今天學弟的比賽剛結束，因為有考到相關內容所以到現在才釋出分數模板。
因為每種題目要求輸出的分數格式都不一樣，所以不提供輸出的程式，必定約到最簡分數。
以下提供模板:

	/************************
	n為分子，d為分母
	若分數為0則n=0,d=1
	若為負數則負號加在分子
	必定約到最簡分數
	************************/
	#ifndef SUNMOON_FRACTION
	#define SUNMOON_FRACTION
	#include<algorithm>
	template<typename T>
	struct fraction{
	T n,d;
	fraction(const T &_n=0,const T &_d=1):n(_n),d(_d){
	T t=std::__gcd(n,d);
	n/=t,d/=t;
	if(d<0)n=-n,d=-d;
	}
	fraction operator-()const{
	return fraction(-n,d);
	}
	fraction operator+(const fraction &b)const{
	return fraction(n*b.d+b.n*d,d*b.d);
	}
	fraction operator-(const fraction &b)const{
	return fraction(n*b.d-b.n*d,d*b.d);
	}
	fraction operator*(const fraction &b)const{
	return fraction(n*b.n,d*b.d);
	}
	fraction operator/(const fraction &b)const{
	return fraction(n*b.d,d*b.n);
	}
	fraction operator+=(const fraction &b){
	return *this=fraction(n*b.d+b.n*d,d*b.d);
	}
	fraction operator-=(const fraction &b){
	return *this=fraction(n*b.d-b.n*d,d*b.d);
	}
	fraction operator*=(const fraction &b){
	return *this=fraction(n*b.n,d*b.d);
	}
	fraction operator/=(const fraction &b){
	return *this=fraction(n*b.d,d*b.n);
	}
	bool operator <(const fraction &b)const{
	return n*b.d<b.n*d;
	}
	bool operator >(const fraction &b)const{
	return n*b.d>b.n*d;
	}
	bool operator ==(const fraction &b)const{
	return n*b.d==b.n*d;
	}
	bool operator <=(const fraction &b)const{
	return n*b.d<=b.n*d;
	}
	bool operator >=(const fraction &b)const{
	return n*b.d>=b.n*d;
	}
	};
	#endif

view raw fraction.cpp hosted with ❤ by GitHub

[ shortest path algorithm ] floyd,dijkstra,bellman,spfa 最短路徑模板

Floyd-Warshall:

	int n;
	int d[105][105];
	int g[105][105];
	/*************************
	點數
	答案
	圖，如果兩點間不存在邊則為INF
	*************************/
	inline void floyd(){
	static int i,j,k;
	for(i=0;i<n;++i){
	for(j=0;j<n;++j){
	d[i][j]=g[i][j];
	}
	}
	for(k=0;k<n;++k){
	for(i=0;i<n;++i){
	for(j=0;j<n;++j){
	if(d[i][j]>d[i][k]+d[k][j]){
	d[i][j]=d[i][k]+d[k][j];
	}
	}
	}
	}
	}

view raw Floyd-Warshall.cpp hosted with ❤ by GitHub

Dijkstra:

	const int MAXN=50005;
	const long long INF=0x3f3f3f3f3f3f3f3f;
	typedef pair<long long ,int > PLI;
	typedef pair<int ,long long > PIL;
	int n,m;//點數、邊數
	vector<PIL> g[MAXN];//圖
	long long d[MAXN];//答案
	inline bool dijkstra(int from,int to){
	priority_queue<PLI,vector<PLI>,greater<PLI> >q;
	for(int i=1;i<=n;++i)d[i]=INF;
	d[from]=0;
	q.push(make_pair(d[from],from));
	while(q.size()){
	PLI u=q.top();
	q.pop();
	int x=u.second;
	if(d[x]<u.first)continue;
	for(size_t i=0;i<g[x].size();++i){
	if(d[g[x][i].first]>d[x]+g[x][i].second){
	d[g[x][i].first]=d[x]+g[x][i].second;
	q.push(make_pair(d[g[x][i].first],g[x][i].first));
	}
	}
	}
	return d[to]==INF?0:1;//有/無解
	}

view raw Dijkstra.cpp hosted with ❤ by GitHub

Bellman-Ford:

	#define x first
	#define y second
	const long long INF=0x3f3f3f3f3f3f3f3f;
	const int MAXN=1005,MAXM=100005;
	int n,m;//點、邊的個數
	pair<pair<int,int > ,int > e[MAXM];//所有邊
	long long d[MAXN];//答案
	int bellman(int from,int to){
	for(int i=0;i<n;++i)d[i]=INF;
	d[from]=0;
	for(int t=1;t<=n;++t){
	bool fix=0;
	for(int i=0;i<m;++i){
	if(d[e[i].x.x]!=INF&&d[e[i].x.y]>d[e[i].x.x]+e[i].y){
	if(e[i].x.y==to&&t==n)return -1;//從from到to會經過負環
	d[e[i].x.y]=d[e[i].x.x]+e[i].y;
	fix=1;
	}
	}
	if(!fix)break;
	}
	return d[to]!=INF;//有/無解
	}

view raw Bellman-Ford.cpp hosted with ❤ by GitHub

SPFA:

	const int MAXN=50005;
	const long long INF=0x3f3f3f3f3f3f3f3f;
	typedef pair<int,long long > PIL;
	int n,m;//點數、邊數
	vector<PIL> g[MAXN];//圖
	long long d[MAXN];//答案
	bool inq[MAXN];//是否在queue裡面
	int inq_cnt[MAXN];//進queue的次數
	inline int spfa(int from,int to){
	queue<int > q;
	for(int i=1;i<=n;++i)d[i]=INF,inq[i]=inq_cnt[i]=0;
	d[from]=0;
	inq[from]=inq_cnt[from]=1;
	q.push(from);
	while(q.size()){
	int o=q.front();
	q.pop();
	if(inq_cnt[o]>n)return -1;//存在負環
	inq[o]=0;
	for(size_t i=0;i<g[o].size();++i){
	if(d[g[o][i].first]>d[o]+g[o][i].second){
	d[g[o][i].first]=d[o]+g[o][i].second;
	if(!inq[g[o][i].first]){
	inq[g[o][i].first]=1;
	++inq_cnt[g[o][i].first];
	q.push(g[o][i].first);
	}
	}
	}
	}
	return d[to]!=INF;//有/無解
	}

view raw SPFA.cpp hosted with ❤ by GitHub

[ Ternary Search ] 三分搜尋法模板

對於一個U型或倒U型函數，我們可以利用三分搜尋法找出其最低/最高點，時間複雜度為$\ord{log \; n}$,這裡n是$10^{(整數位數+浮點位數)}$

關於原理就不多說了，以下提供對於U型函數的模板:

	/*
	f是傳入的函數，eps表示精準度，通常設為1e-5
	l,r為邊界，所以必須先預估出極值點的位置
	*/
	template<typename _F>
	inline double ternary_search(_F f,double l,double r,double eps){
	static double m1,m2;
	while(r-l>eps){
	m1=l+(r-l)/3;
	m2=r-(r-l)/3;
	if(f(m1)>f(m2))l=m1;/*如果對象是倒U型函數，則改成 < */
	else r=m2;
	}
	return (l+r)/2;
	}

view raw ternary_search.cpp hosted with ❤ by GitHub

用法:

	inline double f(double x){
	return x*x-2*x+1;
	}
	/*假設邊界為-10~10*/
	printf("%.5f\n",ternary_search(f,-10,10,1e-6));

view raw use_way.cpp hosted with ❤ by GitHub

[ Lowest Common Ancestor , LCA Euler Tour Technique ] 最近共同祖先樹壓平轉RMQ算法

想法就是把樹壓平，然後紀錄當前節點u的深度在dep[u]，紀錄時間戳記，進入的時候要紀錄一次，出來的時候也要記錄，把當前點u進入的時間戳紀錄在in[u]中，時間戳的大小為2*n，所以必須用兩倍的記憶體來維護。之後就轉成RMQ問題
查詢a,b兩點的LCA=dep[in[a]]到dep[in[b]]中，最小深度的時間戳經過的那個點

	#define MAXN 100000
	typedef vector<int >::iterator VIT;
	int dep[MAXN+5],in[MAXN+5];
	int vs[2*MAXN+5];
	int cnt;/*時間戳*/
	vector<int >G[MAXN+5];
	void dfs(int x,int pa){
	in[x]=++cnt;
	vs[cnt]=x;
	for(VIT i=G[x].begin();i!=G[x].end();++i){
	if(*i==pa)continue;
	dep[*i]=dep[x]+1;
	dfs(*i,x);
	vs[++cnt]=x;
	}
	}

view raw Euler Tour Technique.cpp hosted with ❤ by GitHub

查詢:

	inline int find_lca(int a,int b){
	if(in[a]>in[b])swap(a,b);
	return RMQ(in[a],in[b]);
	}

view raw query.cpp hosted with ❤ by GitHub

RMQ的部分就要自己維護了，使用線段樹效能會非常差，在n,q<=100000的情況下使用稀疏表的效能與使用樹鏈剖分的效能差不多

[ Lowest Common Ancestor , LCA Doubling Search ] 最近共同祖先倍增演算法

首先建立一張倍增表，pa[i][x]表示x的第$2^i$的祖先，若不存在則為$-1$，可以在DFS或BFS的過程中求得，複雜度為$\ord{N*MAX\_LOG}$

對於每筆詢問，我們可以利用剛剛求得的倍增表，在$\ord{MAX\_LOG}$的時間在線查詢
先將要查詢的a,b兩點較深的點移至較淺的點的深度，之後二分搜LCA的位置
附上模板:

	#define MAXN 100000
	#define MAX_LOG 17
	typedef vector<int >::iterator VIT;
	int pa[MAX_LOG+1][MAXN+5],dep[MAXN+5];
	vector<int >G[MAXN+5];
	void dfs(int x,int p){
	pa[0][x]=p;
	for(int i=0;i+1<MAX_LOG;++i){
	pa[i+1][x]=pa[i][pa[i][x]];
	}
	for(VIT i=G[x].begin();i!=G[x].end();++i){
	if(*i==p)continue;
	dep[*i]=dep[x]+1;
	dfs(*i,x);
	}
	}
	inline int find_lca(int a,int b){
	if(dep[a]>dep[b])swap(a,b);
	for(int i=dep[b]-dep[a],k=0;i;i/=2){
	if(i&1)b=pa[k][b];
	++k;
	}
	if(a==b)return a;
	for(int i=MAX_LOG;i>=0;--i){
	if(pa[i][a]!=pa[i][b]){
	a=pa[i][a];
	b=pa[i][b];
	}
	}
	return pa[0][a];
	}

view raw Doubling Search.cpp hosted with ❤ by GitHub

[ Lowest Common Ancestor , LCA Tarjan's Algorithm ] 最近共同祖先tarjan演算法

tarjan算法的步驟是(當dfs到節點u時):

在並查集中建立僅有u的集合,設置該集合的祖先為u
對u的每個孩子v:

1. 遞迴處理v
2. 合併v的集合到父節點u的集合,確保集合的祖先是u

設置u為已遍歷
處理關於u的查詢,若查詢(u,x)中的x已遍歷過,則LCA(u,x) = x所在的集合的祖先

因為只進行過一次DFS，所以複雜度為$\ord{(N+Q)*\alpha(N)}$，為了優化並查集，這裡採用啟發式合併，因此為了記錄集合的祖先而外增加了一個數組head
其雖為線性，但是效率並不怎麼好，跟倍增法的效能差不多

	#define MAXN 100000
	#define MAXQ 100000
	typedef vector<int >::iterator VIT;
	typedef vector<pair<int,int> >::iterator VPIT;
	int tr[MAXN+5],siz[MAXN+5];/*並查集，dfs前必須初始化*/
	int head[MAXN+5],ans[MAXQ+5];
	bool v[MAXN+5];/*走訪標記*/
	vector<int >G[MAXN+5];
	vector<pair<int,int> >Q[MAXN+5];
	inline void init(int n){/*初始化並查集*/
	for(int i=1;i<=n;++i){
	tr[i]=i;
	siz[i]=1;
	}
	}
	int find(int x){
	return tr[x]==x?x:tr[x]=find(tr[x]);
	}
	inline void merge(int a,int b){
	a=find(a);
	b=find(b);
	if(a==b)return;
	if(siz[a]<siz[b])swap(a,b);
	tr[b]=a;
	siz[a]+=siz[b];
	}
	void dfs(int u,int pa){
	head[u]=u;
	for(VIT i=G[u].begin();i!=G[u].end();++i){
	if(*i==pa)continue;
	dfs(*i,u);
	merge(u,*i);
	head[find(u)]=u;/*記錄u所在集合的頭為u*/
	}
	v[u]=1;
	for(VPIT i=Q[u].begin();i!=Q[u].end();++i){
	if(v[i->first])ans[i->second]=head[find(i->first)];
	}
	}

view raw LCA Tarjan.cpp hosted with ❤ by GitHub

使用方法也很簡單，只需要將每個詢問的兩點記錄起來即可

	int n,m,a,b;
	scanf("%d%d",&n,&m);
	for(int i=1;i<n;++i){
	scanf("%d%d",&a,&b);
	G[a].push_back(b);
	G[b].push_back(a);
	}
	for(int i=0;i<m;++i){/*紀錄詢問*/
	scanf("%d%d",&a,&b);
	Q[a].push_back(make_pair(b,i));
	Q[b].push_back(make_pair(a,i));
	}
	init(n);
	dfs(1,-1);/*假設編號為1~n，root=1*/
	for(int i=0;i<m;++i)printf("%d\n",ans[i]);

view raw query.cpp hosted with ❤ by GitHub

[ Mo's algorithm ]莫隊算法

對於一些詢問區間答案的資料結構題，我們已經知道了區間[l,r]的答案，且能在極少的時間( $\ord 1$ or $\ord{logN}$ )內得到區間[l+1,r]、[l-1,r]、[l,r+1]、[l,r-1]的答案，題目也准許離線，只要滿足這些性質就可以使用莫隊。

首先我們須將所有詢問記錄下來，將序列分成sqrt(n)塊，找出每筆詢問l所在的塊，每塊內r由小排到大

	#define MAXQ 100000
	struct Q{
	int l,r,i,block;
	inline bool operator<(const Q &b)const{
	return block==b.block?r<b.r:block<b.block;
	}
	}q[MAXQ+5];

view raw struct.cpp hosted with ❤ by GitHub

	lim=1+(int)sqrt(n);
	for(int i=0;i<m;++i){
	scanf("%d%d",&q[i].l,&q[i].r);
	q[i].block=q[i].l/lim;
	q[i].i=i;
	}
	sort(q,q+m);

view raw scanf query.cpp hosted with ❤ by GitHub

這裡n是序列的長度

之後就按造順序一個一個將序列元素在一些維護費用很小的資料結構進行插入刪除，直到等於當前區間為止，此時就可以把這個資料結構的答案存到ans陣列裡，然後結束後一次輸出

	for(int i=0,L=1,R=0;i<m;++i){/*這是閉區間的寫法*/
	while(R<q[i].r)add(s[++R]);
	while(L>q[i].l)add(s[--L]);
	while(R>q[i].r)sub(s[R--]);
	while(L<q[i].l)sub(s[L++]);
	ans[q[i].i]=cur_ans;
	}
	for(int i=0;i<m;++i)printf("%d\n",ans[i]);

view raw find query.cpp hosted with ❤ by GitHub

這裡進行一下簡單的證明:

• 假設序列長度為N，我們將序列每K個分成一塊，在排完序後，左界屬於同一塊的詢問，每次都不會移動操過K次；而左界屬於不同塊的操作，從第i塊移動到第i+1塊最多要移動2*K次，總共分成N/K塊，故移動次數加起來不操過2*K*(N/K)=2*N次，因此左界移動的複雜度為$\ord{Q*K}+\ord N$，Q為詢問次數
• 對於右界，在同一塊內因為是遞增的關係所以移動次數不會操過N次，而最多只有N/K塊，因此複雜度為$\ord{N*N/K}$
• 這個時候我們可以取K=sqrt(N)得到一個複雜度為$\ord{Q*K+N+N*N/K}=\ord{(Q+N)*sqrt(N)}$的方法

這是區間眾數的模板(對於每個詢問l,r輸出區間[l,r]眾數的個數):

	#include<bits/stdc++.h>
	using namespace std;
	#define MAXN 100000
	#define MAXQ 100000
	struct Q{
	int l,r,i,block;
	inline bool operator<(const Q &b)const{
	return block==b.block?r<b.r:block<b.block;
	}
	}q[MAXQ+5];
	int s[MAXN+5],lim;
	int ans[MAXQ+5];
	int u[MAXN+5],cnt[MAXN+5],cur_ans;
	inline void add(int x){
	x=++u[x];
	--cnt[x-1];
	++cnt[x];
	cur_ans=max(cur_ans,x);
	}
	inline void sub(int x){
	x=--u[x];
	--cnt[x+1];
	++cnt[x];
	if(!cnt[cur_ans])--cur_ans;
	}
	int n,m;
	int p[MAXN+5],nn;
	int main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;++i){
	scanf("%d",&s[i]);
	p[i-1]=s[i];
	}
	sort(p,p+n);
	nn=unique(p,p+n)-p;
	for(int i=1;i<=n;++i)s[i]=lower_bound(p,p+nn,s[i])-p;
	lim=1+(int)sqrt(n);
	for(int i=0;i<m;++i){
	scanf("%d%d",&q[i].l,&q[i].r);
	q[i].block=q[i].l/lim;
	q[i].i=i;
	}
	sort(q,q+m);
	for(int i=0,L=1,R=0;i<m;++i){/*這是閉區間的寫法*/
	while(R<q[i].r)add(s[++R]);
	while(L>q[i].l)add(s[--L]);
	while(R>q[i].r)sub(s[R--]);
	while(L<q[i].l)sub(s[L++]);
	ans[q[i].i]=cur_ans;
	}
	for(int i=0;i<m;++i)printf("%d\n",ans[i]);
	return 0;
	}

view raw range mode.cpp hosted with ❤ by GitHub

[ sparse table ] RMQ問題的稀疏表算法

RMQ問題(Range Minimum/Maximum Query,區間最值問題)，通常是給定一條串列s，有若干次的詢問查詢區間L至區間R的最大(小)值。
相信關於RMQ問題使用線段樹的做法很多人都有一定的了解，這邊提供另外一種在線不可修改的做法，稱為稀疏表算法，其預處理時間為$\ord{NlogN}$，N為串列長度，
查詢時間為$\ord 1$，是一種高效的演算法。

關於此算法的說明請看http://noalgo.info/489.html
由連結提供的程式碼效能是十分差的，因此這邊提供我的模板(以區間最小值為例):

	#define MAXN 100000
	#define MAX_LOG 17
	int n,s[MAXN+5];
	int st[MAX_LOG+1][MAXN+5];
	inline void init(){/*假設區間由[0~n-1]*/
	for(int i=0;i<n;++i)st[0][i]=s[i];
	for(int j=1;(1<<j)<=n;++j)
	for(int i=0;i+(1<<j)<=n;++i)
	st[j][i]=min(st[j-1][i],st[j-1][i+(1<<(j-1))]);
	}
	inline int RMQ(int a,int b){/*用<algorithm>內建函式求log*/
	int k=std::__lg(b-a+1);
	return min(st[k][a],st[k][b-(1<<k)+1]);
	}

view raw Sparse Table.cpp hosted with ❤ by GitHub

[ Palindromic Tree ] 回文樹(回文自動機)

回文樹是去年新出的資料結構，由一個俄羅斯codeforces的紅人Mikhail Rubinchik在2014年發明的(http://codeforces.com/blog/entry/13959)，雖然名子直接翻譯叫回文樹，但他長的比較像一個自動機，所以也有很多人翻回文自動機。
在構造它之前，我們必須先證明一個長度n的字串其不同的回文子字串個數<=n，什麼是不同的回文子字串?就是長的不同的回文子字串，及把出現位置不同但形狀一樣的回文子字串當做是同樣的子字串。

證明:

1. 由左往右一個一個增加字符，則新產生的回文子字串其結尾一定是當前增加的字符x
2. 考慮以當前位置為結尾的最長回文x......x，顯而易見的，若產生除了x......x之外其它的回文子串，則必定被x......x所包含
3. 若增加字符x後，產生了除了x......x之外其它的回文子串，則根據回文的性質，其一定早已在x......的部分理出現(例:假設S是對稱點x...S.xax，因為回文的性質所以S的另一邊也會出現相同的字串xax.S.xax)
4. 因此每增加一個字符最多只會產生一個新的回文子字串，而長度為n的字符串最多只會產生n種不同的回文子字串

接下來進行回文自動機元素的詳細解說:

回文自動機包含以下元素:

1. 狀態St，所有節點的集合，一開始兩個節點，0表示偶數長度串的根和1表示奇數長度串的根
2. last 新增加一個字符後所形成的最長回文串的節點編號
3. s 當前的字符串(一開始設s[0]=-1(可以是任意一個在串S中不會出現的字符))
4. n 表示添加的字符個數

每個節點代表一個不同的回文子字串，我們在每個節點會儲存一些數值:

1. len 表示所代表的回文子字串長度
2. next[c] 表示所代表的回文子字串在頭尾各增加一個字符c後的回文字串其節點編號
3. fail 表示所代表的回文子字串不包括本身的最長後綴回文子串的節點編號
4. cnt(非必要) 表示所代表的回文子字串在整體字串出現的次數(在建構完成後呼叫count()才能計算)
5. num(非必要) 表示所代表的回文子字串其後綴為回文字串的個數

注意一開始必須先將狀態初始化St[0].len=0，St[1].len=-1，last=0，s[0]=-1，n=0

關於構造方法及圖片說明請參考:http://blog.csdn.net/u013368721/article/details/42100363
概念請參考:http://blog.csdn.net/MaticsL/article/details/43651169
證明請參考:http://yqcmmd.com/index.php/2015/03/30/746/

最後提供模板(使用STL vector)，因為St的元素在解題的過程中常會被用到所以就不封裝了:

	#ifndef PALINDROMIC_TREE
	#define PALINDROMIC_TREE
	#include<vector>
	struct palindromic_tree{
	struct node{
	int next[26],fail,len;/*這些是必要的元素*/
	int cnt,num;/*這些是額外維護的元素*/
	node(int l=0):fail(0),len(l),cnt(0),num(0){
	for(int i=0;i<26;++i)next[i]=0;
	}
	};
	std::vector<node >St;
	std::vector<char >s;
	int last,n;
	palindromic_tree():St(2),last(1),n(0){
	St[0].fail=1;
	St[1].len=-1;
	s.push_back(-1);
	}
	inline void clear(){
	St.clear();
	s.clear();
	last=1;
	n=0;
	St.push_back(0);
	St.push_back(-1);
	St[0].fail=1;
	s.push_back(-1);
	}
	inline int get_fail(int x){
	while(s[n-St[x].len-1]!=s[n])x=St[x].fail;
	return x;
	}
	inline void add(int c){
	s.push_back(c-='a');
	++n;
	int cur=get_fail(last);
	if(!St[cur].next[c]){
	int now=St.size();
	St.push_back(St[cur].len+2);
	St[now].fail=St[get_fail(St[cur].fail)].next[c];
	/*不用擔心會找到空節點，由證明的過程可知*/
	St[cur].next[c]=now;
	St[now].num=St[St[now].fail].num+1;
	}
	last=St[cur].next[c];
	++St[last].cnt;
	}
	inline void count(){/*cnt必須要在構造完後呼叫count()去計算*/
	std::vector<node>::reverse_iterator i=St.rbegin();
	for(;i!=St.rend();++i){
	St[i->fail].cnt+=i->cnt;
	}
	}
	inline int size(){/*傳回其不同的回文子串個數*/
	return St.size()-2;
	}
	};
	#endif

view raw Palindromic Tree.cpp hosted with ❤ by GitHub

[ Heavy-Light Decomposition ] 樹鏈剖分模板+LCA

在某些狀況下，我們希望對樹上a點到b點的路徑進行修改或查詢的動作，而且這些動作用一般的樹壓平RMQ或是離線處理無法完成的，這個時候我們可以利用特殊的方法將樹拆成長鏈。
長鏈可以搭配良好的資料結構。只要找出樹上所有長鏈，每條長鏈套用線段樹、BIT、Sparse Table、BST、Heap，就能降低時間複雜度。

可以利用簡單的做法將dfs序切成許多長鏈:

1. 首先先設定其中一點為根，算出每個子樹有多少節點
2. 樹上每個節點各自連向最大的子樹，相連的部分會自然形成鏈

通常步驟二會進行一次dfs，每次朝著最大的子樹走，每個節點記錄其時間戳記及所在的鏈(通常是記鏈首)，這樣就可以將樹壓平的結果切成許多鏈，使得樹鏈剖分同時可以做到樹壓平的操作。

我們可以進行一個簡單的證明:

• 假設節點數為V，由根往任意葉節點走，一旦遇到新鏈，新鏈的子樹必小於等於原鏈子樹，剩下的節點數量不到一半，所以沿途最多遇到 logV 條鏈。

因為一條路徑藉由 LCA 拆成兩段，沿途最多遇到 2*logV 條鏈，所以在任意兩點a,b間最多只需要進行2*logV次的區間操作，若是套線段樹之類的結構可以做到$\ord{log^2V}$
這些操作都可以在求LCA的過程中完成:

	#include<vector>
	#define MAXN 100005
	typedef std::vector<int >::iterator VIT;
	int siz[MAXN],max_son[MAXN],pa[MAXN],dep[MAXN];
	/*節點大小、節點大小最大的孩子、父母節點、深度*/
	int link_top[MAXN],link[MAXN],cnt;
	/*每個點所在鏈的鏈頭、樹鏈剖分的DFS序、時間戳*/
	std::vector<int >G[MAXN];/*用vector存樹*/
	void find_max_son(int x){
	siz[x]=1;
	max_son[x]=-1;
	for(VIT i=G[x].begin();i!=G[x].end();++i){
	if(*i==pa[x])continue;
	pa[*i]=x;
	dep[*i]=dep[x]+1;
	find_max_son(*i);
	if(max_son[x]==-1||siz[*i]>siz[max_son[x]])max_son[x]=*i;
	siz[x]+=siz[*i];
	}
	}
	void build_link(int x,int top){
	link[x]=++cnt;/*記錄x點的時間戳*/
	link_top[x]=top;
	if(max_son[x]==-1)return;
	build_link(max_son[x],top);/*優先走訪最大孩子*/
	for(VIT i=G[x].begin();i!=G[x].end();++i){
	if(*i==max_son[x]||*i==pa[x])continue;
	build_link(*i,*i);
	}
	}
	inline int find_lca(int a,int b){
	/*求LCA，可以在過程中對區間進行處理*/
	int ta=link_top[a],tb=link_top[b];
	while(ta!=tb){
	if(dep[ta]<dep[tb]){
	std::swap(ta,tb);
	std::swap(a,b);
	}
	//這裡可以對a所在的鏈做區間處理
	//區間為(link[ta],link[a])
	ta=link_top[a=pa[ta]];
	}
	/*最後a,b會在同一條鏈，若a!=b還要在進行一次區間處理*/
	return dep[a]<dep[b]?a:b;
	}

view raw Heavy-Light Decomposition.cpp hosted with ❤ by GitHub

若樹鏈剖分的目的只是為了求LCA，build_link()及find_lca()有更好的做法:

	inline void build_link(int root,int n){
	pa[root]=-1;
	for(int i=1;i<=n;++i){/*假設節點編號是1~n*/
	if(~pa[i]&&max_son[pa[i]]==i)continue;
	for(int j=i;j!=-1;j=max_son[j])link_top[j]=i;
	}
	}
	inline int find_lca(int a,int b){
	while(link_top[a]!=link_top[b]){
	if(dep[link_top[a]]>dep[link_top[b]])a=pa[link_top[a]];
	else b=pa[link_top[b]];
	}
	return dep[a]<dep[b]?a:b;
	}

view raw Heavy-Light Decomposition LCA.cpp hosted with ❤ by GitHub

這就只需進行一次DFS

[ Suffix Array Longest Common Prefix (LCP) ] 後綴數組之最長共同前綴(高度數組)

不解釋了，可以在$\ord N$時間內得到，N是字串長度
以下提供模板:

	/*
	h:高度數組
	sa:後綴數組
	rank:排名
	*/
	inline void suffix_array_lcp(const char *s,int len,int *h,int *sa,int *rank){
	for(int i=0;i<len;++i)rank[sa[i]]=i;
	for(int i=0,k=0;i<len;++i){
	if(rank[i]==0)continue;
	if(k)--k;
	while(s[i+k]==s[sa[rank[i]-1]+k])++k;
	h[rank[i]]=k;
	}
	h[0]=0;
	}

view raw LCP.cpp hosted with ❤ by GitHub

[ Suffix Array SA-IS Algorithm ] 後綴數組線性SA-IS算法

SA-IS是我目前看過最快的線性後綴數組演算法，但是做為競賽用途而進行簡化後他的效率在某些硬體上會比DC3慢，不過記憶體使用量是DC3的1/3~1/2倍，而最短的實現code也比DC3短很多，因此我認為這是十分優秀的算法

因為某些SA-IS的實作方式會利用傳入的陣列s的空間進行計算，因此傳入的陣列s必須是int範圍，而傳入的字串的尾端字元必須是最小的而且在之前完全沒有出現過，因此必須給字串先做以下處理:

假設要傳入字串 mississippi
先在尾端加入字元'#': mississippi#
算法結束後陣列sa為:11 10 7 4 1 0 9 8 6 3 5 2
將第一個數字11去除後剩下的數字即為字串mississippi的後綴數組:
10 7 4 1 0 9 8 6 3 5 2
關於SA-IS算法的論文請看這裡:
Two Efficient Algorithms for Linear Time Suffix Array Construction

關於SA-IS算法的教學(專利申請文)請看這裡:
后缀数组构造方法  CN 102081673 A

SA-IS的教學投影片:
https://www.cs.helsinki.fi/u/tpkarkka/opetus/11s/spa/lecture11.pdf

SA-IS中文教學:
https://riteme.github.io/blog/2016-6-19/sais.html

如果想對Suffix Array算法進行測試請使用這個Online Judge:
http://www.spoj.com/problems/SARRAY/

以下將會提供一些SA-IS的實做模板

1.台大黑暗code界的黑暗codebook:

	#include<string.h>
	#define MAGIC(XD){\
	memset(sa,0,sizeof(int)*n);\
	memcpy(x,c,sizeof(int)*z);\
	XD\
	memcpy(x+1,c,sizeof(int)*(z-1));\
	for(i=0;i<n;++i){\
	if(sa[i]&&!t[sa[i]-1])sa[x[s[sa[i]-1]]++]=sa[i]-1;\
	}\
	memcpy(x,c,sizeof(int)*z);\
	for(i=n-1;i>=0;--i){\
	if(sa[i]&&t[sa[i]-1])sa[--x[s[sa[i]-1]]]=sa[i]-1;\
	}\
	}
	void sais(int *s,int *sa,int *p,bool *t,int *c,int n,int z){
	bool neq=t[n-1]=1;
	int nn=0,nmxz=-1,*q=sa+n,*ns=s+n,*x=c+z,lst=-1,i,j;
	memset(c,0,sizeof(int)*z);
	for(i=0;i<n;++i)++c[s[i]];
	for(i=0;i<z-1;++i)c[i+1]+=c[i];
	for(i=n-2;i>=0;i--)t[i]=(s[i]==s[i+1]?t[i+1]:s[i]<s[i+1]);
	MAGIC(
	for(i=1;i<n;++i){
	if(t[i]&&!t[i-1])sa[--x[s[i]]]=p[q[i]=nn++]=i;
	}
	);
	for(i=0;i<n;++i)if((j=sa[i])&&t[j]&&!t[j-1]){
	neq=lst<0||memcmp(s+j,s+lst,(p[q[j]+1]-j)*sizeof(int));
	ns[q[lst=j]]=nmxz+=neq;
	}
	if(nmxz==nn-1)for(i=0;i<nn;++i)q[ns[i]]=i;
	else sais(ns,q,p+nn,t+n,c+z,nn,nmxz+1);
	MAGIC(
	for(i=nn-1;i>=0;--i)sa[--x[s[p[q[i]]]]]=p[q[i]];
	);
	}
	#undef MAGIC
	/*****************這些是需要用到的陣列大小**************/
	static const int MXN=10000000;
	int s[MXN*2+5],sa[MXN*2+5],p[MXN+5],c[MXN*2+5];
	bool t[MXN*2+5];

view raw SA-IS black magic.cpp hosted with ❤ by GitHub

其實這段code本來是被壓縮得更短，而且用到非常多的記憶體，不過經過卦長的改良後成功減少記憶體的使用並將他排版成正常人能看得懂的樣子
而這裡的MXN則是字串的最長長度(假設字元數<字串長度)

2.卦長自行實作的模板(記憶體用量少):

	#define pushS(x) sa[--bkt[s[x]]] = x
	#define pushL(x) sa[bkt[s[x]]++] = x
	#define induce_sort(v){\
	fill_n(sa,n,0);\
	copy_n(_bkt,A,bkt);\
	for(i=n1-1;~i;--i)pushS(v[i]);\
	copy_n(_bkt,A-1,bkt+1);\
	for(i=0;i<n;++i)if(sa[i]&&!t[sa[i]-1])pushL(sa[i]-1);\
	copy_n(_bkt,A,bkt);\
	for(i=n-1;~i;--i)if(sa[i]&&t[sa[i]-1])pushS(sa[i]-1);\
	}
	template<typename T>
	void sais(const T s,int n,int *sa,int *_bkt,int *p,bool *t,int A){
	int *rnk=p+n,*s1=p+n/2,*bkt=_bkt+A;
	int n1=0,i,j,x=t[n-1]=1,y=rnk[0]=-1,cnt=-1;
	for(i=n-2;~i;--i)t[i]=(s[i]==s[i+1]?t[i+1]:s[i]<s[i+1]);
	for(i=1;i<n;++i)rnk[i]=t[i]&&!t[i-1]?(p[n1]=i,n1++):-1;
	fill_n(_bkt,A,0);
	for(i=0;i<n;++i)++_bkt[s[i]];
	for(i=1;i<A;++i)_bkt[i]+=_bkt[i-1];
	induce_sort(p);
	for(i=0;i<n;++i)if(~(x=rnk[sa[i]]))
	j=y<0||memcmp(s+p[x],s+p[y],(p[x+1]-p[x])*sizeof(s[0]))
	,s1[y=x]=cnt+=j;
	if(cnt+1<n1)sais(s1,n1,sa,bkt,rnk,t+n,cnt+1);
	else for(i=0;i<n1;++i)sa[s1[i]]=i;
	for(i=0;i<n1;++i)s1[i]=p[sa[i]];
	induce_sort(s1);
	}
	/*****************這些是需要用到的陣列大小**************/
	const int MAXN=10000005,MAXA='z'+1;
	int sa[MAXN],bkt[MAXN+MAXA],p[MAXN*2];
	bool t[MAXN*2];
	char s[MAXN];

view raw SA-IS.cpp hosted with ❤ by GitHub

為了簡化實作方法及減少記憶體的使用，因此將計算後剩下的空間進行重複利用，壓縮後只需要這些記憶空間，而傳入的陣列s可以直接傳入char陣列，因此對使用者來說是非常方便的一份code
這裡的MXN則是字串的最長長度(假設字元數<字串長度)

3.論文提供的實作code:

	/*SA-IS Algorithm*/
	const unsigned char mask[]={0x80,0x40,0x20,0x10,0x08,0x04,0x02,0x01};
	#define tget(i) ( (t[(i)/8]&mask[(i)%8])?1:0)
	#define tset(i, b) t[(i)/8]=(b)?(mask[(i)%8]|t[(i)/8]):((~mask[(i)%8])&t[(i)/8])
	#define chr(i) (cs==sizeof(int)?((int*)s)[i]:((unsigned char *)s)[i])
	#define isLMS(i) (i>0&&tget(i)&&!tget(i-1))
	// find the start or end of each bucket
	inline void getBuckets(unsigned char *s,int *bkt,int n,int K,int cs,bool end){
	int i,sum=0;
	// clear all buckets
	for(i=0;i<=K;++i)bkt[i]=0;
	// compute the size of each bucket
	for(i=0;i<n;++i)++bkt[chr(i)];
	for(i=0;i<=K;++i){
	sum+=bkt[i];
	bkt[i]=end?sum:sum-bkt[i];
	}
	}
	// compute SAl
	inline void induceSAl(unsigned char *t,int *SA,unsigned char *s,int *bkt,int n,int K,int cs,bool end){
	int i,j;
	// find starts of buckets
	getBuckets(s,bkt,n,K,cs,end);
	for(i=0;i<n;++i){
	j=SA[i]-1;
	if(j>=0&&!tget(j))SA[bkt[chr(j)]++]=j;
	}
	}
	// compute SAs
	inline void induceSAs(unsigned char *t,int *SA,unsigned char *s,int *bkt,int n,int K,int cs,bool end){
	int i,j;
	// find ends of buckets
	getBuckets(s,bkt,n,K,cs,end);
	for(i=n-1;i>=0;--i){
	j=SA[i]-1;
	if(j>=0&&tget(j))SA[--bkt[chr(j)]]=j;
	}
	}
	// find the suffix array SA of s[0..n-1] in {1..K}?n
	// require s[n-1]=0 (the sentinel!), n>=2
	// use a working space (excluding s and SA) of
	// at most 2.25n+O(1) for a constant alphabet
	void SA_IS(unsigned char *s,int *SA,int n,int K,short cs){
	// LS-type array in bits
	unsigned char *t=(unsigned char *)malloc(n/8+1);
	int i,j;
	// classify the type of each character
	// the sentinel must be in s1, important!!!
	tset(n-2,0);tset(n-1,1);
	for(i=n-3;i>=0;--i)
	tset(i,(chr(i)<chr(i+1)||(chr(i)==chr(i+1)&&tget(i+1)==1))?1:0);
	// stage 1: reduce the problem by at least 1/2
	// sort all the S-substrings
	// bucket array
	int *bkt=(int *)malloc(sizeof(int)*(K+1));
	// find ends of buckets
	getBuckets(s,bkt,n,K,cs,true);
	for(i=0;i<n;++i)SA[i]=-1;
	for(i=1;i<n;++i)if(isLMS(i))SA[--bkt[chr(i)]]=i;
	induceSAl(t,SA,s,bkt,n,K,cs,false);
	induceSAs(t,SA,s,bkt,n,K,cs,true);
	free(bkt);
	// compact all the sorted substrings into
	// the first n1 items of SA
	// 2*n1 must be not larger than n (proveable)
	int n1=0;
	for(i=0;i<n;++i)if(isLMS(SA[i]))SA[n1++]=SA[i];
	// find the lexicographic names of substrings
	// init the name array buffer
	for(i=n1;i<n;++i)SA[i]=-1;
	int name=0,prev=-1;
	for(i=0;i<n1;++i){
	int pos=SA[i];bool diff=false;
	for(int d=0;d<n;++d)
	if(prev==-1||chr(pos+d)!=chr(prev+d)||tget(pos+d)!=tget(prev+d)){
	diff=true;break;
	}else if(d>0&&(isLMS(pos+d)||isLMS(prev+d)))break;
	if(diff){
	++name;prev=pos;
	}
	pos=(pos%2==0)?pos/2:(pos-1)/2;
	SA[n1+pos]=name-1;
	}
	for(i=n-1,j=n-1;i>=n1;--i)if(SA[i]>=0)SA[j--]=SA[i];
	// stage 2: solve the reduced problem
	// recurse if names are not yet unique
	int *SA1=SA,*s1=SA+n-n1;
	if(name<n1)SA_IS((unsigned char*)s1,SA1,n1,name-1,sizeof(int));
	else // generate the suffix array of s1 directly
	for(i=0;i<n1;++i)SA1[s1[i]]=i;
	// stage 3: induce the result for
	// the original problem
	// bucket array
	bkt=(int *)malloc(sizeof(int)*(K+1));
	// put all the LMS characters into their buckets
	// find ends of buckets
	getBuckets(s,bkt,n,K,cs,true);
	for(i=1,j=0;i<n;++i)if(isLMS(i))s1[j++]=i; // get p1
	// get index in s
	for(i=0;i<n1;++i)SA1[i]=s1[SA1[i]];
	// init SA[n1..n-1]
	for(i=n1;i<n;++i)SA[i]=-1;
	for(i=n1-1;i>=0;--i){
	j=SA[i];SA[i]=-1;
	SA[--bkt[chr(j)]]=j;
	}
	induceSAl(t,SA,s,bkt,n,K,cs,false);
	induceSAs(t,SA,s,bkt,n,K,cs,true);
	free(bkt);free(t);
	}

view raw SA-IS paper code.c hosted with ❤ by GitHub

這是其中一個比DC3快的code，而且記憶體使用量是最少的，但是長度很長就是了，不適合在比賽時使用

4.超快記憶體使用超少的模板庫code:
https://gist.github.com/jacky860226/1d33adad858eef71bfe18120d8d69e6d#file-sa-is-very-fast-cpp
因為長度太長所以就直接貼上網址了，沒有人會在比賽時寫這種東西

[ Suffix Array DC3 algoruthm ] 後綴數組線性DC3演算法

這邊提供後綴數組DC3算法的模板，他是一個能在線性時間$\ord N$內求得後綴數組的演算法，注意這邊傳入的陣列s及sa長度必須是字串長的3倍以上，因為DC3會用到大量的記憶體，而且會利用傳入的陣列s的空間進行計算，因此傳入的陣列s必須是int範圍，而傳入的字串的尾端字元必須是最小的而且在之前完全沒有出現過，因此必須給字串先做以下處理:

假設要傳入字串 mississippi
先在尾端加入字元'#': mississippi#
算法結束後陣列sa為:11 10 7 4 1 0 9 8 6 3 5 2
將第一個數字11去除後剩下的數字即為字串mississippi的後綴數組:
10 7 4 1 0 9 8 6 3 5 2

因為是遞迴的關係，所以常數很大，但在大數據的時候比倍增法快
關於DC3演算法的介紹請看這兩篇文章:
http://ching.im/post/algorithm/difference-cover-modulo-algorithm
《后缀数组——处理字符串的有力工具》

如果想對Suffix Array算法進行測試請使用這個Online Judge:
http://www.spoj.com/problems/SARRAY/

在傳入的參數中最後一個A為最大的字元+1通常設為'z'+1

以下提供模板(其實code挺短的):

	#define F(x) (x)/3+((x)%3==1?0:tb)
	#define G(x) (x)<tb?(x)*3+1:((x)-tb)*3+2
	int wa[3*maxn+5],wb[3*maxn+5],wv[3*maxn+5],c[maxn+5];
	inline bool c0(int *s,int a,int b){
	return s[a]==s[b]&&s[a+1]==s[b+1]&&s[a+2]==s[b+2];
	}
	inline bool c12(int k,int *s,int a,int b){
	if(k==2)return s[a]<s[b]||s[a]==s[b]&&c12(1,s,a+1,b+1);
	else return s[a]<s[b]||s[a]==s[b]&&wv[a+1]<wv[b+1];
	}
	inline void radix_sort(int *s,int *a,int *b,int len,int A){
	static int i;
	for(i=0;i<len;++i)wv[i]=s[a[i]];
	for(i=0;i<A;++i)c[i]=0;
	for(i=0;i<len;++i)++c[wv[i]];
	for(i=1;i<A;++i)c[i]+=c[i-1];
	for(i=len-1;i>=0;--i)b[--c[wv[i]]]=a[i];
	}
	void dc3(int *s,int *sa,int len,int A){
	int i,j,*san=sa+len,ta=0,tb=(len+1)/3,tbc=0,p;
	int *rn=s+len;
	s[len]=s[len+1]=0;
	for(i=0;i<len;++i){
	if(i%3!=0)wa[tbc++]=i;
	}
	radix_sort(s+2,wa,wb,tbc,A);
	radix_sort(s+1,wb,wa,tbc,A);
	radix_sort(s,wa,wb,tbc,A);
	for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;++i){
	rn[F(wb[i])]=c0(s,wb[i-1],wb[i])?p-1:p++;
	}
	if(p<tbc)dc3(rn,san,tbc,p);
	else for(i=0;i<tbc;i++)san[rn[i]]=i;
	for(i=0;i<tbc;++i){
	if(san[i]<tb)wb[ta++]=san[i]*3;
	}
	if(len%3==1)wb[ta++]=len-1;
	radix_sort(s,wb,wa,ta,A);
	for(i=0;i<tbc;++i)wv[wb[i]=G(san[i])]=i;
	for(i=0,j=0,p=0;i<ta&&j<tbc;++p){
	sa[p]=c12(wb[j]%3,s,wa[i],wb[j])?wa[i++]:wb[j++];
	}
	for(;i<ta;++p)sa[p]=wa[i++];
	for(;j<tbc;++p)sa[p]=wb[j++];
	}
	#undef F
	#undef G
	/*
	DC3的輸入陣列s必須能儲存0~len的數值，所以設為int
	s及sa陣列的大小必須大於3*maxn
	*/

view raw DC3.cpp hosted with ❤ by GitHub

[ Suffix Array Prefix doubling algorithms ] 後綴數組倍增算法

後綴數組(又稱尾碼陣列)是一個十分強大的字串處理武器，大部分的問題都可以用它來解決，它可以幾乎做到所有後綴樹(Suffix Tree)能做到的事，所以這邊就不介紹後綴樹了

因為後綴數組可以由後綴樹進行遍歷轉換而來，而構造後綴樹僅需花費線性的時間，所以構造後綴數組的時間可為線性$\ord N$，但是後綴數組本身就是為了減少構造後綴樹的空間與代碼量而被發明出來的，直接由後綴樹轉換是沒意義的
但是仍然有其他線性構造後綴數組的方法，像是DC3、SA-IS等會在下一篇介紹，這次要講的是比較簡單常用的方法-----倍增法
關於後綴數組的使用說明可以參考《后缀数组——处理字符串的有力工具》
關於倍增法的說明可以參考 演算法筆記-SuffixArray 的部分
這邊提供$\ord{N*logN*logN}$及$\ord{NlogN}$的模板

$\ord{N*logN*logN}$:

	#include<algorithm>
	struct CMP{
	int len,k,*rank,a,b;
	inline bool operator()(int i,int j){
	if(rank[i]!=rank[j])return rank[i]<rank[j];
	a=(i+=k)<len?rank[i]:-1;
	b=(j+=k)<len?rank[j]:-1;
	return a<b;
	}
	};
	inline void suffix_array(const char *s,int len,int *sa,int *rank,int *tmp){
	for(int i=0;i<len;++i){
	sa[i]=i;
	rank[i]=s[i];
	}
	CMP cmp={len,1};
	for(;;cmp.k<<=1){
	cmp.rank=rank;
	std::sort(sa,sa+len,cmp);
	tmp[sa[0]]=0;
	for(int i=1;i<len;++i)tmp[sa[i]]=tmp[sa[i-1]]+cmp(sa[i-1],sa[i]);
	if(tmp[sa[len-1]]==len-1)break;
	std::swap(rank,tmp);
	}
	}

view raw doubling algorithm NlogNlogN.cpp hosted with ❤ by GitHub

$\ord{NlogN}$:

	#include<algorithm>
	#define radix_sort(x,y){\
	for(i=0;i<A;++i)c[i]=0;\
	for(i=0;i<len;++i)c[x[y[i]]]++;\
	for(i=1;i<A;++i)c[i]+=c[i-1];\
	for(i=len-1;i>=0;--i)sa[--c[x[y[i]]]]=y[i];\
	}
	#define sac(r,a,b) r[a]!=r[b]||a+k>=len||r[a+k]!=r[b+k]
	inline void suffix_array(const char *s,int len,int *sa,int *rank,int *tmp,int *c){
	int A='z'+1,i,k,id=0;
	for(i=0;i<len;++i)rank[tmp[i]=i]=s[i];
	radix_sort(rank,tmp);
	for(k=1;id<len-1;k<<=1){
	for(id=0,i=len-k;i<len;++i)tmp[id++]=i;
	for(i=0;i<len;++i)if(sa[i]>=k)tmp[id++]=sa[i]-k;
	radix_sort(rank,tmp);
	std::swap(rank,tmp);
	for(rank[sa[0]]=id=0,i=1;i<len;++i)
	rank[sa[i]]=id+=sac(tmp,sa[i-1],sa[i]);
	A=id+1;
	}
	}
	#undef radix_sort
	#undef sac

view raw doubling algorithm NlogN.cpp hosted with ❤ by GitHub

注意此方法必須要在字元集大小為常數的情況下有效，否則必須離散化

所需的陣列長度只需要與字串陣列相同即可
當然$N*logN*logN$的做法會比較值觀，$NlogN$的方法則是利用radix_sort進行的，radix_sort本來在倍增的時候要先排序第一關鍵字跟第二關鍵字，但是第二關鍵字排序的結果可以用已經求好的SA直接求出來
對radix_sort還不了解的人請看這個網頁:https://www.cs.usfca.edu/~galles/visualization/RadixSort.html
如果想對Suffix Array算法進行測試請使用這個Online Judge:
http://www.spoj.com/problems/SARRAY/

英文字母大小寫轉換特殊做法

假設有一個題目是這樣的:
        給定一串英文字母，請將大寫的部分轉成小寫，小寫的部分轉成大寫並輸出

一般我們會用if或是三元運算子做判斷，但是這太麻煩了
經過觀察發現，摁合一個小寫字母ascii與其對應的大寫字母ascii相差皆為32，
而其二進位編碼剛好允許透過 xor 32的方式進行轉換，但是32這個數字不好記，又可以發現ascii 32='空白'，而以下的code就可以將大小寫互換:

1 2 3 4 5 6 7 8

#include<stdio.h> char s[1000005]; int main(){ scanf("%s",s); for(int i=0;s[i];++i)putchar(s[i]^' ');puts(""); for(int i=0;s[i];++i)putchar(s[i]^32);puts(""); return 0; }

兩種寫法會有相同的效果

[ Persistent Binary Index Tree ] 持久化BIT

(二分索引樹)BIT是一種容易實現的資料結構，但在持久化的應用上是非常困難的，但我們可以犧牲些許的複雜度來簡化其實現方式

一般區間和的BIT長這樣:

	int n;
	int tree[100005];
	#define lowbit(x) x&(-x)
	inline void modify(int x,int d){
	for(;x<=n;x+=lowbit(x))tree[x]+=d;
	}
	inline int query(int x){
	int ans=0;
	for(;x;x-=lowbit(x))ans+=tree[x];
	return ans;
	}

view raw BIT.cpp hosted with ❤ by GitHub

這是持久化的版本(也可以用pair來實作):

	#include<vector>
	#include<algorithm>
	#define lowbit(x) x&(-x)
	struct P{
	int data,id;
	P(int d=0,int i=0):data(d),id(i){}
	inline friend bool operator<(const P &a,const P &b){
	return a.id<b.id;
	}
	};
	int n,now;
	std::vector<P >tree[100005];
	inline void init(){
	now=0;
	for(int i=1;i<=n;++i)tree[i].clear(),tree[i].push_back(P());
	}
	inline void modify(int x,int d){
	for(;x<=n;x+=lowbit(x))tree[x].push_back(P(tree[x].back().data+d,now));
	++now;
	}
	inline int query(int x,int id){/*查詢第id次操作的區間和，id從0開始計算*/
	int ans=0;
	std::vector<P >::iterator a;
	for(;x;x-=lowbit(x)){
	a=std::upper_bound(tree[x].begin(),tree[x].end(),P(0,id))-1;
	ans+=a->data;
	}
	return ans;
	}

view raw Persistent BIT.cpp hosted with ❤ by GitHub

可以知道其更新的複雜度為$\ord{logN}$，查詢複雜度為$\ord{logN*logQ}$，Q是查詢次數

[ Suffix Automaton ] 後綴自動機

定義字串S的後綴自動機是一種自動機可以接受S的所有後綴，其狀態最多為2*strlen(s)-1個，雖然理論複雜，但其構造方式十分簡單，且它可以用來處理一些後綴數組或事後綴樹的問題，因此是一個十分強大的資料結構。
後綴自動機大致的構造方式是將字串擔一字元按順序進行插入，插入的均攤時間為$\ord 1$，關於詳細的情形可以在 http://blog.sina.com.cn/s/blog_70811e1a01014dkz.html 這個網站中得到
而關於複雜度及詳細的證明請參考陈立杰的簡報
關於狀態數為線性的證明在19~21頁，有詳細的圖文解說
這個網站包含其一些延伸資料結構http://maskray.me/blog/2013-05-10-suffix-automaton
為了減少記憶體的使用，以vector代替指標維護狀態
以下提供模板:

	#ifndef SUFFIX_AUTOMATON
	#define SUFFIX_AUTOMATON
	#include<vector>
	template<char L='a',char R='z'>
	struct suffix_automaton{
	struct node{
	int ch[R-L+1],fa,len;
	node(int l=0):fa(0),len(l){
	for(int i=0;i<=R-L;++i)ch[i]=0;
	}
	};
	std::vector<node >S;
	int tail,u,v,np,vp;
	suffix_automaton():S(2),tail(1){S[0].fa=0;}
	inline void add(int c){
	c-=L;
	u=tail;
	S.push_back(node(S[u].len+1));
	tail=np=S.size()-1;
	for(;u&&!S[u].ch[c];u=S[u].fa)S[u].ch[c]=np;
	if(!u)S[np].fa=1;
	else{
	v=S[u].ch[c];
	if(S[v].len==S[u].len+1)S[np].fa=v;
	else{
	S.push_back(S[v]);
	vp=S.size()-1;
	S[vp].len=S[u].len+1;
	S[v].fa=S[np].fa=vp;
	for(;u&&S[u].ch[c]==v;u=S[u].fa)S[u].ch[c]=vp;
	}
	}
	}
	inline int size(){return S.size()-2;}
	};
	#endif

view raw Suffix Automaton.cpp hosted with ❤ by GitHub

[ Aho-Corasick Automaton ] AC自動機

Aho-Corasick Automaton，該算法在1975年產生於貝爾實驗室，是著名的多模匹配算法之一
在了解AC自動機之前，必須先學會字典樹Trie的使用，簡單的來說，他是將KMP的fail函數應用在Trie做為失敗邊。
有$n$個字串$B_0,B_1,...,B_{n-1}$ 要匹配字串$A$，則一般用KMP的算法會得到$\ord{\sum_{i=0}^{n-1}(|A|+|B_i|)}$，使用AC自動機的演算法其複雜度為$\ord{|A|+\sum_{i=0}^{n-1}|B_i|}$，必須用DP

關於AC自動機的介紹請看這裡

以下為AC自動機的模板

	#ifndef SUNMOON_AHO_CORASICK_AUTOMATON
	#define SUNMOON_AHO_CORASICK_AUTOMATON
	#include<queue>
	#include<vector>
	template<char L='a',char R='z'>
	class ac_automaton{
	private:
	struct joe{
	int next[R-L+1],fail,efl,ed,cnt_dp,vis;
	joe():ed(0),cnt_dp(0),vis(0){
	for(int i=0;i<=R-L;++i)next[i]=0;
	}
	};
	public:
	std::vector<joe> S;
	std::vector<int> q;
	int qs,qe,vt;
	ac_automaton():S(1),qs(0),qe(0),vt(0){}
	inline void clear(){
	q.clear();
	S.resize(1);
	for(int i=0;i<=R-L;++i)S[0].next[i]=0;
	S[0].cnt_dp=S[0].vis=qs=qe=vt=0;
	}
	inline void insert(const char *s){
	int o=0;
	for(int i=0,id;s[i];++i){
	id=s[i]-L;
	if(!S[o].next[id]){
	S.push_back(joe());
	S[o].next[id]=S.size()-1;
	}
	o=S[o].next[id];
	}
	++S[o].ed;
	}
	inline void build_fail(){
	S[0].fail=S[0].efl=-1;
	q.clear();
	q.push_back(0);
	++qe;
	while(qs!=qe){
	int pa=q[qs++],id,t;
	for(int i=0;i<=R-L;++i){
	t=S[pa].next[i];
	if(!t)continue;
	id=S[pa].fail;
	while(~id&&!S[id].next[i])id=S[id].fail;
	S[t].fail=~id?S[id].next[i]:0;
	S[t].efl=S[S[t].fail].ed?S[t].fail:S[S[t].fail].efl;
	q.push_back(t);
	++qe;
	}
	}
	}
	/*DP出每個前綴在字串s出現的次數並傳回所有字串被s匹配成功的次數O(N+M)*/
	inline int match_0(const char *s){
	int ans=0,id,p=0,i;
	for(i=0;s[i];++i){
	id=s[i]-L;
	while(!S[p].next[id]&&p)p=S[p].fail;
	if(!S[p].next[id])continue;
	p=S[p].next[id];
	++S[p].cnt_dp;/*匹配成功則它所有後綴都可以被匹配(DP計算)*/
	}
	for(i=qe-1;i>=0;--i){
	ans+=S[q[i]].cnt_dp*S[q[i]].ed;
	if(~S[q[i]].fail)S[S[q[i]].fail].cnt_dp+=S[q[i]].cnt_dp;
	}
	return ans;
	}
	/*多串匹配走efl邊並傳回所有字串被s匹配成功的次數O(N*M^1.5)*/
	inline int match_1(const char *s)const{
	int ans=0,id,p=0,t;
	for(int i=0;s[i];++i){
	id=s[i]-L;
	while(!S[p].next[id]&&p)p=S[p].fail;
	if(!S[p].next[id])continue;
	p=S[p].next[id];
	if(S[p].ed)ans+=S[p].ed;
	for(t=S[p].efl;~t;t=S[t].efl){
	ans+=S[t].ed;/*因為都走efl邊所以保證匹配成功*/
	}
	}
	return ans;
	}
	/*枚舉(s的子字串∩A)的所有相異字串各恰一次並傳回次數O(N*M^(1/3))*/
	inline int match_2(const char *s){
	int ans=0,id,p=0,t;
	++vt;
	/*把戳記vt+=1，只要vt沒溢位，所有S[p].vis==vt就會變成false
	這種利用vt的方法可以O(1)歸零vis陣列*/
	for(int i=0;s[i];++i){
	id=s[i]-L;
	while(!S[p].next[id]&&p)p=S[p].fail;
	if(!S[p].next[id])continue;
	p=S[p].next[id];
	if(S[p].ed&&S[p].vis!=vt){
	S[p].vis=vt;
	ans+=S[p].ed;
	}
	for(t=S[p].efl;~t&&S[t].vis!=vt;t=S[t].efl){
	S[t].vis=vt;
	ans+=S[t].ed;/*因為都走efl邊所以保證匹配成功*/
	}
	}
	return ans;
	}
	/*把AC自動機變成真的自動機*/
	inline void evolution(){
	for(qs=1;qs!=qe;){
	int p=q[qs++];
	for(int i=0;i<=R-L;++i)
	if(S[p].next[i]==0)S[p].next[i]=S[S[p].fail].next[i];
	}
	}
	};
	#endif

view raw Aho-Corasick automaton.cpp hosted with ❤ by GitHub

[ Knuth-Morris-Pratt Algorithm ] 克努斯-莫里斯-普拉特(KMP)算法

這是一個大家很常聽到的字串匹配演算法，較Z Algorithm複雜，但更為通用，C++的string就有內建(string::find)。假設要以B匹配A，則會先建立B的部分匹配表(又叫做失配函數)，其定義如下:
     fail(i)=max(滿足B[0,k]=B[i-k,i]的所有k)
同樣的也可以在線性時間來完成
對於KMP演算法的詳細情況請參考這裡
以下提供fail function陣列及匹配的實作

	/*產生fail function*/
	inline void kmp_fail(char *s,int len,int *fail){
	int id=-1;
	fail[0]=-1;
	for(int i=1;i<len;++i){
	while(~id&&s[id+1]!=s[i])id=fail[id];
	if(s[id+1]==s[i])++id;
	fail[i]=id;
	}
	}
	/*以字串B匹配字串A，傳回匹配成功的數量*/
	inline int kmp_match(char *A,int lenA,char *B,int lenB,int *fail){
	int id=-1,ans=0;
	for(int i=0;i<lenA;++i){
	while(~id&&B[id+1]!=A[i])id=fail[id];
	if(B[id+1]==A[i])++id;
	if(id==lenB-1){/*匹配成功*/
	++ans;
	id=fail[id];
	}
	}
	return ans;
	}

view raw KMP.cpp hosted with ❤ by GitHub

[ Manacher's algorithm ] Linear time Longest palindromic substring 線性最長回文子串

Manacher演算法是Z algorithm的變種，可以說是雙向的Z algorithm，複雜度也是$\ord N$
其Z function定義如下:
     Z(i)=以位置i為中心的最長回文半徑
但是這樣只能處理回文長度是奇數的子串，因此必須幫字串做一些修改

假設原來的字串是: "asdsasdsa"
修改後變成: "@#a#s#d#s#a#s#d#s#a#"
其中'@'及'#'為不同字元且皆未在原字串中出現過
其Z function陣列為: 0 1 2 1 2 1 6 1 2 1 10 1 2 1 6 1 2 1 2 1
則 最大的數字-1 (10-1=9)就是我們要的答案

這裡提供Manacher algorithm產生Z function陣列的實作:

	/*
	原字串: asdsasdsa
	先把字串變成這樣: @#a#s#d#s#a#s#d#s#a#
	*/
	inline void manacher(char *s,int len,int *z){
	int l=0,r=0;
	for(int i=1;i<len;++i){
	z[i]=r>i?min(z[2*l-i],r-i):1;
	while(s[i+z[i]]==s[i-z[i]])++z[i];
	if(z[i]+i>r)r=z[i]+i,l=i;
	}
	}

view raw Manacher.cpp hosted with ❤ by GitHub

[ Z algorithm ] Linear-time pattern matching 線性字串匹配 Z演算法

定義一個字串S的Z function:
    Z(i)=0 如果i=0 or S[i]!=S[0]     否則
    Z(i)=max(所有的k滿足 S[0,k-1]=S[i,i+k-1])
從Z function的定義，假設要在字串A裡面匹配字串B
可以先建構字串S=B+(A和B都沒有的字元)+A的Z function陣列
若Z[i]=lenB則表示在A[i-lenB]有一個完全匹配
這裡提供兩個可以$\ord{N}$時間求出Z function陣列的方法:

	inline void z_alg1(char *s,int len,int *z){
	int l=0,r=0;
	z[0]=len;
	for(int i=1;i<len;++i){
	z[i]=r>i?min(r-i+1,z[z[l]-(r-i+1)]):0;
	while(i+z[i]<len&&s[z[i]]==s[i+z[i]])++z[i];
	if(i+z[i]-1>r)r=i+z[i]-1,l=i;
	}
	}
	inline void z_alg2(char *s,int len,int *z){
	int l=0,r=0;
	z[0]=len;
	for(int i=1;i<len;++i){
	z[i]=i>r?0:(i-l+z[i-l]<z[l]?z[i-l]:r-i+1);
	while(i+z[i]<len&&s[i+z[i]]==s[z[i]])++z[i];
	if(i+z[i]-1>r)r=i+z[i]-1,l=i;
	}
	}

view raw Z.cpp hosted with ❤ by GitHub

對一個特殊構建的數據strlen(s)=10000000其效率平均為
z_alg1:   0.106s
z_alg2:   0.089s
這兩種方法的想法是相同的
我們可以知道第一個方法是較為直觀好寫的，但是運算量較大且調用min函數的效率並不高
但仍在可接受的範圍

[ skew heap ] implementation 斜堆 實作

斜堆(Skew heap)也叫自適應堆(self-adjusting heap)，它是(重量)左偏樹的一個變種。
相較於(重量)左偏樹，斜堆不需記錄其高度或是節點個數(size)等附加域，其合併(merge)的方式也較為簡潔，效率也較高。和(重量)左偏樹一樣，斜堆的push、pop的時間複雜度為$\ord{log \; n}$、查詢最值為$\ord 1$。
關於複雜度分析
以下提供斜堆的實作，使用方法與STL priority_queue相似

	#include<functional>
	#ifndef SKEW_HEAP
	#define SKEW_HEAP
	template<typename T,typename _Compare=std::less<T> >
	class skew_heap{
	private:
	struct node{
	T data;
	node *l,*r;
	node(const T&d):data(d),l(0),r(0){}
	}*root;
	int _size;
	_Compare cmp;
	node *merge(node *a,node *b){
	if(!a||!b)return a?a:b;
	if(cmp(a->data,b->data))return merge(b,a);
	node *t=a->r;
	a->r=a->l;
	a->l=merge(b,t);
	return a;
	}
	void _clear(node *&o){
	if(o)_clear(o->l),_clear(o->r),delete o;
	}
	public:
	skew_heap():root(0),_size(0){}
	~skew_heap(){_clear(root);}
	inline void clear(){
	_clear(root);root=0;_size=0;
	}
	inline void join(skew_heap &o){
	root=merge(root,o.root);
	o.root=0;
	_size+=o._size;
	o._size=0;
	}
	inline void swap(skew_heap &o){
	node *t=root;
	root=o.root;
	o.root=t;
	int st=_size;
	_size=o._size;
	o._size=st;
	}
	inline void push(const T&data){
	_size++;
	root=merge(root,new node(data));
	}
	inline void pop(){
	if(_size)_size--;
	node *tmd=merge(root->l,root->r);
	delete root;
	root=tmd;
	}
	inline const T& top(){return root->data;}
	inline int size(){return _size;}
	inline bool empty(){return !_size;}
	};
	#endif

view raw skew heap.cpp hosted with ❤ by GitHub

[ basic operation of naive order statistic tree ] 樸素二元搜尋樹 名次樹的基本操作

一般的名次樹支援以下幾種操作:

1.插入
2.刪除
3.查找
4.求節點大小
5.求元素名次
6.求第k大
7.求前驅
8.求後繼

之前寫的模板已經有了前六種操作，因此我會在本模板中附上求前驅跟後繼的操作
前驅與後繼並不是名次樹才有的操作，在不記錄節點大小的情形下也能求得
關於求前驅跟後繼操作的說明可以參考 随机平衡二叉查找树Treap的分析与应用 這篇文章

以下提供模板:

	#ifndef NAIVE_BINARY_SEARCH_TREE
	#define NAIVE_BINARY_SEARCH_TREE
	template<typename T>
	class bs_tree{
	private:
	struct node{
	node *ch[2];
	int s;
	T data;
	node(const T&d):s(1),data(d){ch[0]=ch[1]=0;}
	}*root;
	inline int size(node *o){return o?o->s:0;}
	void insert(node *&o,const T&data){
	if(!o)o=new node(data);
	else o->s++,insert(o->ch[o->data<data],data);
	}
	inline void success_erase(node *&o){
	node *t=o;
	o=o->ch[0]?o->ch[0]:o->ch[1];
	delete t;
	}
	bool erase(node *&o,const T &data){
	if(!o)return 0;
	if(o->data==data){
	if(!o->ch[0]||!o->ch[1])success_erase(o);
	else{
	o->s--;
	node **t=&o->ch[1];
	for(;(*t)->ch[0];t=&(*t)->ch[0])(*t)->s--;
	o->data=(*t)->data;
	success_erase(*t);
	}return 1;
	}else if(erase(o->ch[o->data<data],data)){
	o->s--;return 1;
	}return 0;
	}
	node *preorder(node *&x,node *y,const T&data){
	if(!x)return y;
	if(x->data<data)return preorder(x->ch[1],x,data);
	else return preorder(x->ch[0],y,data);
	}
	node *successor(node *&x,node *y,const T&data){
	if(!x)return y;
	if(data<x->data)return successor(x->ch[0],x,data);
	else return successor(x->ch[1],y,data);
	}
	void clear(node *&p){
	if(p)clear(p->ch[0]),clear(p->ch[1]),delete p;
	}
	public:
	bs_tree():root(0){}
	~bs_tree(){clear(root);}
	inline void clear(){clear(root),root=0;}
	inline void insert(const T &data){
	insert(root,data);
	}
	inline bool erase(const T &data){
	return erase(root,data);
	}
	inline bool find(const T &data){
	for(node *o=root;o;)
	if(o->data==data)return 1;
	else o=o->ch[o->data<data];
	return 0;
	}
	inline const T&preorder(const T&data){
	node *o=preorder(root,0,data);
	if(o)return o->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *o=successor(root,0,data);
	if(o)return o->data;
	return data;
	}
	inline int size(){return size(root);}
	inline int rank(const T &data){
	int cnt=0;
	for(node *o=root;o;)
	if(o->data<data)cnt+=size(o->ch[0])+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T &kth(int k){
	for(node *o=root;;)
	if(k<=size(o->ch[0]))o=o->ch[0];
	else if(k==size(o->ch[0])+1)return o->data;
	else k-=size(o->ch[0])+1,o=o->ch[1];
	}
	inline const T &operator[](int k){
	return kth(k);
	}
	};
	#endif

view raw bst.cpp hosted with ❤ by GitHub

[ AA tree ] AA樹

AA樹是紅黑樹的一種變種，是Arne Andersson教授在1993年年在他的論文"Balanced search trees made simple"中介紹，設計的目的是減少紅黑樹考慮的不同情況。區別於紅黑樹的是，AA樹的紅結點只能作為右葉子。另外AA樹為實現方便，不再使用紅黑兩種顏色，而是用level標記結點，結點中的level相當於紅黑樹中結點的黑高度。

其主要透過skew和split兩種旋轉來進行平衡的，插入時只能為紅節點，故會出現不符合其規定則進行平衡操作

若想進一步了解其平衡操作請參考維基百科
在實作中，skew為右璇，split為左旋。本模板以rotate(o,0)、rotate(o,1)代表

以下提供模板:

	#ifndef ARNE_ANDERSSON_TREE
	#define ARNE_ANDERSSON_TREE
	template<typename T>
	class aa_tree{
	private:
	struct node{
	node *ch[2];
	int s,level;
	T data;
	node(const T&d):s(1),level(1),data(d){}
	node():s(0),level(0){ch[0]=ch[1]=this;}
	}*nil,*root;
	inline void rotate(node *&a,bool d){
	if(a->level==a->ch[d]->ch[d]->level){
	node *b=a;
	a=a->ch[d];
	b->ch[d]=a->ch[!d];
	a->ch[!d]=b;
	if(d)a->level++;
	a->s=b->s;
	b->s=b->ch[0]->s+b->ch[1]->s+1;
	}
	}
	inline void erase_fix(node *&o){
	if(o ->ch[0]->level<o->level-1||o->ch[1]->level<o->level-1){
	if(o->ch[1]->level>--o->level)o->ch[1]->level=o->level;
	rotate(o,0);
	if(o->ch[1]->s)rotate(o->ch[1],0);
	if(o->ch[1]->ch[1]->s)rotate(o->ch[1]->ch[1],0);
	rotate(o,1);
	if(o->ch[1]->s)rotate(o->ch[1],1);
	}
	}
	void insert(node *&o,const T&data){
	if(!o->s){
	o=new node(data);
	o->ch[0]=o->ch[1]=nil;
	}else{
	o->s++;
	insert(o->ch[o->data<data],data);
	rotate(o,0);
	rotate(o,1);
	}
	}
	bool erase(node *&o,const T&data){
	if(!o->s)return 0;
	if(o->data==data){
	if(!o->ch[0]->s||!o->ch[1]->s){
	node *t=o;
	o=o->ch[0]->s?o->ch[0]:o->ch[1];
	delete t;
	}else{
	o->s--;
	node *t=o->ch[1];
	while(t->ch[0]->s)t=t->ch[0];
	o->data=t->data;
	erase(o->ch[1],t->data);
	erase_fix(o);
	}return 1;
	}else if(erase(o->ch[o->data<data],data)){
	o->s--,erase_fix(o);
	return 1;
	}else return 0;
	}
	void clear(node *&o){
	if(o->s)clear(o->ch[0]),clear(o->ch[1]),delete(o);
	}
	public:
	aa_tree():nil(new node){
	root=nil->ch[0]=nil->ch[1]=nil;
	}
	~aa_tree(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T&data){
	insert(root,data);
	}
	inline bool erase(const T&data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	node *o=root;
	while(o->s&&o->data!=data)o=o->ch[o->data<data];
	return o->s?1:0;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->s;}
	};
	#endif

view raw aa_tree.cpp hosted with ❤ by GitHub

[ symmetric binary B-tree ( red black tree ) ] 對稱二叉B樹 ( 紅黑樹 )

紅黑樹是一種有最壞情況擔保的平衡樹，其插入最多會用到2個旋轉，刪除最多用到3個旋轉，其犧牲些許的平衡來加快插入刪除的速度，並有高度上h<=2*log(n+1)的保證。
故其插入刪除之常數會比較小，而查詢之常數會較大(相對於一般高度平衡樹)，而一般的函式庫裡的"關聯數組"通常適用紅黑樹實做的
其平衡原因及方法請參考維基百科
因其操作複雜故不常在競賽中被使用，且刪除速度較size balanced tree慢約2倍

以下提供模板:

	/* 0 => black , 1 => red */
	#ifndef SYMMETRIC_BINARY_B_TREE
	#define SYMMETRIC_BINARY_B_TREE
	template<typename T>
	class rb_tree{
	private:
	struct node{
	node *ch[2],*pa;
	int s;
	bool c;
	T data;
	node(const T&d,node*p):pa(p),s(1),c(1),data(d){}
	node():s(0),c(0){ch[0]=ch[1]=pa=this;}
	}*nil,*root;
	inline void rotate(node*o,bool d){
	if(!o->pa->s)root=o->ch[!d];
	else o->pa->ch[o==o->pa->ch[1]]=o->ch[!d];
	o->ch[!d]->pa=o->pa;
	o->pa=o->ch[!d];
	o->ch[!d]=o->pa->ch[d];
	o->pa->ch[d]=o;
	if(o->ch[!d]->s)o->ch[!d]->pa=o;
	o->pa->s=o->s;
	o->s=o->ch[0]->s+o->ch[1]->s+1;
	}
	inline node *find(node*o,const T&data){
	while(o->s&&o->data!=data)o=o->ch[o->data<data];
	return o->s?o:0;
	}
	inline void insert_fix(node*&o){
	while(o->pa->c){
	bool d=o->pa==o->pa->pa->ch[0];
	node *uncle=o->pa->pa->ch[d];
	if(uncle->c){
	uncle->c=o->pa->c=0;
	o->pa->pa->c=1;
	o=o->pa->pa;
	}else{
	if(o==o->pa->ch[d]){
	o=o->pa;
	rotate(o,!d);
	}
	o->pa->c=0;
	o->pa->pa->c=1;
	rotate(o->pa->pa,d);
	}
	}
	root->c=0;
	}
	inline void erase_fix(node*&x){
	while(x!=root&&!x->c){
	bool d=x==x->pa->ch[0];
	node *w=x->pa->ch[d];
	if(w->c){
	w->c=0;
	x->pa->c=1;
	rotate(x->pa,!d);
	w=x->pa->ch[d];
	}
	if(!w->ch[0]->c&&!w->ch[1]->c){
	w->c=1,x=x->pa;
	}else{
	if(!w->ch[d]->c){
	w->ch[!d]->c=0;
	w->c=1;
	rotate(w,d);
	w=x->pa->ch[d];
	}
	w->c=x->pa->c;
	w->ch[d]->c=x->pa->c=0;
	rotate(x->pa,!d);
	break;
	}
	}
	x->c=0;
	}
	void clear(node*&o){
	if(o->s)clear(o->ch[0]),clear(o->ch[1]),delete(o);
	}
	public:
	rb_tree():nil(new node),root(nil){}
	~rb_tree(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T&data){
	node *o=root;
	if(!o->s){
	root=new node(data,nil);
	root->ch[0]=root->ch[1]=nil;
	}else{
	for(;;){
	o->s++;
	bool d=o->data<data;
	if(o->ch[d]->s)o=o->ch[d];
	else{
	o->ch[d]=new node(data,o),o=o->ch[d];
	o->ch[0]=o->ch[1]=nil;
	break;
	}
	}
	}
	insert_fix(o);
	}
	inline bool erase(const T&data){
	node *o=find(root,data);
	if(!o)return 0;
	node *t=o,*p;
	if(o->ch[0]->s&&o->ch[1]->s){
	t=o->ch[1];
	while(t->ch[0]->s)t=t->ch[0];
	}
	p=t->ch[!t->ch[0]->s];
	p->pa=t->pa;
	if(!t->pa->s)root=p;
	else t->pa->ch[t->pa->ch[1]==t]=p;
	if(t!=o)o->data=t->data;
	for(o=t->pa;o->s;o=o->pa)o->s--;
	if(!t->c)erase_fix(p);
	delete t;
	return 1;
	}
	inline bool find(const T&data){
	return find(root,data);
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->s;}
	};
	#endif

view raw red black tree.cpp hosted with ❤ by GitHub

[ scapegoat tree ] 替罪羊樹

替罪羊樹論文其中一個作者，Ronald Linn Rivest，同時也是RSA加密算法的R
替罪羊樹是一種不需要旋轉的平衡樹，它會要求使用者給定一個$\alpha$值，其值介於0.5到1之間，$\alpha$越大，插入就越快，查詢就越慢，實驗發現$\alpha$介於0.55~0.75是最佳的，但是具體的值要看使用者的要求
其平衡滿足:
1.
    $\alpha*size(o) \le size(o \to left)$
    $\alpha*size(o) \le size(o \to right)$
2.
    $deep(o) \le log_{1/\alpha} (size(tree))$
當插入節點導致條件2.不滿足時，會往上回溯到第一個不滿足的節點，然後重建這顆節點的子樹成為完美平衡的二元樹
刪除可以用一般的BST刪除法，當不滿足條件時往上回溯到最先不滿足的節點，重建子樹，也可以懶惰刪除，先設立一個標記， 要刪除這個節點就把它做記號，當$\alpha$*全部的節點數量>實際上存在的節點數量，就重建這顆子樹
本文提供一般刪除的方法

其實替罪羊樹的節點是不需要記錄其size的，因為若需要重建，求得size的時間複雜度與重建的時間複雜度相等，這裡為了支援其它有如rank等操作才將size域附加在節點中

插入刪除的均攤時間複雜度為$\ord{logN}$，相較之下朝鮮樹只是他的劣質仿冒品，其複雜度證明請參考此論文

以下為模板:

	#ifndef SCAPEGOAT_TREE
	#define SCAPEGOAT_TREE
	#include<cmath>
	#include<algorithm>
	template<typename T>
	class scapegoat_tree{
	private:
	struct node{
	node *ch[2];
	int s;
	T data;
	node(const T&d):s(1),data(d){}
	node():s(0){ch[0]=ch[1]=this;}
	}*nil,*root,t;
	const double alpha,loga;
	int maxn;
	inline bool isbad(node*o){
	return o->ch[0]->s>alpha*o->s||o->ch[1]->s>alpha*o->s;
	}
	node *flatten(node *x,node *y){
	if(!x->s)return y;
	x->ch[1]=flatten(x->ch[1],y);
	return flatten(x->ch[0],x);
	}
	node *build(node *x,int n){
	if(!n){
	x->ch[0]=nil;
	return x;
	}
	node *y=build(x,n/2);
	node *z=build(y->ch[1],n-1-n/2);
	y->ch[1]=z->ch[0];
	z->ch[0]=y;
	y->s=y->ch[0]->s+y->ch[1]->s+1;
	return z;
	}
	bool insert(node *&o,const T&data,int dp){
	if(!o->s){
	o=new node(data);
	o->ch[0]=o->ch[1]=nil;
	return dp<=0;
	}
	++o->s;
	if(insert(o->ch[o->data<data],data,--dp)){
	if(!isbad(o))return 1;
	o=build(flatten(o,&t),o->s)->ch[0];
	}return 0;
	}
	inline void success_erase(node *&o){
	node *t=o;
	o=o->ch[0]->s?o->ch[0]:o->ch[1];
	delete t;
	}
	bool erase(node *&o,const T&data){
	if(!o->s)return 0;
	if(o->data==data){
	if(!o->ch[0]->s||!o->ch[1]->s)success_erase(o);
	else{
	--o->s;
	node **t=&o->ch[1];
	for(;(*t)->ch[0]->s;t=&(*t)->ch[0])(*t)->s--;
	o->data=(*t)->data;
	success_erase(*t);
	}return 1;
	}else if(erase(o->ch[o->data<data],data)){
	--o->s;return 1;
	}else return 0;
	}
	void clear(node *&p){
	if(p->s)clear(p->ch[0]),clear(p->ch[1]),delete p;
	}
	public:
	scapegoat_tree(const double&d=0.75):nil(new node),root(nil),alpha(d),loga(log2(1.0/d)),maxn(1){}
	~scapegoat_tree(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;maxn=1;}
	inline void insert(const T&data){
	insert(root,data,std::__lg(maxn)/loga);
	if(root->s>maxn)maxn=root->s;
	}
	inline bool erase(const T&data){
	bool d=erase(root,data);
	if(root->s<alpha*maxn)rebuild();
	return d;
	}
	inline bool find(const T&data){
	node *o=root;
	while(o->s&&o->data!=data)o=o->ch[o->data<data];
	return o->s;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline void rebuild(){root=build(flatten(root,&t),maxn=root->s)->ch[0];}
	inline int size(){return root->s;}
	};
	#endif

view raw scapegoat tree2.cpp hosted with ❤ by GitHub

這裡也提供不需要記錄節點大小(size)的平衡方式，相對的一些需要用到附加size域的函數將不能被使用，但插入刪除的效率不變

	#ifndef NOSIZE_SCAPEGOAT_TREE
	#define NOSIZE_SCAPEGOAT_TREE
	#include<cmath>
	#include<algorithm>
	template<typename T>
	class scapegoat_tree{
	private:
	struct node{
	node *ch[2];
	T data;
	node(const T&d):data(d){ch[0]=ch[1]=0;}
	node(){}
	}*root,t;
	const double alpha,loga;
	int maxn,n;
	int size(node *o){
	return o?size(o->ch[0])+size(o->ch[1])+1:0;
	}
	node *flatten(node *x,node *y){
	if(!x)return y;
	x->ch[1]=flatten(x->ch[1],y);
	return flatten(x->ch[0],x);
	}
	node *build(node *x,int n){
	if(!n){
	x->ch[0]=0;
	return x;
	}
	node *y=build(x,n/2);
	node *z=build(y->ch[1],n-1-n/2);
	y->ch[1]=z->ch[0];
	z->ch[0]=y;
	return z;
	}
	int insert(node *&o,const T&data,int dp){
	if(!o){
	o=new node(data);
	return dp<=0;
	}
	bool d=o->data<data;
	int ds1=insert(o->ch[d],data,--dp);
	if(ds1){
	int ds2=size(o->ch[!d]),n=ds1+ds2+1;
	if(ds1<=alpha*n&&ds2<=alpha*n)return n;
	o=build(flatten(o,&t),n)->ch[0];
	}return 0;
	}
	inline void success_erase(node *&o){
	node *t=o;
	o=o->ch[0]?o->ch[0]:o->ch[1];
	delete t;
	}
	bool erase(node *&o,const T&data){
	if(!o)return 0;
	if(o->data==data){
	if(!o->ch[0]||!o->ch[1])success_erase(o);
	else{
	node **t=&o->ch[1];
	for(;(*t)->ch[0];t=&(*t)->ch[0]);
	o->data=(*t)->data;
	success_erase(*t);
	}return 1;
	}else return erase(o->ch[o->data<data],data);
	}
	void clear(node *&p){
	if(p)clear(p->ch[0]),clear(p->ch[1]),delete p;
	}
	public:
	scapegoat_tree(const double&d=0.75):root(0),alpha(d),loga(log2(1.0/d)),maxn(1),n(0){}
	~scapegoat_tree(){clear(root);}
	inline void clear(){clear(root),root=0;maxn=1,n=0;}
	inline void insert(const T&data){
	insert(root,data,std::__lg(maxn)/loga);
	if(++n>maxn)maxn=n;
	}
	inline bool erase(const T&data){
	bool d=erase(root,data);
	if(d&&--n<alpha*maxn)rebuild();
	return d;
	}
	inline bool find(const T&data){
	node *o=root;
	while(o&&o->data!=data)o=o->ch[o->data<data];
	return o;
	}
	inline const T& min(){
	node *o=root;
	while(o->ch[0])o=o->ch[0];
	return o->data;
	}
	inline const T& max(){
	node *o=root;
	while(o->ch[1])o=o->ch[1];
	return o->data;
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline void rebuild(){root=build(flatten(root,&t),maxn=n)->ch[0];}
	inline int size(){return n;}
	};
	#endif

view raw scapegoat tree1.cpp hosted with ❤ by GitHub

如果編譯發現沒有__lg()這個函數的話，請在程式碼開頭附上這份code:

	namespace std{
	inline int __lg(int n){//Precondition: n >= 0
	int k=0;
	for(;n!=0;n>>=1)++k;
	return k?--k:1;
	}
	}

view raw __lg.cpp hosted with ❤ by GitHub

[ korea tree ] 朝鮮樹

某天albus YY出來了一個沒有旋轉操作的平衡樹，由於albus的真名(华中科技大学附属中学 金正中)聽起來很像朝鮮主席，所以就叫朝鮮樹了！

朝鮮樹是一種不需旋轉的 自平衡二元搜尋樹，利用當深度大於max_deep(通常取sqrt(N))時暴力重建的方式讓平均時間為$\ord{N*sqrt(N)}$，是一種很糟糕的平衡樹，但是我們還是要理解他，因為這樣暴力重建的概念會在其他樹上用到。
若題目的要求數據量<=500000時可以使用

所有操作都與本篇其他樹差不多只多了個rebuild函數，它會暴力將整顆樹重建成深度logN的平衡樹

對於重建的部分則使用了拍扁重建法，將樹先拍扁成鍊表後重建，只需要花費$\ord 1$的空間複雜度

以下提供模板:

	#ifndef KOREA_TREE
	#define KOREA_TREE
	template<typename T>
	class korea_tree{
	private:
	struct node{
	node *ch[2];
	T data;
	int s;
	node(const T&d):data(d),s(1){}
	node():s(0){ch[0]=ch[1]=this;}
	}*nil,*root,t;
	int max_deep,d;
	void insert(node *&o,const T&data){
	if(!o->s){
	o=new node(data);
	o->ch[0]=o->ch[1]=nil;
	}else o->s++,d++,insert(o->ch[o->data<data],data);
	}
	inline void success_erase(node *&o){
	node *t=o;
	o=o->ch[0]->s?o->ch[0]:o->ch[1];
	delete t;
	}
	bool erase(node *&o,const T&data){
	if(!o->s)return 0;
	if(o->data==data){
	if(!o->ch[0]->s||!o->ch[1]->s)success_erase(o);
	else{
	o->s--;
	node **t=&o->ch[1];
	for(;(*t)->ch[0]->s;t=&(*t)->ch[0])(*t)->s--;
	o->data=(*t)->data;
	success_erase(*t);
	}return 1;
	}else if(erase(o->ch[o->data<data],data)){
	o->s--;return 1;
	}else return 0;
	}
	node *flatten(node *x,node *y){
	if(!x->s)return y;
	x->ch[1]=flatten(x->ch[1],y);
	return flatten(x->ch[0],x);
	}
	node *build(node *x,int n){
	if(!n){
	x->ch[0]=nil;
	return x;
	}
	node *y=build(x,n/2);
	node *z=build(y->ch[1],n-1-n/2);
	y->ch[1]=z->ch[0];
	z->ch[0]=y;
	y->s=y->ch[0]->s+y->ch[1]->s+1;
	return z;
	}
	void clear(node *&p){
	if(p->s)clear(p->ch[0]),clear(p->ch[1]),delete p;
	}
	public:
	korea_tree(int d=1000):nil(new node),root(nil),max_deep(d){}
	~korea_tree(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T&data){
	d=0;
	insert(root,data);
	if(d>=max_deep)rebuild();
	}
	inline bool erase(const T&data){
	return erase(root,data);
	}
	inline void rebuild(){
	root=build(flatten(root,&t),root->s)->ch[0];
	}
	inline bool find(const T&data){
	node *o=root;
	while(o->s&&o->data!=data)o=o->ch[o->data<data];
	return o->s;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->s;}
	};
	#endif

view raw korea_tree.cpp hosted with ❤ by GitHub

[ Splay Tree ] 伸展樹模板

最近研究了伸展樹，它可以在均攤$\ord{logN}$的時間完成伸展操作，所謂的伸展操作就是將某個節點旋轉(rotate)到根(root)的過程，因為我們大部分在搜尋資料時，常常都只會用到特定的資料，其他大多數的資料是較少被搜尋的，因此將常被搜尋的資料靠近根是很好的策略，伸展樹由此被發明。
因為伸展樹可以分裂合併，所以可以取代分裂合併式treap或是分裂合併式randomize binary tree的功能，可以在競賽中被使用，而普通的伸展樹效率不高，故本文只介紹作區間處理的伸展樹
之後的link-cut tree會用到他，所以必須要能理解
例題: UVA 11922

定義伸展樹節點:(T為資料型態)

	struct node{
	T data;
	node *ch[2];
	int s;
	node(const T&d):data(d), s(1){}
	node():s(0){}
	/*cmp是很重要的一個函式一定要理解*/
	char cmp(int v){
	int d=v-ch[0]->s;
	if(d==1)return -1;
	return d>0;
	}
	void down(){
	/*可以加一些下推懶惰標記*/
	}
	void up(){
	s=ch[0]->s+ch[1]->s+1;
	/*可以加一些上推懶惰標記*/
	}
	}*root,*nil;

view raw splay_tree_node.cpp hosted with ❤ by GitHub

初始化時nil的l,r要指向自己

	nil=new node();
	root=nil->ch[0]=nil->ch[1]=nil;

view raw splay_tree_nil.cpp hosted with ❤ by GitHub

旋轉(跟其他樹差不多)
d=0左旋,1右旋

	void rotate(node *&a,bool d){
	node *b=a;
	a=a->ch[!d];
	b->ch[!d]=a->ch[d];
	a->ch[d]=b;
	b->up();
	a->up();
	}

view raw splay_tree_rotate.cpp hosted with ❤ by GitHub

接下來是最重要的伸展操作
這是按劉汝佳在橘色那本書的寫法實現的
主要是將第k個節點伸展到root的過程
注意!! k必須大於0

	void splay(node *&o,int k){
	o->down();
	char d1=o->cmp(k);
	if(~d1){
	if(d1)k-=o->ch[0]->s+1;
	node *p=o->ch[d1];
	p->down();
	char d2=p->cmp(k);
	if(~d2){
	int k2=d2?k-p->ch[0]->s-1:k;
	splay(p->ch[d2],k2);
	if(d1==d2)rotate(o,d1^1);
	else rotate(o->ch[d1],d1);
	}
	rotate(o,d1^1);
	}
	}

view raw splay_tree_splay.cpp hosted with ❤ by GitHub

有了splay操作，接下來是分裂&合併

	/* l不能為nil */
	node *merge(node *l,node *r){
	splay(l,l->s);
	l->ch[1]=r;
	l->up();
	return l;
	}
	/* k必須 > 0 */
	void split(node *o,node *&l,node *&r,int k){
	splay(o,k);
	l=o;
	r=o->ch[1];
	o->ch[1]=nil;
	l->up();
	}

view raw splay_tree_merge_split.cpp hosted with ❤ by GitHub

當我們要維護一個區間1~n時，可以建立0~n的splay tree，0是廢棄節點，較為方便使用

最後提供新增節點的方法:

	inline node *new_node(const T&d){
	node *o=new node(d);
	o->l=o->r=nil;
	}

view raw splay_tree_new_node.cpp hosted with ❤ by GitHub

這樣所有的操作就齊全了
剩下的一些懶惰標記跟treap一樣，可以做區間翻轉和線段樹的功能
但是請小心指標的使用
(nil為衛兵指標，方便用來編輯，若不想用可以將程式碼略為修改)

[ Big Interger ] 大數模板

今天在家裡寫了一整天的大數，好不容易加減乘除都有了，但是乘法的部分FFT還不會寫所以先用做基本的n^2乘法(聽說有一個奇怪的方法叫做Karatsuba演算法也能做大數乘法，還蠻快的)

定義一個大數:
bigN a;//定義大數

以下是大數模板:

	#include<algorithm>
	#include<iostream>
	#include<sstream>
	#include<iomanip>
	#include<vector>
	#include<string>
	#include<cmath>
	using namespace std;
	template<typename T>
	inline string to_string(const T& x){
	stringstream ss;
	return ss<<x,ss.str();
	}
	typedef long long LL;
	struct bigN:vector<LL>{
	const static int base=1000000000,width=log10(base);
	bool negative;
	bigN(const_iterator a,const_iterator b):vector<LL>(a,b){}
	bigN(string s){
	if(s.empty())return;
	if(s[0]=='-')negative=1,s=s.substr(1);
	else negative=0;
	for(int i=int(s.size())-1;i>=0;i-=width){
	LL t=0;
	for(int j=max(0,i-width+1);j<=i;++j)
	t=t*10+s[j]-'0';
	push_back(t);
	}
	trim();
	}
	template<typename T>
	bigN(const T &x):bigN(to_string(x)){}
	bigN():negative(0){}
	void trim(){
	while(size()&&!back())pop_back();
	if(empty())negative=0;
	}
	void carry(int _base=base){
	for(size_t i=0;i<size();++i){
	if(at(i)>=0&&at(i)<_base)continue;
	if(i+1u==size())push_back(0);
	int r=at(i)%_base;
	if(r<0)r+=_base;
	at(i+1)+=(at(i)-r)/_base;
	at(i)=r;
	}
	}
	int abscmp(const bigN &b)const{
	if(size()>b.size())return 1;
	if(size()<b.size())return -1;
	for(int i=int(size())-1;i>=0;--i){
	if(at(i)>b[i])return 1;
	if(at(i)<b[i])return -1;
	}
	return 0;
	}
	int cmp(const bigN &b)const{
	if(negative!=b.negative)return negative?-1:1;
	return negative?-abscmp(b):abscmp(b);
	}
	bool operator<(const bigN&b)const{return cmp(b)<0;}
	bool operator>(const bigN&b)const{return cmp(b)>0;}
	bool operator<=(const bigN&b)const{return cmp(b)<=0;}
	bool operator>=(const bigN&b)const{return cmp(b)>=0;}
	bool operator==(const bigN&b)const{return !cmp(b);}
	bool operator!=(const bigN&b)const{return cmp(b)!=0;}
	bigN abs()const{
	bigN res=*this;
	return res.negative=0, res;
	}
	bigN operator-()const{
	bigN res=*this;
	return res.negative=!negative,res.trim(),res;
	}
	bigN operator+(const bigN &b)const{
	if(negative)return -(-(*this)+(-b));
	if(b.negative)return *this-(-b);
	bigN res=*this;
	if(b.size()>size())res.resize(b.size());
	for(size_t i=0;i<b.size();++i)res[i]+=b[i];
	return res.carry(),res.trim(),res;
	}
	bigN operator-(const bigN &b)const{
	if(negative)return -(-(*this)-(-b));
	if(b.negative)return *this+(-b);
	if(abscmp(b)<0)return -(b-(*this));
	bigN res=*this;
	if(b.size()>size())res.resize(b.size());
	for(size_t i=0;i<b.size();++i)res[i]-=b[i];
	return res.carry(),res.trim(),res;
	}
	bigN operator*(const bigN &b)const{
	bigN res;
	res.negative=negative!=b.negative;
	res.resize(size()+b.size());
	for(size_t i=0;i<size();++i)
	for(size_t j=0;j<b.size();++j)
	if((res[i+j]+=at(i)*b[j])>=base){
	res[i+j+1]+=res[i+j]/base;
	res[i+j]%=base;
	}//乘法用carry會溢位
	return res.trim(),res;
	}
	/* 用二分搜做除法很慢
	bigN operator/(const bigN &b)const{
	bigN res;
	if(size()<b.size())return res;
	if(negative||b.negative){
	res=abs()/b.abs();
	res.negative=negative!=b.negative;
	return res.trim(),res;
	}
	res.resize(size()-b.size()+1);
	for(int i=int(res.size())-1;i>=0;--i){
	int l=0,r=base-1,ans=l,mid;
	while(l<=r){
	res[i]=mid=(l+r)/2;
	if((res*b).abscmp(*this)>0)r=mid-1;
	else l=mid+1,ans=mid;
	}
	res[i]=ans;
	}
	return res.carry(),res.trim(),res;
	}
	/*/
	bigN operator/(const bigN &b)const{
	int norm=base/(b.back()+1);
	bigN x=abs()*norm;
	bigN y=b.abs()*norm;
	bigN q,r;
	q.resize(x.size());
	for(int i=int(x.size())-1;i>=0;--i){
	r=r*base+x[i];
	int s1=r.size()<=y.size()?0:r[y.size()];
	int s2=r.size()<y.size()?0:r[y.size()-1];
	int d=(LL(base)*s1+s2)/y.back();
	r=r-y*d;
	while(r.negative)r=r+y,--d;
	q[i]=d;
	}
	q.negative=negative!=b.negative;
	return q.trim(),q;
	}
	//*/
	bigN operator%(const bigN &b)const{
	return *this-(*this/b)*b;
	}
	friend istream& operator>>(istream &ss,bigN &b){
	string s;
	return ss>>s, b=s, ss;
	}
	friend ostream& operator<<(ostream &ss,const bigN &b){
	if(b.negative)ss<<'-';
	ss<<(b.empty()?0:b.back());
	for(int i=int(b.size())-2;i>=0;--i)
	ss<<setw(width)<<setfill('0')<<b[i];
	return ss;
	}
	template<typename T>
	operator T(){
	stringstream ss;
	ss<<*this;
	T res;
	return ss>>res,res;
	}
	};

view raw bigN.cpp hosted with ❤ by GitHub

[ AVL Tree ] 優化AVL樹

AVL樹是以子樹的高度做為平衡的條件
其平衡條件為:
     左子樹的高度與右子樹的高度差不能超過1

因為每次的插入及刪除都有可能壞平衡，所以必須要進行旋轉來修正其平衡性。
共有4種旋轉方法，一序為 : 左左、左右、右右、右左，若想知道關於如何進行平衡的介紹請點擊這裡，其插入刪除常數較大，故其平均效率比Size Balanced Tree要差。

因為在網路上看到很多的AVL樹，其寫法都過於冗長，所以我這邊發一個AVL樹的模板，其code複雜度跟Size Balanced Tree差不多。
主要是將相同的平衡操作合併成一個balanced()函數，插入和刪除完後呼叫即可。

以下為模板:

	#ifndef AVL_TREE
	#define AVL_TREE
	template<typename T>
	class avl_tree{
	private:
	struct node{
	node *ch[2];
	int size;
	char h;
	T data;
	inline void up(){h=ch[ch[0]->h<ch[1]->h]->h+1;}
	node(const T&d):size(1),h(1),data(d){}
	node():size(0),h(0){ch[0]=ch[1]=this;}
	}*nil,*root;
	inline void rotate(node *&a,bool d){
	node *b=a;
	a=a->ch[!d];
	b->ch[!d]=a->ch[d];
	a->ch[d]=b;
	a->size=b->size;
	b->size=b->ch[0]->size+b->ch[1]->size+1;
	b->up(),a->up();
	}
	inline void balanced(node *&o,bool d){
	if(o->ch[d]->h-o->ch[!d]->h>1){
	node *&t=o->ch[d];
	if(t->ch[d]->h>t->ch[!d]->h)rotate(o,!d);
	else rotate(o->ch[d],d),rotate(o,!d);
	}else o->up();
	}
	void insert(node *&o,const T &data){
	if(!o->size){
	o=new node(data);
	o->ch[0]=o->ch[1]=nil;
	}else{
	o->size++;
	bool d=o->data<data;
	insert(o->ch[d],data);
	balanced(o,d);
	}
	}
	bool erase(node *&o,const T &data){
	if(!o->size)return 0;
	if(o->data==data){
	if(!o->ch[0]->size||!o->ch[1]->size){
	node *t=o;
	o=o->ch[0]->size?o->ch[0]:o->ch[1];
	delete t;
	}else{
	o->size--;
	node *tmd=o->ch[1];
	while(tmd->ch[0]->size)tmd=tmd->ch[0];
	o->data=tmd->data;
	erase(o->ch[1],tmd->data);
	balanced(o,0);
	}return 1;
	}
	bool d=o->data<data;
	if(erase(o->ch[d],data)){
	o->size--,balanced(o,!d);
	return 1;
	}else return 0;
	}
	void clear(node *&o){
	if(o->size)clear(o->ch[0]),clear(o->ch[1]),delete(o);
	}
	public:
	avl_tree():nil(new node),root(nil){}
	~avl_tree(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T &data){
	insert(root,data);
	}
	inline bool erase(const T &data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	node *o=root;
	while(o->size&&o->data!=data)o=o->ch[o->data<data];
	return o->size;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->size;)
	if(o->data<data)cnt+=o->ch[0]->size+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->size)o=o->ch[0];
	else if(k==o->ch[0]->size+1)return o->data;
	else k-=o->ch[0]->size+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->size)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->size)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->size;}
	inline int height(){return root->h;}
	};
	#endif

view raw AVL tree.cpp hosted with ❤ by GitHub

Treap 模板

這是 旋轉式Treap 的模板，是目前所有平衡樹中code最為簡便的一種。
插入和刪除的時間複雜度跟Randomized binary search tree一樣，是在演算法競賽中最適合使用的平衡樹。
因為 分裂/合併式 Treap 可以完全被 分裂/合併式 Randomized binary search tree取代，故不提供

做法1，完全依靠旋轉進行平衡(刪除速度較慢):

	#ifndef TREAP
	#define TREAP
	template<typename T>
	class treap{
	private:
	struct node{
	T data;
	unsigned fix;
	int size;
	node *ch[2];
	node(const T&d):data(d),size(1){}
	node():size(0){ch[0]=ch[1]=this;}
	}*nil,*root;
	unsigned x;
	inline unsigned ran(){
	return x=x*0xdefaced+1;
	}
	inline void rotate(node *&a,bool d){
	node *b=a;
	a=a->ch[!d];
	a->size=b->size;
	b->ch[!d]=a->ch[d];
	a->ch[d]=b;
	b->size=b->ch[0]->size+b->ch[1]->size+1;
	}
	void insert(node *&o,const T &data){
	if(!o->size){
	o=new node(data),o->fix=ran();
	o->ch[0]=o->ch[1]=nil;
	}else{
	o->size++;
	bool d=o->data<data;
	insert(o->ch[d],data);
	if(o->ch[d]->fix>o->fix)rotate(o,!d);
	}
	}
	bool success_erase(node *&o){
	if(!o->ch[0]->size||!o->ch[1]->size){
	node *t=o;
	o=o->ch[0]->size?o->ch[0]:o->ch[1];
	delete t;
	}else{
	o->size--;
	bool d=o->ch[0]->fix>o->ch[1]->fix;
	rotate(o,d);
	success_erase(o->ch[d]);
	}
	return 1;
	}
	bool erase(node *&o,const T &data){
	if(!o->size)return 0;
	if(o->data==data)return success_erase(o);
	if(erase(o->ch[o->data<data],data)){
	o->size--;return 1;
	}else return 0;
	}
	void clear(node *&o){
	if(o->size)clear(o->ch[0]),clear(o->ch[1]),delete o;
	}
	public:
	treap(unsigned s=20150119):nil(new node),root(nil),x(s){}
	~treap(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T &data){
	insert(root,data);
	}
	inline bool erase(const T &data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	for(node *o=root;o->size;)
	if(o->data==data)return 1;
	else o=o->ch[o->data<data];
	return 0;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->size;)
	if(o->data<data)cnt+=o->ch[0]->size+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->size)o=o->ch[0];
	else if(k==o->ch[0]->size+1)return o->data;
	else k-=o->ch[0]->size+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->size)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->size)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->size;}
	};
	#endif

view raw treap1.cpp hosted with ❤ by GitHub

作法2，刪除時合併左右子樹(刪除速度較快):

	#ifndef TREAP
	#define TREAP
	template<typename T>
	class treap{
	private:
	struct node{
	T data;
	unsigned fix;
	int s;
	node *ch[2];
	node(const T&d):data(d),s(1){}
	node():s(0){ch[0]=ch[1]=this;}
	}*nil,*root;
	unsigned x;
	inline unsigned ran(){
	return x=x*0xdefaced+1;
	}
	inline void rotate(node *&a,bool d){
	node *b=a;
	a=a->ch[!d];
	a->s=b->s;
	b->ch[!d]=a->ch[d];
	a->ch[d]=b;
	b->s=b->ch[0]->s+b->ch[1]->s+1;
	}
	void insert(node *&o,const T &data){
	if(!o->s){
	o=new node(data),o->fix=ran();
	o->ch[0]=o->ch[1]=nil;
	}else{
	o->s++;
	bool d=o->data<data;
	insert(o->ch[d],data);
	if(o->ch[d]->fix>o->fix)rotate(o,!d);
	}
	}
	node *merge(node *a,node *b){
	if(!a->s||!b->s)return a->s?a:b;
	if(a->fix>b->fix){
	a->ch[1]=merge(a->ch[1],b);
	a->s=a->ch[0]->s+a->ch[1]->s+1;
	return a;
	}else{
	b->ch[0]=merge(a,b->ch[0]);
	b->s=b->ch[0]->s+b->ch[1]->s+1;
	return b;
	}
	}
	bool erase(node *&o,const T &data){
	if(!o->s)return 0;
	if(o->data==data){
	node *t=o;
	o=merge(o->ch[0],o->ch[1]);
	delete t;
	return 1;
	}
	if(erase(o->ch[o->data<data],data)){
	o->s--;return 1;
	}else return 0;
	}
	void clear(node *&o){
	if(o->s)clear(o->ch[0]),clear(o->ch[1]),delete o;
	}
	public:
	treap(unsigned s=20150119):nil(new node),root(nil),x(s){}
	~treap(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T &data){
	insert(root,data);
	}
	inline bool erase(const T &data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	for(node *o=root;o->s;)
	if(o->data==data)return 1;
	else o=o->ch[o->data<data];
	return 0;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->s;}
	};
	#endif

view raw treap2.cpp hosted with ❤ by GitHub

[ 線性同餘方法 ] 亂數產生器 實做

今天來講幾個簡單的亂數產生方法
主要是利用線性同餘方法來處理的

它的首要條件便是必須符合平均分配率Uniform distribution，也就是在$0 \sim m-1$之間的每個數字, 出現的機率必須相等。
Linear Congruential Method雖是常用法，但也有它的條件:
$X(n+1) = (a \times X(n) + c) \%m$，其中$\%$是取餘數的意思
$X(n+1)$為新的亂數，$X(n)$為前一次產生的亂數。則各參數必須符合下列條件:

1. $c$與$m$必須互質。
2. 對於任何$m$的質因數$p$，也必須為$a-1$的質因數。
3. 若$m$為$4$的倍數，則$a-1$也必須為4的倍數。

如此才能確保產生出來的亂數符合平均分配率，同時也具有最長的重複週期。

以下是在網路上找到的一些不錯的亂數產生方法:

	//線性同餘方法
	//這個方法數字很好記
	unsigned random1(){
	static unsigned x=20150118;
	return x=x*0xdefaced+1;
	}
	int random2(){
	static unsigned x=time(0);
	return x=x*16807;//7^5
	}
	// 網路上看到的，用來處理treap的亂數
	int random3(){
	static int x=123456789;
	return x+=(x<<2)+1;
	}
	long long random4(){
	static long long x=323232;
	return x+=x<<2|1;
	}
	//5、6是較快的算法，利用簡單的移位處理
	inline int random5(){
	static int seed=20160424;
	return seed+=(seed<<16)+0x1db3d743;
	}
	inline long long random6(){
	static long long seed=20160424;
	return seed+=(seed<<32)+0xdb3d742c265539d;
	}
	//要記數字太麻煩了而且比較慢，比賽時不要用這種
	int random7(){
	static int x=1;
	return x=x*1103515245+12345;
	}

view raw lehmer random.cpp hosted with ❤ by GitHub

其中random1是Brain大神在處理randomize binary tree的亂數(0xdefaced是質數，defaced是英文單字比較好記)
random2聽說是stdlib裡面的rand()原型，而16807是7的5次方，random2可以產生較rand()大的亂數
random4不是線性同於方法，個人覺得他比較像Subtract with carry的方法
在做[分裂/合併式]randomize binary tree時，其實可以利用seed++的方式當作亂數用(連結中有教學)

c++ rope 基本應用

rope(粗繩)是一個很強大的字串處理物件，聽說是為了與string(細繩)做對比才叫做rope的，其底層是一棵持久化平衡樹，因此在做合併或在中間插入的期望時間是logn

再使用rope前必須在前置作以下的處理

#include<ext/rope>
using namespace __gnu_cxx;

定義一個元素為char的rope:

rope<char> r;

其實rope本身為了方便已經將rope<char>重新定義為crope
因此也可以這樣寫:

crope r;

當然rope也可以存取非char的元素
只是有些操作會不支援(例如輸出操作)

對於rope的一些常用操作請參考由SGI網站翻譯的rope成員詳細說明
裡面整理了rope的各種函式及元素
且已經將其翻譯成中文方便閱讀

Randomized Binary Search Tree 平均深度證明

一棵Randomized Binary Search Tree ，平均深度為2logn
其每個節點成為該顆子樹的根的機率有1/(該顆子樹的大小+1)
對於[分裂/合併式][Split / Merge] Randomized Binary Search Tree 亦同
而其等價於在BST中進行隨機插入
這裡提供嚴謹的證明 連結
以下簡略的證明原自於陳柏叡大大(數奧國手)親筆所寫:

現在總共有n個數，要證說對他隨機排序之後一個一個插入BST得到的平均最大深度=2logn。
那看任意一個數x，他的祖先個數就是它的深度，所以只要算x的祖先個數的期望直。
這只要算，對x以外的每個數y，y當x的祖先的機率是多少，再把他們加起來即可。

但注意到，假設在這n個數(由小到大)排序好的序列裡面，x 排在第 i 個， y 排在第 j 個(  y 可以 <x，所以 j 可以 < i  )，那麼落在 [ i , j ] ( 或 [ j , i ] )這個區間裡的所有數， y 是第一個被插入BST的(否則會和 y 是 x 的祖先矛盾)，而這件事情發生的機率為 1/(  | i - j | +1 ) 。
也就是，如果一個數 y 和 x 在排序好的序列裡面的坐標差 = d 的話，那麼 y 當 x 的祖先的機率就是 1/(d+1)。

所以如果 x 是排序好的序列裡的第 i 個，那所求是
1 / i + 1 / ( i - 1 ) + ... + 1/2 + 1/2 + 1/3 + ... + 1/( n - i + 1 )
<= 2 ( 1/2 + 1/3 + ... + 1/n )
約 = 2logn

[smart pointer] [reference counter] 智能指標 引用計數

鑒於有些題目需要記憶體的限制，因此必須要考慮到記憶體的管理。
而智能指標可以幫助我們解決這個問題

在C++11，STL已經有內件shared_ptr，但是速度很慢，而大部分在比賽時會用到的只有引用計數的概念(可以參考Brain大神的blog)，因此我將引用計數(參考自 C++ Template 侯捷譯)部分特別獨立出來，做了一份模板

以下為模板:

	#ifndef REFERENCE_POINTER
	#define REFERENCE_POINTER
	template<typename T>
	struct _RefCounter{
	T data;
	int ref;
	_RefCounter(const T&d=0):data(d),ref(0){}
	};
	template<typename T>
	struct reference_pointer{
	_RefCounter<T> *p;
	T *operator->(){return &p->data;}
	T &operator*(){return p->data;}
	operator _RefCounter<T>*(){return p;}
	reference_pointer &operator=(const reference_pointer &t){
	if(p&&!--p->ref)delete p;
	p=t.p;
	p&&++p->ref;
	return *this;
	}
	reference_pointer(_RefCounter<T> *t=0):p(t){
	p&&++p->ref;
	}
	reference_pointer(const reference_pointer &t):p(t.p){
	p&&++p->ref;
	}
	~reference_pointer(){
	if(p&&!--p->ref)delete p;
	}
	};
	template<typename T>
	inline reference_pointer<T> new_reference(const T&nd){
	return reference_pointer<T>(new _RefCounter<T>(nd));
	}
	#endif

view raw reference_counting.cpp hosted with ❤ by GitHub

用法:

reference_pointer<int> a;//建立一個int的引用指標a
a = new_reference(5);//a指向新增int動態變數，其值為5
a = new_reference<int>(5);//同上，只是定義較為嚴謹
a = new_reference((int)5);//同上，只是定義較為嚴謹
reference_pointer<int> b = a;//將b指向a所指向之物

struct P{
     int a,b;
     P(int _a,int _b):a(_a),b(_b){}
}p(2,3);
reference_pointer<P> a;//建立一個P的引用指標a
c = new_reference(P(1,2));//指向新增P動態變數，其值為1,2
c = new_reference<P>(P(1,2));//同上，只是定義較為嚴謹
c = new_reference(p);//指向新增P動態變數，其值為p

其他的用法可以由下面的代碼看出來(HOJ Problem : 226 - K. CP AC代碼):

	#include<bits/stdc++.h>
	using namespace std;
	#ifndef REFERENCE_POINTER
	#define REFERENCE_POINTER
	template<typename T>
	struct _RefCounter{
	T data;
	int ref;
	_RefCounter(const T&d=0):data(d),ref(0){}
	};
	template<typename T>
	struct reference_pointer{
	_RefCounter<T> *p;
	T *operator->(){return &p->data;}
	T &operator*(){return p->data;}
	operator _RefCounter<T>*(){return p;}
	reference_pointer &operator=(const reference_pointer &t){
	if(p&&!--p->ref)delete p;
	p=t.p;
	p&&++p->ref;
	return *this;
	}
	reference_pointer(_RefCounter<T> *t=0):p(t){
	p&&++p->ref;
	}
	reference_pointer(const reference_pointer &t):p(t.p){
	p&&++p->ref;
	}
	~reference_pointer(){
	if(p&&!--p->ref)delete p;
	}
	};
	template<typename T>
	inline reference_pointer<T> new_reference(const T&nd){
	return reference_pointer<T>(new _RefCounter<T>(nd));
	}
	#endif
	struct node;
	typedef reference_pointer<node> pt;
	struct node{
	char data;
	int size;
	pt l,r;
	node(const node&t):data(t.data),size(t.size),l(t.l),r(t.r){}
	node(const char &d):data(d),size(1){}
	inline void up(){
	size=1;
	if(l)size+=l->size;
	if(r)size+=r->size;
	}
	};
	inline int size(pt &o){return o?o->size:0;}
	void split(pt o,pt &a,pt &b,int k){
	if(!o)a=b=0;
	else{
	o=new_reference(*o);
	if(k<=size(o->l)){
	b=o;
	split(o->l,a,b->l,k);
	}else{
	a=o;
	split(o->r,a->r,b,k-size(o->l)-1);
	}
	o->up();
	}
	}
	pt merge(pt a,pt b){
	if(!a||!b)return a?a:b;
	static int x;
	if(x++%(a->size+b->size)<a->size){
	a->r=merge(a->r,b);
	a->up();
	return a;
	}else{
	b->l=merge(a,b->l);
	b->up();
	return b;
	}
	}
	pt build(char *s,int l,int r){
	int mid=(l+r)>>1;
	pt k=new_reference(node(s[mid]));
	if(l<=mid-1)k->l=build(s,l,mid-1);
	if(mid+1<=r)k->r=build(s,mid+1,r);
	k->up();
	return k;
	}
	void dfs(pt&t){
	if(!t)return;
	dfs(t->l);
	putchar(t->data);
	dfs(t->r);
	}
	char s[1000005];
	pt root;
	pt a,b,c,d,e,f;
	int n,m;
	int main(){
	scanf("%d%s%d",&m,s,&n);
	root=build(s,0,strlen(s)-1);
	while(n--){
	int l,r,x;
	scanf("%d%d%d",&l,&r,&x);
	split(root,b,c,++r);
	split(b,a,b,l);
	split(root,d,e,x);
	root=merge(merge(d,b),e);
	split(root,root,f,m);
	}
	dfs(root);puts("");
	return 0;
	}

view raw reference_counting_example.cpp hosted with ❤ by GitHub

其使用了[持久化][分裂/合併式]隨機二分搜尋樹
持久化的部分會在之後介紹

[分裂/合併式] 隨機二分查找樹 [Split / Merge] Randomized Binary Search Tree

今天來講解可以做區間處理的Randomized Binary Search Tree
好比treap，可以利用拆分跟合併的方式來做區間處理
感謝伯恩大神的指導，請參考其網站

寫起來並不會比需要旋轉的難
先定義他的節點:

	struct node{
	T data;
	int s;
	//可以再加入其他維護值(max,min等等)做區間處理用
	node *l,*r;
	node(const T &d):data(d),s(1),l(0),r(0){}
	inline void up(){
	s=1;
	if(l)s+=l->s;
	if(r)s+=r->s;
	}
	inline void down(){
	//這邊是放作區間維護需要的更新(懶惰標記下推)
	//(min,max,是否反轉區間....等等)
	//本範例只提供最基本的處理
	}
	}*root;

view raw rbst_node.cpp hosted with ❤ by GitHub

可以在裡面增加一些懶惰標記的更新(min,max,rev...等等)

計算節點大小

	inline int size(node *o){
	return o?o->s:0;
	}

view raw rbst_size.cpp hosted with ❤ by GitHub

主要利用[分裂/合併]來進行區間處理
接下來講解分裂(參見 treap分裂):

	void split(node *o,node *&a,node *&b,int k){
	if(!o)a=b=0;
	else{
	o->down();
	if(k<=size(o->l)){
	b=o;
	split(o->l,a,b->l,k);
	}else{
	a=o;
	split(o->r,a->r,b,k-size(o->l)-1);
	}
	o->up();
	}
	}

view raw rbst_split.cpp hosted with ❤ by GitHub

這是將o按中序將前k個節點分到a，剩下的分到b

合併(參見 treap合併):

	node *merge(node *a,node *b){
	if(!a||!b)return a?a:b;
	static int x;
	if(x++%(a->s+b->s)<a->s){
	a->down();
	a->r=merge(a->r,b);
	a->up();
	return a;
	}else{
	b->down();
	b->l=merge(a,b->l);
	b->up();
	return b;
	}
	}

view raw rbst_merge.cpp hosted with ❤ by GitHub

這就是隨機的精華所在
每次merge時，a有size(a)/(size(a)+size(b))的機率將a作為root
b反之亦然
而本來random是這樣寫的:

	inline int ran(){
	static int x=20160425;
	return (x=0xdefaced*x+1)&INT_MAX;
	}

view raw rbst_ran.cpp hosted with ❤ by GitHub

但是看了丁安立魔法師的blog，發現只要x++就好了，也不知道為甚麼

如果只是要做普通的二分搜尋樹，可以利用排名來求出其在整顆樹中的名次，然後再相對應處做插入

	inline int rank(const T &data){
	//使用前必須保證整顆樹是BST
	//否則會run time error
	node *p=root;
	int cnt=0;
	while(p){
	if(data<=p->data)p=p->l;
	else cnt+=size(p->l)+1,p=p->r;
	}
	return cnt;
	}
	inline void insert(const T &data,int k){
	//在第k個位置之前插入data
	node *a,*b,*now;
	split(root,a,b,k);
	now=new node(data);
	a=merge(a,now);
	root=merge(a,b);
	}

view raw rbst_rank_insert.cpp hosted with ❤ by GitHub

插入data可以寫成:
bst.insert(data,bst.rank(data));
刪除也可以用類似的方法處理

旋轉式的randomize bst就可以做這些處理了(其實分裂合併式bst速度有點慢)，那為何還需要 分裂/合併 呢?
我們可以把整顆樹的中序走訪看成是一條陣列
-->因此可以利用它來作區間的維護

像是區間的最大最小值，區間反轉啦，區間移動、刪除等等都可以利用它在logn的時間完成
只須修改其up跟down函數就可以達成目的
凡是線段樹能做到的事情，他都能做到，可以用它來解決很複雜的區間處理問題

注意的是，在有懶惰標記的情況下，如果要詢問一個區間必須在切出那塊區間後先進行down的操作再進行詢問，否則懶惰標記可能沒有被推下去，範例如下:

	/*需要再node裡面增加tag跟ma用來儲存區間增加量跟最大值*/
	int ask(node *&o,int l,int r){
	/*在o子樹詢問區間[l,r]的最大值*/
	node *a,*b,*c;
	split(o,a,b,l-1);
	split(b,b,c,r-l+1);
	b->down();/*記得要先將懶惰標記下推*/
	int ans=b->ma;
	o=merge(a,merge(b,c));
	return ans;
	}
	int add(node *&o,int l,int r,int v){
	/*在o子樹將區間[l,r]的值都加上v*/
	node *a,*b,*c;
	split(o,a,b,l-1);
	split(b,b,c,r-l+1);
	b->tag+=v;
	o=merge(a,merge(b,c));
	}

view raw rbst_ask_add.cpp hosted with ❤ by GitHub

這個範例up()跟down()也要重寫:

	inline void up(){
	size=1;
	ma=data;
	if(l){
	s+=l->s;
	ma=max(ma,l->ma+l->tag);
	}
	if(r){
	s+=r->s;
	ma=max(ma,r->ma+r->tag);
	}
	}
	inline void down(){
	data+=tag;
	if(l)l->tag+=tag;
	if(r)r->tag+=tag;
	tag=0;
	}

view raw rbst_tag_up_down.cpp hosted with ❤ by GitHub

因為應用時分靈活，故不提供樣板，只做介紹

[ Randomized Binary Search Tree ] 隨機二分查找樹

今天來介紹隨機二分查找樹，因為看了很多外國的文章看了半天都看不懂，現在好不容易看懂英文了，就把它用中文的方式來介紹吧!
隨機二分查找樹，是利用隨機數或機率分布來達到平衡的條件，其證明以附在連結中。
利用隨機數的方法包括Treap，但是這太常見了，所以掠過。

這邊要介紹不需額外再節點紀錄隨機值的方式，透過其節點大小來計算其是否為根或是根據節點大小的比例隨機插入。

維持平衡用左旋根右旋(參見:樹旋轉):
接下來來講解一下它的插入
這裡會用到兩個函數:
insert():插入節點
insert_as_root(node *&p,T k):將k插入子樹p並做為p的根節點
為甚麼要這樣做呢?
當在對節點p進行插入時，新插入點根據隨機分布會有1/(size(p)+1)的機率成為這顆子樹的根
因此可以以遞迴處理
若此節點的值>插入值，則插入其左節點，然後右旋，反之亦然
其常數很小，所以常常插入時間比size balanced tree快，但是深度太大，刪除時間較慢
其平均深度等價於在一棵普通的BST進行隨機插入最終的平均深度，因此為2nlogn

而對於刪除節點，和Treap一樣，有兩種不同的做法:

第一種是透過與左偏樹類似的方法刪除，不需要旋轉，通常效率交高
合併左右子樹，然後刪除其根節點
合併的方式利用了隨機分布的概念
假設需將a,b合併
則有size(a)/(size(a)+size(b))的機率將b合併到a子樹
反之有size(b)/(size(a)+size(b))的機率將a合併到b子樹

第二種作法完全依賴於旋轉
在刪除時如果要被刪除的節點不是一條鏈的話，就按造隨機的概念進行旋轉
有size(左子樹)/(size(左子樹)+size(右子樹))的機率進行右旋
反之有size(右子樹)/(size(左子樹)+size(右子樹))的機率進行左旋

以下為第一種作法(效率較高):

	#ifndef RANDOMIZED_BINARY_SEARCH_TREE
	#define RANDOMIZED_BINARY_SEARCH_TREE
	template<typename T>
	class random_bst{
	private:
	struct node{
	T data;
	unsigned size;
	node *ch[2];
	node(const T &d):data(d),size(1){ch[0]=ch[1]=0;}
	inline void up(){
	size=(ch[0]?ch[0]->size:0)+(ch[1]?ch[1]->size:0)+1;
	}
	}*root;
	unsigned seed;
	inline int size(node *o){return o?o->size:0;}
	inline void rotate(node *&a,bool d){
	node *b=a;
	a=a->ch[!d];
	a->size=b->size;
	b->ch[!d]=a->ch[d];
	a->ch[d]=b;
	b->up();
	}
	inline int ran(){
	return seed=0xdefaced*seed+1;
	}
	void clear(node *&p){
	if(p)clear(p->ch[0]),clear(p->ch[1]),delete p;
	}
	void insert_as_root(node *&p,const T &data){
	if(!p)p=new node(data);
	else{
	p->size++;
	insert_as_root(p->ch[p->data<data],data);
	rotate(p,!(p->data<data));
	}
	}
	void insert(node *&p,const T &data){
	if(ran()%(size(p)+1)==0){
	insert_as_root(p,data);
	}else{
	p->size++;
	insert(p->ch[p->data<data],data);
	}
	}
	node *merge(node *a,node *b){
	if(!a||!b)return a?a:b;
	if(ran()%(a->size+b->size)<a->size){
	a->ch[1]=merge(a->ch[1],b);
	a->up();
	return a;
	}else{
	b->ch[0]=merge(a,b->ch[0]);
	b->up();
	return b;
	}
	}
	bool erase(node *&o,const T &data){
	if(!o)return 0;
	if(o->data==data){
	node *t=o;
	o=merge(o->ch[0],o->ch[1]);
	delete t;
	return 1;
	}else if(erase(o->ch[o->data<data],data)){
	o->size--;return 1;
	}return 0;
	}
	public:
	random_bst(unsigned x=20150112):root(0),seed(x){}
	~random_bst(){clear(root);}
	inline void clear(){clear(root),root=0;}
	inline void insert(const T &data){
	insert(root,data);
	}
	inline bool erase(const T &data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	for(node *o=root;o;)
	if(o->data==data)return 1;
	else o=o->ch[o->data<data];
	return 0;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o;)
	if(o->data<data)cnt+=size(o->ch[0])+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=size(o->ch[0]))o=o->ch[0];
	else if(k==size(o->ch[0])+1)return o->data;
	else k-=size(o->ch[0])+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return size(root);}
	};
	#endif

view raw Randomized Binary Search Tree1.cpp hosted with ❤ by GitHub

以下為第二種作法:

	#ifndef RANDOMIZED_BINARY_SEARCH_TREE
	#define RANDOMIZED_BINARY_SEARCH_TREE
	template<typename T>
	class random_bst{
	private:
	struct node{
	T data;
	unsigned s;
	node *ch[2];
	node(const T&d):data(d),s(1){}
	node():s(0){ch[0]=ch[1]=this;}
	}*nil,*root;
	unsigned x;
	inline unsigned ran(){
	return x=x*0xdefaced+1;
	}
	inline void rotate(node *&a,bool d){
	node *b=a;
	a=a->ch[!d];
	a->s=b->s;
	b->ch[!d]=a->ch[d];
	a->ch[d]=b;
	b->s=b->ch[0]->s+b->ch[1]->s+1;
	}
	void insert_as_root(node *&p,const T &data){
	if(!p->s){
	p=new node(data);
	p->ch[0]=p->ch[1]=nil;
	}else{
	++p->s;
	insert_as_root(p->ch[p->data<data],data);
	rotate(p,!(p->data<data));
	}
	}
	void insert(node *&p,const T &data){
	if(ran()%(p->s+1)==0){
	insert_as_root(p,data);
	}else{
	++p->s;
	insert(p->ch[p->data<data],data);
	}
	}
	bool erase(node *&o,const T &data){
	if(!o->s)return 0;
	if(o->data==data){
	if(!o->ch[0]->s||!o->ch[1]->s){
	node *t=o;
	o=o->ch[0]->s?o->ch[0]:o->ch[1];
	delete t;
	}else{
	--o->s;
	bool d=ran()%(o->ch[0]->s+o->ch[1]->s)<o->ch[0]->s;
	rotate(o,d);
	erase(o->ch[d],data);
	}
	return 1;
	}else if(erase(o->ch[o->data<data],data)){
	--o->s;return 1;
	}else return 0;
	}
	void clear(node *&o){
	if(o->s)clear(o->ch[0]),clear(o->ch[1]),delete o;
	}
	public:
	random_bst(unsigned s=20150112):nil(new node),root(nil),x(s){}
	~random_bst(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T &data){
	insert(root,data);
	}
	inline bool erase(const T &data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	for(node *o=root;o->s;)
	if(o->data==data)return 1;
	else o=o->ch[o->data<data];
	return 0;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->s;}
	};
	#endif

view raw Randomized Binary Search Tree2.cpp hosted with ❤ by GitHub

字典樹 Trie

今天完成了字典樹 Trie的模板，可以在O(N) (N=字串長度)的時間搜尋一個字串
插入的時間亦是O(N)，在字串匹配上有很大的幫助(AC字動機)，但是因為每個節點都要記錄所有字元，是非常耗空間的
通常將根節點設為"沒有字元"，若是想進一步了解Trie請自行google吧!
以下為模板:

	#ifndef TRIE
	#define TRIE
	template<typename T,char L='A',char R='z'>
	class trie{
	private:
	struct node{
	node *tr[R-L+1],*pa;
	T data;
	int ref;
	bool is;
	node():ref(0),is(0){
	for(int i=0;i<=R-L;i++)tr[i]=NULL;
	}
	~node(){
	for(int i=0;i<=R-L;i++){
	if(tr[i])delete tr[i];
	}
	}
	}root;
	int _size;
	public:
	class point_iterator{
	friend trie;
	private:
	node *p;
	public:
	point_iterator(){p=NULL;}
	inline operator bool (){return p?1:0;}
	inline T &operator *(){return p->data;}
	inline void operator =(T d){p->data=d;}
	};
	inline void clear(){
	for(int i=0;i<=R-L;i++)if(root.tr[i]){
	delete root.tr[i],root.tr[i]=NULL;
	}
	}
	trie(){root.ref=0;clear();_size=0;}
	~trie(){clear();}
	inline int size(){return _size;}
	inline point_iterator find(const char *s){
	point_iterator it;
	node *now=&root;
	for(int i=0;s[i];i++){
	if(!now->tr[s[i]-L])return it;
	now=now->tr[s[i]-L];
	}
	if(now->is)it.p=now;
	return it;
	}
	inline point_iterator insert(const char *s,T data){
	point_iterator it=find(s);
	if(it){it=data;return it;}
	node *now=&root;
	_size++;
	for(int i=0;s[i];i++){
	if(!now->tr[s[i]-L]){
	now->tr[s[i]-L]=new node;
	now->tr[s[i]-L]->pa=now;
	}
	now=now->tr[s[i]-L],now->ref++;
	}
	now->data=data;now->is=1;
	it.p=now;
	return it;
	}
	inline T &operator[](const char *s){
	point_iterator it=find(s);
	if(!it)it=insert(s,T());
	return it.p->data;
	}
	inline bool erase(const char *s){
	if(!find(s))return 0;
	node *now=&root;_size--;
	for(int i=0;s[i];i++){
	if(i&&!now->ref){
	now->pa->tr[s[i-1]-L]=NULL;
	delete now;return 1;
	}
	now=now->tr[s[i]-L],now->ref--;
	}
	now->is=0;
	return 1;
	}
	};
	#endif

view raw trie.cpp hosted with ❤ by GitHub

跟stl map的用法類似
trie<int >t;  //定義一個存int的Trie
trie<int ,'a','z'>t  //定義一個字元範圍是'a'到'z'的Trie(預設是'A'~'z')
trie<int >::point_iterator it;  //定義一個節點迭代器
t.insert("asd",1);  //插入1到"asd"的位置(若原來已有值則覆蓋)，傳回插入節點的point_iterator
t.find("asd");  //傳回"asd"位置的point_iterator，若無資料point_iterator為NULL，請注意
t["asf"];  //若"asd"存在資料，傳回其引用位置，若不存在則插入並傳回其引用位置
t.erase("asd");  //刪除"asd"節點，若成功刪除則傳回true，失敗(無結點)則傳回false
t.size(); //傳回字串的總數
t.clear();  //清空trie

[ size balanced tree ] 節點大小平衡樹 ( 傻B樹 ) 複雜度證明(論文)

找了好久終於找到中国广东中山纪念中学的陈启峰size balanced tree的原作論文中文翻譯版
若需參考要者請點擊下載
原文連結

重量左偏樹 weight-biased leftist tree

之前提到的左偏樹是深度左偏樹，雖然複雜度也是logn，但是重量左偏樹的速度在實作上是明顯較深度左偏樹快的
若是想了解深度左偏樹 height-biased leftist tree 或是左偏樹的理論請參考這篇文章
以下提供模板(使用方法與深度左偏樹相同):

	#include<functional>
	#ifndef WEIGHT_BIASED_LEFTIST_TREE
	#define WEIGHT_BIASED_LEFTIST_TREE
	template<typename T,typename _Compare=std::less<T> >
	class wb_leftist_tree{
	private:
	struct node{
	T data;
	int size;
	node *l,*r;
	node(const T&d):data(d),size(1),l(0),r(0){}
	}*root;
	_Compare cmp;
	inline int size(node *o){return o?o->size:0;}
	node *merge(node *a,node *b){
	if(!a||!b)return a?a:b;
	if(cmp(a->data,b->data)){
	node *t=a;a=b;b=t;
	}
	a->r=merge(a->r,b);
	if(size(a->l)<size(a->r)){
	node *t=a->l;a->l=a->r;a->r=t;
	}
	a->size=size(a->l)+size(a->r)+1;
	return a;
	}
	void _clear(node *&o){
	if(o)_clear(o->l),_clear(o->r),delete o;
	}
	public:
	wb_leftist_tree():root(0){}
	~wb_leftist_tree(){_clear(root);}
	inline void clear(){
	_clear(root),root=0;
	}
	inline void join(wb_leftist_tree &o){
	root=merge(root,o.root);
	o.root=0;
	}
	inline void swap(wb_leftist_tree &o){
	node *t=root;
	root=o.root;
	o.root=t;
	}
	inline void push(const T&data){
	root=merge(root,new node(data));
	}
	inline void pop(){
	if(!root)return;
	node *tmd=merge(root->l,root->r);
	delete root;
	root=tmd;
	}
	inline const T& top(){return root->data;}
	inline int size(){return size(root);}
	inline bool empty(){return !size(root);}
	};
	#endif

view raw weight_biased_leftist_tree.cpp hosted with ❤ by GitHub

在函式中傳入使用者自定函式的方法(C++ Functor 仿函式)

一般來說在函式中傳入使用者自定函式通常會用函數指標
但是C++ template提供了一個較為便利的函數參數，以物件的方式傳入
可支援一般函數或利用重載()運算子所定義的物件
詳情請看代碼:

	#include<stdio.h>
	template<typename T,typename _Compare>
	T min_element(T first,T last,_Compare cmp){
	//這種做法可以將函式當成物件傳入
	T ans=first;
	while(++first!=last){
	if(cmp(*first,*ans))ans=first;
	}
	return ans;
	}
	//盡量使用const避免出錯，一般STL的傳入函式也請這樣使用
	bool const cmp(const int &a,const int &b){
	return a<b;
	}
	struct _cmp{
	bool operator ()(const int &a,const int &b)const{
	return a>b;
	}
	};
	template<typename T>
	struct __cmp{
	bool operator ()(const T &a,const T &b)const{
	return a<b;
	}
	};
	int s[20]={1,5,6,-91,2,8,4,5,-5,-9,-4,-1,5,8,6,3,7,-9};
	//max:8 min:-91
	int main(){
	_cmp my_cmp;//創立一個比較物件
	printf("%d\n",*min_element(s,s+18,my_cmp));//8
	printf("%d\n",*min_element(s,s+18,cmp));//-91
	printf("%d\n",*min_element(s,s+18,_cmp()));//8
	printf("%d\n",*min_element(s,s+18,__cmp<int >()));//-91
	return 0;
	}

view raw Functor.cpp hosted with ❤ by GitHub

歐拉函數 乘法模逆元

今天來講一下歐拉函數及乘法模逆元
首先是歐拉函數
定義: 歐拉函數$φ(N)$是小于或等于$N$的正整數中與$N$互質的數的個數$(N>0)$
其值為
$φ(N)=N \times (1-1/P_1) \times (1-1/P_2) \times (1-1/P_3) \times ... \times (1-1/P_k)$
其中$P_i$為$N$的質因數，總共有$k個$

因此可以利用質數篩法在線性或趨近線性的時間內完成某一區間的歐拉函數:

	#define maxn 10000000
	int euler[maxn+5];
	bool is_prime[maxn+5];
	inline void init_euler(){
	is_prime[1]=1;//一不是質數
	for(int i=1;i<=maxn;i++)euler[i]=i;
	for(int i=2;i<=maxn;i++){
	if(!is_prime[i]){//是質數
	euler[i]--;
	for(int j=i<<1;j<=maxn;j+=i){
	is_prime[j]=1;
	euler[j]=euler[j]/i*(i-1);
	}
	}
	}
	}

view raw range phi.cpp hosted with ❤ by GitHub

計算單一數的歐拉函數值可利用平方根質數法求:

	inline int Euler(int n){
	int ans=n;
	for(int i=2;i*i<=n;i++){
	if(n%i==0){
	ans=ans/i*(i-1);
	while(n%i==0)n/=i;
	}
	}
	if(n>1)ans=ans/n*(n-1);
	return ans;
	}

view raw phi.cpp hosted with ❤ by GitHub

接下來介紹乘法模逆元
一整數$A$對同餘$B$之乘法模逆元是指滿足以下公式的整數$R$ $$A^{-1}≡R \; (mod \; B)$$
也可以寫成 $$AR≡1 \; (mod \; B)$$
根據歐拉定理: 當$gcd(A,B)=1$，$A^{φ(B)} ≡ 1 \; mod \;B$
(若$B$是質數，則$φ(B)=B-1$，競賽題目$B$大多為質數，不需求歐拉函數)
為了方便，我們定義 $X\%Y$ 表示$X$除以$Y$的餘數
顯然 $R=A^{φ(B)-1} \;\% \; B$ 是$A$的模$B$乘法模逆元，可以由次方快速冪再$\ord{logn}$的時間得到:

	inline long long pow_mod(int a,int b){//條件:gcd(a,b)=1
	long long mod=b,tmd=a,ans=1;
	for(b=Euler(b)-1;b;b>>=1){
	if(b&1)ans=ans*tmd%mod;
	tmd=tmd*tmd%mod;
	}
	return ans;
	}

view raw mod inverse.cpp hosted with ❤ by GitHub

另一種模逆元的解法:
計算$A$的模$B$乘法模逆元，可由擴展歐幾里得演算法得到
設$exgcd$擴展歐幾里得演算法的函數，它接受兩個整數$A,B$，輸出三個整數$g,x,y$。
$g,x,y$滿足等式$A \times x+B \times y=g$，且$g=gcd(A,B)$。
當$B=0$時，有$x=1,\; y=0$使等式成立
當$B>0$，在歐基里德算法的基礎上，已知 $$gcd(A,B)=gcd(B,A\%B)$$
先遞迴求出$x^{'},y^{'}$滿足 $$Bx^{'}+(A\%B)y^{'}=gcd(B,A\%B)=gcd(A,B)$$
然後可以將上面的式子化簡得 $$Bx^{'}+(A-(A/B) \times B)y^{'}=gcd(A,B)$$ $$Ay^{'}+Bx^{'}-(A/B) \times By^{'}=gcd(A,B)$$
這裡的除法是整除法，會自動無條件捨去。把含$B$的因式提取一個$B$，可得 $$Ay^{'}+B(x^{'}-(A/B) \times y^{'})=gcd(A,B)$$
故$x=y^{'},\; y=x^{'}-(A/B) \times y^{'}$

若$g=1$，則$A$的模$B$乘法模逆元為$B+x$
$(A \times x+B \times y)\%B=1 \longrightarrow B \times y是B的倍數 \longrightarrow A(x+B)\%B=1 \; (因為x可能小於0)$
總複雜度和歐基里德算法相同，皆為$\ord{logn}$

擴展歐幾里得演算法函數實做:

	template<typename T>
	T exgcd(T a,T b,T &x,T &y){//a*x+b*y=(return)
	if(b){
	T tmd=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return tmd;
	}
	x=1,y=0;
	return a;
	}

view raw exgcd.cpp hosted with ❤ by GitHub

[ size balanced tree ] 節點大小平衡樹 ( 傻B樹 )

它是由中国广东中山纪念中学的陈启峰發明的[已新增其中文論文連結]
跟一般的平衡樹一樣，靠左旋和右旋維持其平衡
他的精華來自於maintain函數
聽說是均攤$\ord 1$的，傳說中速度僅次於紅黑樹
coding複雜度跟Treap差不多，比AVL Tree好寫
平衡調件為: 每顆子樹的大小不小於其兄弟子樹的大小
我只提供模板，這裡有詳細的介紹

	#ifndef SIZE_BALANCED_TREE
	#define SIZE_BALANCED_TREE
	template<typename T>
	class sb_tree{
	private:
	struct node{
	node *ch[2];
	int s;
	T data;
	node(const T&d):s(1),data(d){}
	node():s(0){ch[0]=ch[1]=this;}
	}*nil,*root;
	inline void rotate(node *&a,bool d){
	node *b=a->ch[!d];
	a->ch[!d]=b->ch[d];
	b->ch[d]=a;
	b->s=a->s;
	a->s=a->ch[0]->s+a->ch[1]->s+1;
	a=b;
	}
	void maintain(node *&o,bool d){
	if(o->ch[d]->ch[d]->s>o->ch[!d]->s){
	rotate(o,!d);
	}else if(o->ch[d]->ch[!d]->s>o->ch[!d]->s){
	rotate(o->ch[d],d);
	rotate(o,!d);
	}else return;
	maintain(o->ch[1],1);
	maintain(o->ch[0],0);
	maintain(o,1);
	maintain(o,0);
	}
	void insert(node *&o,const T&data){
	if(!o->s){
	o=new node(data);
	o->ch[0]=o->ch[1]=nil;
	}else{
	++o->s;
	bool d=o->data<data;
	insert(o->ch[d],data);
	maintain(o,d);
	}
	}
	inline void success_erase(node *&o){
	node *t=o;
	o=o->ch[0]->s?o->ch[0]:o->ch[1];
	delete t;
	}
	bool erase(node *&o,const T&data){
	if(!o->s)return 0;
	if(o->data==data){
	if(!o->ch[0]->s||!o->ch[1]->s)success_erase(o);
	else{
	--o->s;
	node **t=&o->ch[1];
	for(;(*t)->ch[0]->s;t=&(*t)->ch[0])(*t)->s--;
	o->data=(*t)->data;
	success_erase(*t);
	}return 1;
	}else if(erase(o->ch[o->data<data],data)){
	--o->s;return 1;
	}else return 0;
	}
	void clear(node *&o){
	if(o->s)clear(o->ch[0]),clear(o->ch[1]),delete o;
	}
	public:
	sb_tree():nil(new node),root(nil){}
	~sb_tree(){clear(root),delete nil;}
	inline void clear(){clear(root),root=nil;}
	inline void insert(const T&data){
	insert(root,data);
	}
	inline bool erase(const T&data){
	return erase(root,data);
	}
	inline bool find(const T&data){
	node *o=root;
	while(o->s&&o->data!=data)o=o->ch[o->data<data];
	return o->s;
	}
	inline int rank(const T&data){
	int cnt=0;
	for(node *o=root;o->s;)
	if(o->data<data)cnt+=o->ch[0]->s+1,o=o->ch[1];
	else o=o->ch[0];
	return cnt;
	}
	inline const T&kth(int k){
	for(node *o=root;;)
	if(k<=o->ch[0]->s)o=o->ch[0];
	else if(k==o->ch[0]->s+1)return o->data;
	else k-=o->ch[0]->s+1,o=o->ch[1];
	}
	inline const T&operator[](int k){
	return kth(k);
	}
	inline const T&preorder(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(x->data<data)y=x,x=x->ch[1];
	else x=x->ch[0];
	if(y)return y->data;
	return data;
	}
	inline const T&successor(const T&data){
	node *x=root,*y=0;
	while(x->s)
	if(data<x->data)y=x,x=x->ch[0];
	else x=x->ch[1];
	if(y)return y->data;
	return data;
	}
	inline int size(){return root->s;}
	};
	#endif

view raw size balanced tree.cpp hosted with ❤ by GitHub

因為自己寫的速度有點慢，可以參考這邊一些大神寫的code