[ Cyclic Longest Common Subsequence ] 環狀最長共同子序列

相信會看我發的文的人應該多少都知道LCS(最長共同子序列)是什麼吧？對於這樣的問題很早以前就有人提出$\ord{nm}$的方法了。
現在如果將問題中的字串改成頭尾相接環狀的字串，我們稱其最長共同子序列為CLCS，例如：

設A,B為頭尾相接環狀字串

A = "aabaa" , B = "aaabb"

CLCS(A,B) = "aaab"

像這樣的問題要怎麼解呢？剛好codeforces上面就有一題：http://codeforces.com/gym/100203/problem/B

為了解決這一題，查閱了不少奇怪的文章，最後給我找到了這篇論文：Solving Cyclic Longest Common Subsequence in Quadratic Time

對於CLCS提出了$\ord{nm}$的演算法，在這裡附上該演算法的實作：

	#include<cstring>
	#include<algorithm>
	using std::max;
	const int MAXN = 1505;
	enum traceType{LEFT,DIAG,UP};
	int dp[MAXN*2][MAXN], pa[MAXN*2][MAXN];
	char AA[MAXN*2];
	void LCS(const char *a, const char *b, int m, int n){
	for(int i=1; i<=m; ++i)
	for(int j=1; j<=n; ++j){
	if(a[i-1]==b[j-1]) dp[i][j]=dp[i-1][j-1]+1;
	else dp[i][j]=max(dp[i][j-1], dp[i-1][j]);
	if(dp[i][j]==dp[i][j-1]) pa[i][j]=LEFT;
	else if(a[i-1]==b[j-1]) pa[i][j]=DIAG;
	else pa[i][j]=UP;
	}
	}
	int trace(int m, int n){
	int res = 0;
	while(m&&n){
	if(pa[m][n]==LEFT) --n;
	else if(pa[m][n]==UP) --m;
	else --m, --n, ++res;
	}
	return res;
	}
	void reRoot(int root,int m, int n){
	int i=root, j=1;
	while(j<=n&&pa[i][j]!=DIAG) ++j;
	if(j>n) return;
	pa[i][j] = LEFT;
	while(i<m&&j<n){
	if(pa[i+1][j]==UP) pa[++i][j]=LEFT;
	else if(pa[i+1][j+1]==DIAG)
	pa[++i][++j]=LEFT;
	else ++j;
	}
	while(i<m&&pa[++i][j]==UP) pa[i][j]=LEFT;
	}
	int CLCS(const char *a, const char *b){
	int m=strlen(a), n=strlen(b);
	strcpy(AA,a); strcpy(AA+m,a);
	LCS(AA,b,m*2,n);
	int ans = dp[m][n];
	for(int i=1; i<m; ++i){
	reRoot(i,m*2,n);
	ans=max(ans,trace(m+i,n));
	}
	return ans;
	}

view raw reRoot.cpp hosted with ❤ by GitHub

作法大致上是將b字串變成原來的兩倍，然後做一般的LCS，但在過程中要注意更新答案的順序。接著修改DP來源表格，每次將來源表格的第一個row刪掉(reRoot函數)，然後計算當前的答案。

當我AC這一題後，我看到了一個非常驚人的作法，也是用DP複雜度是$\ord{nm}$，但是我目前還無法完全理解為何這樣是可行的，而且這樣的做法我也不知道如何回溯找出一組解，只能知道CLCS的長度，希望路過的大大們能提供想法：

	#include<cstring>
	#include<algorithm>
	const int MAXN=1505;
	int h[MAXN][MAXN*2], v[MAXN][MAXN*2];
	char B[MAXN*2];
	int CLCS(const char *a, const char *b){
	int m = strlen(a), n = strlen(b);
	strcpy(B,b); strcpy(B+n,b);
	for(int j=0; j<n*2; ++j) h[0][j] = j+1;
	for(int i=0; i<m; ++i){
	for(int j=0; j<n*2; ++j){
	if(a[i]==B[j]){
	h[i+1][j] = v[i][j];
	v[i][j+1] = h[i][j];
	}else{
	h[i+1][j] = std::max(h[i][j], v[i][j]);
	v[i][j+1] = std::min(h[i][j], v[i][j]);
	}
	}
	}
	int ans=0;
	for(int i=0; i<n; ++i){
	int sum = 0;
	for(int j=0; j<n; ++j)
	if (h[m][j+i] <= i) ++sum;
	ans = std::max(ans,sum);
	}
	return ans;
	}

view raw HV.cpp hosted with ❤ by GitHub

[ lower convex hull self-balancing binary search tree ] 下凸包平衡樹

下凸包平衡樹是一個動態維護下凸包的資料結構，他至少能在$\ord{\log n}$時間完成以下兩個操作：

1. 查詢 point(x,y) 是否被包含於下凸包中
2. 將 point(x,y) 加入至下凸包中，然後將不位於下凸包頂點上的點刪除

由於繼承自C++ STL map，因此所有能對map用的操作都能被使用，對於那些會改動map資料的操作請小心使用避免影響到凸包的正確性，例如operator[]就盡量不要使用。

對於某些動態規劃問題需要用凸包斜率優化時應該非常有用，有空會補上專門用來動態做斜率優化的版本

	#include<map>
	using std::map;
	using std::pair;
	#define x first
	#define y second
	template<typename T>
	class convexBST : public map<T,T>{
	typedef const pair<T,T>& point;
	static T cross(point a, point b, point c){
	return (a.x-c.x) * (b.y-c.y) - (b.x-c.x) * (a.y-c.y);
	}
	public:
	bool isIn(point p)const{
	if(this->empty()) return 0;
	if(this->count(p.x)) return p.y >= this->at(p.x);
	if(p.x < this->begin()->x || (--this->end())->x < p.x) return 0;
	auto l = this->lower_bound(p.x), r = l--;
	return cross(*r,p,*l)>=0;
	}
	void add(point p){
	if(isIn(p)) return;
	(*this)[p.x] = p.y;
	auto l = this->upper_bound(p.x), r = l--;
	if(r!=this->end())
	for(auto r2 = r++; r!=this->end(); r2=r++){
	if(cross(*r2,*r,p)>0) break;
	this->erase(r2);
	}
	if(l!=this->begin()&&--l!=this->begin())
	for(auto l2 = l--; l2!=this->begin(); l2=l--){
	if(cross(*l,*l2,p)>0) break;
	this->erase(l2);
	}
	}
	};
	#undef x
	#undef y

view raw convexBST.cpp hosted with ❤ by GitHub

	#include<bits/stdc++.h>
	using namespace std;
	int main(){
	convexBST<int> hull;
	vector<pair<int,int>> pointSet={{0,-1},{0,1},{-1,-5},{6,4},{-8,8}};
	for(auto p:pointSet)
	hull.add(p);
	printf("Points on the lower convex hull:");
	for(auto p:hull)
	printf(" (%d,%d)",p.first,p.second);
	puts("");
	puts(hull.isIn({6,0}) ? "Point(6,0) is in." : "Point(6,0) isn\'t in.");
	puts(hull.isIn({1,1}) ? "Point(1,1) is in." : "Point(1,1) isn\'t in.");
	return 0;
	}

view raw test.cpp hosted with ❤ by GitHub

[ simplex algorithm ] 線性規劃 - 單純形法

線性規劃簡介

線性規劃是在線性約束條件下尋找某個目標線性函數的最大、最小值，例如：
$$\begin{array}{lrl}     \min& -3x_1+2x_2 & \\     \mathrm{subject~to}         &4x_1-6x_2&\le 5\\         &5x_1-2x_2&\ge 7\\         &x_1&\ge 0.     \end{array}$$

在最佳化問題中，線性規劃是非常重要的領域，很多的演算法問題都可以轉換成線性規劃問題來解決，例如：最大網路流、最大匹配等。因此研究如何系統化的解決線性規劃問題是非常重要的事情

線性規劃標準形

所有的線性規劃問題都可以轉換成下列標準型：

$$\begin{array}{lrl}     \max &\sum_{j=1}^n c_j \times x_j &\\     \mathrm{subject~to}         &\sum_{j=1}^n A_{i,j} \times x_j &\le b_j ,\; \forall i=1\sim m\\         &x_j& \geq 0, \; \forall  j = 1\sim n     \end{array}$$

一般來說大多數解線性規劃的演算法都會要求將問題轉換成標準形，因此務必要理解其轉換方法。轉換成標準形的步驟如下：

1. 將最小化改成最大化：
   只要改變目標函數的正負號，即可將最小化問題改成標準型設定的最大化問題，也就是說 $$\min \sum_{j=1}^n c_j \times x_j$$ 等價於 $$\max \sum_{j=1}^n -c_j \times x_j$$
2. 去除等式：
   如果約束條件中存在等式，將其轉化為兩個分別為大於等於及小於等於的不等式取代，例如：
   $$x_1+x_2=5$$ 等價於 $$x_1+x_2 \le 5\\x_1+x_2 \ge 5$$
3. 去除大於等於：
   如果約束條件中存在大於等於約束，將約束兩邊取負即可，例如：
   $$x_1+x_2 \ge 5$$ 等價於 $$-x_1-x_2 \le -5$$
4. 去除自由變數：
   如果未知數 $x_i$ 沒有非負約束，我們說 $x_i$ 是自由變數。加入新變量 $x_i'$，並用 $x_i-x_i'$替換原來的變量 $x_i$，並增加$x_i,x_i' \ge 0$的條件即可

標準化範例：

$$\begin{array}{lrl}     \min& -x_1+2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&= 5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
首先將最小化改成最大化：
$$\begin{array}{lrl}     \max& x_1-2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&= 5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
接著去除等式：
$$\begin{array}{lrl}     \max& x_1-2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&\le 5\\         &x_1+x_2&\ge 5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
去除大於等於：
$$\begin{array}{lrl}     \max& x_1-2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&\le 5\\         &-x_1-x_2&\le -5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
最後去除自由變數，將$x_1$用$x_1-x_3$取代：
$$\begin{array}{lrl}     \max& x_1-2x_2-x_3 & \\     \mathrm{subject~to}         &x_1+x_2-x_3&\le 5\\         &-x_1-x_2+x_3&\le -5\\         &x_1-x_2-x_3&\le 2\\         &x_1,x_2,x_3&\ge 0.     \end{array}$$

單純形法

輸入input

本人實作的單純形法其輸入是一個被稱為單純形表的 $m+1 \times n+1$ 矩陣$A$，標準形中每個項目的值都可以對應到矩陣$A$上：
$$\begin{array}{lrl}     \max &\sum_{j=1}^n A_{0,j} \times x_j &\\     \mathrm{subject~to}         &\sum_{j=1}^n A_{i,j} \times x_j &\le A_{i,0} ,\; \forall i=1\sim m\\         &x_j& \geq 0, \; \forall  j = 1\sim n     \end{array}$$
可以發現：

1. $A_{0,1} \sim A_{0,n}$對應標準形中的$c_1 \sim c_n$
2. $A_{i,0}$對應標準形中的$b_i ,\; \forall i=1\sim m$

為了方便計算請設$A_{0,0}=0$

上面標準化範例寫成單純形表為以下形式：

$$\begin{bmatrix} 0& 1 &-2 &-1 \\ 5 &1& 1 &-1 \\ -5 &-1 &-1 &1 \\ 2 &1 &-1& -1 \end{bmatrix}$$

輸出output

輸出為一個$size=n+1$的vector $ans$，其中：

1. $ans[0]$為線性規劃的解，也就是最大值
2. $ans[1] \sim ans[n]$為變數$x_1 \sim x_n$的其中一組滿足條件的答案
3. 如果無解或是解無限大則$ans$的$size=0$
4. 注意如果答案是 $0$ 由於浮點數誤差的關係可能會產生 $``-0"$ 這種東西，要小心處理

上面標準化範例的解如下：

$$\{0.5,\; 3.5,\; 1.5,\; 0\}$$

意即在 $x_1=3.5,\;x_2=1.5,\;x_3=0$ 的情況下會有最大值 $x_1-2x_2-x_3=0.5$

程式碼

使用C++11

	#include<vector>
	#include<algorithm>
	using namespace std;
	template<class VDB>
	VDB simplex(int m,int n,vector<VDB> a){
	vector<int> left(m+1), up(n+1);
	iota(left.begin(), left.end(), n);
	iota(up.begin(), up.end(), 0);
	auto pivot = [&](int x, int y){
	swap(left[x], up[y]);
	auto k = a[x][y]; a[x][y] = 1;
	vector<int> pos;
	for(int j = 0; j <= n; ++j){
	a[x][j] /= k;
	if(a[x][j] != 0) pos.push_back(j);
	}
	for(int i = 0; i <= m; ++i){
	if(a[i][y]==0 || i == x) continue;
	k = a[i][y], a[i][y] = 0;
	for(int j : pos) a[i][j] -= k*a[x][j];
	}
	};
	for(int x,y;;){
	for(int i=x=1; i <= m; ++i)
	if(a[i][0]<a[x][0]) x = i;
	if(a[x][0]>=0) break;
	for(int j=y=1; j <= n; ++j)
	if(a[x][j]<a[x][y]) y = j;
	if(a[x][y]>=0) return VDB();
	pivot(x, y); // infeasible
	}
	for(int x,y;;){
	for(int j=y=1; j <= n; ++j)
	if(a[0][j] > a[0][y]) y = j;
	if(a[0][y]<=0) break;
	x = -1;
	for(int i=1; i<=m; ++i) if(a[i][y] > 0)
	if(x == -1 || a[i][0]/a[i][y] < a[x][0]/a[x][y])
	x = i;
	if(x == -1) return VDB();
	pivot(x, y); // unbounded
	}
	VDB ans(n + 1);
	for(int i = 1; i <= m; ++i)
	if(left[i] <= n) ans[left[i]] = a[i][0];
	ans[0] = -a[0][0];
	return ans;
	}

view raw simplex.h hosted with ❤ by GitHub

	#include<cstdio>
	#include"simplex.cpp"
	int main(){
	int m = 3, n = 3;
	vector<vector<double>> a=
	{
	{ 0, 1, -2, -1},
	{ 5, 1, 1, -1},
	{-5, -1, -1, 1},
	{ 2, 1, -1, -1}
	};
	auto res=simplex(m,n,a);
	if(res.size()==n+1){
	for(int i=0;i<=n;++i) printf("%f ",res[i]);puts("");
	}
	return 0;
	}

view raw test.cpp hosted with ❤ by GitHub

原理

有空時會補上，先讓我好好休息吧！終於可以正常睡覺了

[ Implement Queue using Stacks ] 用堆疊實作佇列

一般在C++ STL中queue都是用deque實作的，而deque用和vector類似的倍增法來實作(稍微多了一些操作導致常數有點大但是每次操作時間較為平均)，雖然夠用但有些時候這不是一個好選擇。
大二資料結構期中考時，考了一題叫你用stack實作queue的方法，我那時候想了一下想了一個超帥的作法，需要用到兩個stack，我把它成為stack A和B。

• push:
  push時直接把資料push進B中，$\ord{1}$
• pop:
  這是個問題，我們把A的top當成是queue的front，所以直接把A pop就可以了，如果A是empty的話，就把B的資料依序pop然後push進A中，在進行該操作，均攤$\ord{1}$
• front:
  如pop所述，把A的top當成是queue的front，如果A是empty再另做處理，$\ord{1}$
• back:
  和front相反，把B的top當成是queue的back，如果B是empty再另做處理，$\ord{1}$

以下為實作過程，以c++ vector代替stack：

	#include<vector>
	template<typename T>
	class st_queue{
	std::vector<T> A,B;
	public:
	void swap(st_queue& stq){
	A.swap(stq.A);
	B.swap(stq.B);
	}
	bool empty()const{
	return A.empty()&&B.empty();
	}
	size_t size()const{
	return A.size()+B.size();
	}
	T &front(){
	return A.empty() ? B.front() : A.back();
	}
	T &back(){
	return B.empty() ? A.front() : B.back();
	}
	void push(const T& data){
	B.push_back(data);
	}
	void pop(){
	if(A.size()) A.pop_back();
	else{
	for(int i=int(B.size())-1;i>0;--i)
	A.push_back(B[i]);
	B.clear();
	}
	}
	};

view raw st_queue.cpp hosted with ❤ by GitHub