后缀自动机(sam)上的字符串匹配

====

我们把相对较短的模式串构造成sam。

对于P=”abcabcacab”, T[1..i]的后缀，使得它是sam的最长前缀长度:

T: b a b c b a b c a b c a a b c a b c a b c a c a b  c

   1 1 2 3 1 1 2 3 4 5 6 7 1 2 3 4 5 6 7 5 6 7 8 9 10 4

如果最长前缀长度是|P|，则表示T[1..i]的后缀和P匹配。

内存使用

可能多个trans指针同一个节点，因此像删除树那样会引起double-free:

为此我们暂时采用内存池的做法。

如果扩展到包括数字和空格，则需要表示37个转移指针。

KMP算法

====

给定模式串P,文本串T,

假设在s位置已匹配了q个字符, 即P[1,..,q]=T[s+1,..,s+q], 而在P[q+1]不匹配。

strstr()这时会把指针指向s+2,从P[1]重新开始匹配。

当时Knuth,Morris,Pratt就想可不可以把指针再移远一点。

假设有P[1,..,k]=T[s+q+1-k,..,s+q],这时从P[k+1]开始比就行了，显然我们希望k越大越好，对应地指针移动增量=q-k越小，因此应该不会错过某些完全匹配的位置。

我们把上面两个等式合并，得到P[1,..,k]是P[1,..,q]的后缀。问题变成：

对于每个q, 求P[1,..,q]的最长的真前缀（长度记为k)，同时它也是P[1,..,q]的后缀。

我们定义前缀函数pi(q):=k.

如何计算pi(q) ?

====

使用递推的想法，假设我们已经计算好了pi(q)=k。

如果P[k+1] = P[q+1], 则显然有pi(q+1) = k+1；

否则，看作是一个匹配问题, 我们来看pi(q)的含义是与P[1,..,q]末尾匹配的最长前缀长度k，我们就拿这个前缀来匹配，并期望P[k+1]和P[q+1]一样，否则k=pi(k)循环下去。

初始条件:pi(1) = 0， 因为最长真前缀是空串。

当前P[1,..,k]匹配T[q-k+1,q]，而在T[q+1]不匹配，

应用前缀函数的定义，应该从位置s+1-k + q-k = 

一个字符串P[1,..,j]去匹配P[q+1-j,..q+1]的过程。

k=pi(k), 直到P[k] = P[q+1]。

如何做线性的字符串匹配？

====

参照后缀自动机的做法， 我们把pi和P组成一个自动机，T在这个自动机上

走一遍。

前缀函数练习题目：

1. P 在T中的出现次数? 提示：检查pi(PT)

2. (ab)^3 = ababab, 如何求最大的重复因子r=3?

3. 如何在线性时间内判断是否为循环移位，比如arc和car。（这个我还不知道怎么做）

#include 
#include 
#include 

typedef struct State State;
struct State{
	int len;
	State *suffix;
	State *trans[26];
};

static State mem_pool[1000];
static int mem_cnt = 0;

static State *NewState(void)
{
	State *p = &mem_pool[mem_cnt++];
	memset(p, 0, sizeof(*p));
	return p;
}

static void DeleteState(State *p)
{
}

static State* add(State *root, State *last, int x)
{
	State *q, *nq, *p = last, *np = NewState();
	np->len = last->len + 1;
	//printf("np %p for %c\n", np, 'a'+x);
	for(p = last; p && !p->trans[x]; p = p->suffix)p->trans[x] = np;
	if(!p)np->suffix = root;
	else{
		q = p->trans[x];
		if(q->len == p->len + 1)np->suffix = q;
		else{
			nq = NewState();
			//printf("nq %p\n", nq);
			*nq = *q;
			nq->len = p->len + 1;
			q->suffix = np->suffix = nq;
			for(; p && (p->trans[x] == q); p = p->suffix)p->trans[x] = nq;
		}
	}
	last = np;
	return last;
}

State* sam_new(char *str)
{
	State *last = NULL, *root = NULL;
	last = root = NewState();
	for(; *str; ++str)last = add(root, last, *str-'a');
	return root;
}

void sam_delete(State *root)
{
	mem_cnt = (root - mem_pool)/sizeof(State);
}

int sam_match(State *sam, int patlen, char *text)
{
	State *curr = sam;
	int x, len = 0;
	char *ptr = text;
	for(; *ptr; ++ptr){
		x = *ptr-'a';
		while(curr->suffix && !curr->trans[x])
			curr = curr->suffix, len = curr->len;
		if(curr->trans[x]){
			len++;
			curr = curr->trans[x];
		}else curr = sam, len = 0;
		//printf("%d ", len);
		if(len == patlen)printf("match begin at %ld\n", ptr-text-patlen+1);
	}
}

int main(int ac, char **av)
{
	char *pat = "aabbab";
	char *text ="xaabbabdxaabbabdx"; 
	State *sam = sam_new(pat);
	
	sam_match(sam, strlen(pat), text);
	sam_delete(sam);
}

KMP比SAM节省内存：

/*
给定模式串P,文本串T,
假设在s位置已匹配了q个字符, 即P[1,..,q]=T[s+1,..,s+q], 而在P[q+1]不匹配。
strstr()这时会把指针指向s+2,从P[1]重新开始匹配。
当时Knuth,Morris,Pratt就想可不可以把指针再移远一点。
假设有P[1,..,k]=T[s+q+1-k,..,s+q],这时从P[k+1]开始比就行了，显然我们希望k越大越好，对应地指针移动增量=q-k越小，因此应该不会错过某些完全匹配的位置。
我们把上面两个等式合并，得到P[1,..,k]是P[1,..,q]的后缀。问题变成：
对于每个q, 求P[1,..,q]的最长的真前缀（长度记为k)，同时它也是P[1,..,q]的后缀。
我们定义前缀函数pi(q):=k.

如何计算pi(q) ?
====
使用递推的想法，假设我们已经计算好了pi(q)=k。
如果P[k+1] = P[q+1], 则显然有pi(q+1) = k+1；
否则，看作是一个匹配问题, 我们来看pi(q)的含义是与P[1,..,q]末尾匹配的最长前缀长度k，我们就拿这个前缀来匹配，并期望P[k+1]和P[q+1]一样，否则k=pi(k)循环下去。
初始条件:pi(1) = 0， 因为最长真前缀是空串。

当前P[1,..,k]匹配T[q-k+1,q]，而在T[q+1]不匹配，
应用前缀函数的定义，应该从位置s+1-k + q-k = 
一个字符串P[1,..,j]去匹配P[q+1-j,..q+1]的过程。
k=pi(k), 直到P[k] = P[q+1]。

如何做线性的字符串匹配？
====
参照后缀自动机的做法， 我们把pi和P作为一个自动机，T在这个自动机上
走一遍。

前缀函数练习题目：
1. P 在T中的出现次数? 提示：检查pi(PT)
2. (ab)^3 = ababab, 如何求最大的重复因子r=3?
*/
#include 
#include 
#include 


void calc_prefix(int *pi, char *pat, int len)
{
	int q, k, ans;

	pi[0] = -1; //因为可能有多个0,所以需要这样的标记。
	pi[1] = 0;
	for(q = 2; q <= len; ++q){
		ans = 0;
		for(k = pi[q-1]; k >= 0; k = pi[k]){
			if(pat[k] == pat[q-1]){ans = k + 1; break;}
		}
		pi[q] = ans;
	}

	for(q = 1; q <= len; ++q)printf(" %d", pi[q]);
	printf("\n");
}	
void match(int *pi, char *pat, int m, char *text, int n)
{
	int q, k, ans;
	k = 1;
	for(q = 0; q < n; ++q){
		for(; k >= 0; k = pi[k]){
			if(k > 0 && pat[k-1] == text[q]){k++; break;}	
		}
		if(k <= 0)k = 1;
		if(k-1 == m){
			printf("match start at %d\n", q-m+1);
			k = 1;
		}
	}
}	

int main(int ac, char **av)
{
	int k, q;
	char *pat, *text;
	pat = "abcabcacab"; 
	text = "babcbabcabcaabcabcabcacabc";

	//pat = "ababaca"; 
	//text = "bacbababacadababacae"; 
	int len = strlen(pat);
	int *pi = malloc((len+1) * sizeof(int));

	calc_prefix(pi, pat, len);
	match(pi, pat, len, text, strlen(text));
	
	free(pi);
}

/*
$ gcc -g kmp.c && ./a.out 
 0 0 0 1 2 3 4 0 1 2
match start at 16
*/