[ 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