解构公共子串及K阶马尔可夫随机文本

最新推荐文章于 2024-03-26 09:59:44 发布

随心随意随缘

最新推荐文章于 2024-03-26 09:59:44 发布

阅读量9.3k

点赞数 4

分类专栏：人工智能与机器学习算法文章标签： dictionary 算法数据结构 null qt

本文链接：https://blog.csdn.net/jn1158359135/article/details/7228705

版权

算法同时被 2 个专栏收录

7 篇文章 0 订阅

订阅专栏

人工智能与机器学习

1 篇文章 0 订阅

订阅专栏

搜索是一个常谈常新的话题，不仅有对数值型数据的搜索，更多的时候还有对字符类型的数据进行搜索的要求。比如常用的各大搜索引擎都是通过获取用户所输入的字符型数据，之后根据引擎自身所采用的算法将输入与网页数据库中存储的信息匹配，最后将这些信息对应的网页链接返回给用户，本文介绍内容之一的公共子串匹配算法就属于文本匹配算法之一。另外计算机科学在分子生物学中一个重要的应用领域即为比对(alignment)，打个不太恰当的比方，可以把氨基酸看成是字母，蛋白质是一个长句，而结构域就像单词，即具有特别意义、频繁出现在不同句子中的一串字母，结构域之所以令人感兴趣，是因为它们表征了序列内部（蛋白质是氨基酸序列）的结构或功能要素，通过采用公共子串匹配算法我们可以找出最频繁出现的一类结构域，推动分子生物学的研究。而马尔科夫随机文本算法则偏向于人工智能研究领域，在现实生活中的语音识别和数据压缩中随处可见他们的踪影。之所以将这两种算法写在一起，是因为考虑到它们的实现都利用了同一种数据结构，这里暂且留个悬念，我会放在后文揭秘这种神奇的数据结构。另外笔者在思考这篇文章的时候，根据马尔科夫随机文本算法的思想，初步构思出了一个能够实现人机对话的智能机器人系统，如果读者有兴趣，看完本文后不妨自己实现一下 :-)。

解构公共子串
我们首先从单个序列的简单例子入手，建立起一个比较直观的概念：
设序列S为"government of the people, by the people, for the people" ——引自《亚伯拉罕•林肯文集》
我们一眼便能分辨出多次出现（即至少出现两次）的最长公共子串为“ the people,”。然而如果某个序列很长，那么串匹配问题根本不可能通过人工完成，因此我们需要从上述实例中提取出一种方法，使得串匹配的过程能够通过机器自动完成。
以序列S为例放大上述“分辨”过程，我们得到一种最朴素的方法，如下所述：

定义序列的首字符g为当前字符
从当前字符的下一字符开始顺序寻找与当前字符相同的字符
一旦找到，我们定义该位置处的元素为匹配字符，随后将匹配字符与当前字符之后的元素项一一进行匹配直至某一对元素不等，我们记录匹配字符的位置为location，以及匹配子串的长度记为length，接着从匹配字符的下一字符继续顺序寻找与当前字符相同的字符
重复过程3直至序列的最后一个字符，在每一次得到新的匹配子串的长度length'后，判断length'是否大于length，以确定是否需要更新location以及length的值
将当前字符更新为下一个字符，重复执行过程2、3、4，直至当前字符已被更新为序列中的最后一个字符

根据上述描述，得到如下查找最长公共子串的过程：

getLongestSubString(s)  /*the substring appears more than once in the sequence 's'*/
location←0  /*the start of the substring*/
length←0  /*the length of the substring*/
for current←0 to length[s]-2  /*iterate until to the precursor of the last character*/
  do for match←current+1 to length[s]-1
       do if s[match]=s[current]
            then temp←0
                 repeat temp←temp+1
                   /*pay attention that must detect whether overstep the boundary of the string or not*/
                   until s[match+temp]≠NULL && s[match+temp]=s[current+temp]
                 if temp>length  /*update the location and length of the substring*/
                   then length←temp
                        location←match
return length & location

然而由于每次从当前字符开始查找子串时都需要遍历一次余下的整个序列，这使得串匹配过程的运行时间与n²渐进等同。而之所以将其放在前面介绍，是因为遇到上述字符串匹配问题时，这是符合我们直觉的解决方案。
对其优化之前，让我们更深入的思考一下上述方法：
设序列S=A₀A₁A₂A₃A₄A₅A₆A₇A₈A₉，其中A₀=A₃=A₇，A₂=A₆，A₅=A₉
上述串匹配过程可以使用下图形象化的表示

图1

由上图，我们发现字符的比对过程只在4个地方产生了暂停，然而我们将当前字符与余下序列中的字符的比对次数总共达到了(9+8+…+1)=45次。是否存在一种产生较低开销的预处理方法，使得我们可以立刻知道，所有与当前字符相同的字符各自在序列中的位置。
一种比较容易想到的方法预先是准备一组空“桶”，而后遍历整个字符串序列，最终所有相同字符在原序列中的下标以链表的形式链接在同一个“桶”中，但这种方式存在一个问题，即为序列中的当前字符寻找合适的“桶”时，期望的比对次数往往与已装入字符的桶的个数成正比，造成这种现象的原因是序列中的字符与待存入的桶并未形成一一对应的关系，而是基于一种先到先服务(First Come First Serve, FCFS)的原则。这个问题本来很容易解决，因为我们只需要像前文中所介绍的散列一样定义字符与下标之间的关系即可，一种可行的简单定义方式如下：
设序列中的各个不同元素与“桶”的下标的对应关系分别为
<(A₀, 0), (A₁, 1), (A₂, 2), (A₄, 3), (A₅, 4), (A₈, 5)>
根据上述定义方式有：

图2

预处理完成之后，我们找出某个“桶”中的所有节点，其后在原序列中进行匹配，最后得到最长子串的长度及起始位置。这种方法带来的好处是减少了不必要的一级比对过程，但因为常用字符通常占整个的序列绝大部分，这一事实决定了字符的分配必然不太均衡，从而有些桶中的链表很长，而很多桶中的节点寥寥无几甚至根本没有节点，因此这种做法不仅浪费了内存空间，并且在比对过程中造成了额外的运行时间开销，那么是否存在更好的解决方案？

如果更仔细思考上述方案，我们发现将所有相同字符链接在一起的做法相当于“聚集”的过程，然而如果“聚集”的过程不是以装”桶“的方式进行，而是以排序的方式进行呢？这是一种很有趣的思维跳跃，因为其实原先的装”桶“过程我们可以看成是区间排序，虽然已经形成基本有序，但仍不够精确，如果我们使用更加精确的排序方式来完成”聚集“过程，会有什么意想不到的效果呢？仍从一个简单的例子入手：
设序列S为banana，该序列的索引分别为index[0]~index[5]
该序列最终得到如下形成字典序的结果：

index[5]: a
index[3]: ana
index[1]: anana
index[0]: banana
index[4]: na
index[2]: nana

分析上述结果，查找最长子串的方法为：从index[5]开始直至index[2]顺序比较相邻之间的元素，根据当前比较所得结果判断是否更新最长子串的起始位置及长度。比如将元素index[5]与index[3]比较，得到最长公共子串“a”，长度为1，起始位置为5或3，之后将index[3]与index[1]比较，得到最长公共子串“ana”，长度为3，起始位置为3或1，由于当前子串的长度大于先前的子串，因此执行更新操作。依次执行上述操作，最终得到整个序列的最长公共子串即为“ana”。有意思的是，我们并不需要像先前的装“桶”策略一样，将index[5]与index[1]也执行比较操作，因为形成字典序这一精确的“聚集”过程替我们完成了隐式的比较操作。

然而有人可能会发现，这种策略其实仍然存在额外的不必要的比较操作，比如index[1]与index[0]以及index[0]与index[4]之间的比较，我们是否可以通过将已形成字典序的子串对应的索引按序装入桶中以形成有序链表？其实完全没有必要再这样做，因为这种做法占用了更多的内存空间，同时使得整个结构更加复杂，更重要的是，一个序列的绝大多数元素都是常用字符，按照这种做法最后的结果只是将原先的数组“切”成了相对较短的链表而已。基于上述分析最终的实现如下：

getLongestSubString'(s)  /*by forming dictionary order*/
/*create the array where elements are the indexes of sequence*/
/*to which we give the name 'suffix'*/
for i←0 to length[s]-1
  do suffix[i]←i
sort the 'suffix' by dictionary order
length←0
location←suffix[0]
for i←0 to length[suffix]-2
  do temp←-1
     repeat temp←temp+1
       until s[suffix[i]+temp]≠NULL && s[suffix[i]+temp]=s[suffix[i+1]+temp]
     if temp>length
       then length←temp
            location←suffix[i]
return length & location

与相比最初的实现，我们借助一个额外的数组以及排序操作将运行时间从O(n²)降低至O(n㏒n+n)。其中所使用的数组即为上文所提到的“神奇的数据结构”，它有一个专用的名称——“后缀数组(Suffix Array)”。该术语于上世纪90年代提出，但早在70年代就已经被使用了。

上文讨论了多次出现的最长公共子串，但对于“多次”实际上要求得相当宽松，以至我们的实现仅仅只是考虑了序列中重复出现两次的最长公共子串。但有时会有更加苛刻的要求，比如找出一个序列中重复出现超过N次的最长公共子串，如果我们以最原始的方法来实现恐怕将使复杂度以几何倍数上升，并且随着序列的增长，运行时间将越来越难以忍受，然而如果使用形成了字典序的后缀数组，那么这种复杂度将直接被该数据结构的序性质所规避。要理解这句话的含义，我们仍以上述序列“banana”为例，其形成字典序后的结果如下：

index[5]: a
index[3]: ana
index[1]: anana
index[0]: banana
index[4]: na
index[2]: nana

假定需要查找重复出现超过2次的最长子串，那么最直观的方法如下：
首先将index[5]与index[3]比较，我们得到重复出现两次的子串为“a”，该子串做为一个临时对象temp，继续与index[1]进行比较，我们得到重复出现3次的子串为“a”，接着我们从index[3]继续重复执行上述操作，直至index[0]处为止，依次执行上述过程最终得到出现超过2次的最长子串为“a”。仔细思考上述过程，我们发现，出现超过N次的公共子串实际就是原字典序中N+1个连续元素的交集，求解的过程可以形式化的表示为( ( (index[0] ∩ index[1]) ∩ index[2]) …∩ index[N])。
如果深入的观察一下，就会发现字典序其实使原序列形成了“梯度”有序，可以抽象的表示为下图：

图3

正如图中的旁注所示，形成字典序后，要查找出现超过N次的子串只需将索引1与索引N+1处的元素比较即可，更加严谨的证明过程如下：

设索引1~N+1所对应的序列分别为S₁,S₂,S₃,…S_N+1，其中S_N=A_N1A_N2A_N3……
因为S_x<S_y<S_z表示形成字典序，S_x与Sy必然在某个对应元素上产生不等，假定该元素的位置为i，同样的，S_y与S_z产生不等的位置假定为j
因此我们得到A_x1A_x2…A_x(i-1)=A_y1A_y2…A_y(i-1)且A_y1A_y2…A_y(j-1)=A_z1A_z2…A_z(j-1)，A_xi<A_yi，A_yj<A_zj
分两种情况讨论：①i≤j；②i>j
--->① 当i≤j时
设S_x=αA_xiβA_xjγ，S₂=α'A_yiβ'A_yjγ'，S₃=α''A_ziβ''A_zjγ''，其中α=α'，α'A_yiβ'=α''A_ziβ''，因此α=α'=α''，A_xi<A_yi=A_zi
则((S_x∩S_y)∩S_z) = ((αA_xiβA_xjγ∩α'A_yiβ'A_yjγ')∩α''A_ziβ''A_zjγ'') = α∩α''A_ziβ''A_zjγ'' = α
而S_x∩S_z = (αA_xiβA_xjγ∩α''A_ziβ''A_zjγ'') = α
所以((S_x∩S_y)∩S_z)=S_x∩S_z 得证
--->② 当i>j时同理可证
反复运用上述方法递归求解如下表达式可得
((((S₁∩S₂)∩S₃)∩S₄)……∩S_N+1)
=(((S₁∩S₃)∩S₄)……∩S_N+1)
=((S₁∩S₄)……∩S_N+1)
=((S₁∩S_N)∩S_N+1)
=S₁∩S_N+1 <--------等式(*)，注意后文的多个序列中将用到该等式
因此我们只需将索引1与N+1处的元素比较即可得到重复出现超过N次的公共子串，结论得证

根据上述分析，最终出现超过N次的最长公共子串的查找实现如下：

getLSSMultiOccur(s, N)  /*occure more than N times*/
/*assume that we have created the 'suffix' array and formed the dictionary order*/
length←0
location←suffix[0]
for i←0 to length[s]-N
  do temp←-1
     repeat temp←temp+1
       until s[suffix[i]+temp]≠NULL && s[suffix[i]+temp]=s[suffix[i+N]+temp]
     if temp>length
       then length←temp
            location←suffix[i]
return length & location

我们发现上述实现仅仅只是改变了过程getLongestSubString'(s)中的某一处，即由原先的s[suffix[i]+temp]=s[suffix[i+1]+temp]变为s[suffix[i]+temp]=s[suffix[i+N]+temp]。这也就是我在上文中所说的，“查找重复出现超过N次的子串所带来的复杂度直接被该数据结构的序性质所规避”的原因所在。

上文都是从单个序列入手，查找出现超过一定次数的最长子串，然而更多时候我们还有在多个序列中查找最长公共子串的要求。
其实在多个序列中查找最长公共子串的解决方案与单个序列相同，因为如果我们将多个序列顺序拼接成单个序列即从逻辑上将该问题转换为了前文所描述的情况。接着我们便可通过建立后缀数组并对其进行排序以找出最长公共子串。然而这样做可能存在的问题是：有可能找出的最长公共子串只是在单个序列或部分序列中重复出现。一种可行的办法首先在每个序列的尾部建立一个标记，使得将后缀数组中的元素与每个标记比较即可知道该元素所对应的序列，更形象的表示如下：

图4

从逻辑上构建完该序列之后，在实际的实现过程中，我们并不将“逻辑”序列中的标记放入后缀数组中，而是另外创建一个独立的“标记”数组。在后缀数组中的元素能够通过标记知道自身所在的原序列之后，以上问题就已经解决了一半，剩下的问题是如何在已经形成字典序的后缀数组中查找符合要求的最长公共子串。

因为在排序操作后，连续的索引可能仅仅只对应单个或部分序列。但考虑到没法形式化的定义这种序列，因此不失一般性地，可以选择一种比较简单的情况来找出一种能够得到最长公共子串的策略，并且这种策略对于任意其他的情况也是成立的。
设可能存在的索引顺序为 R₁R₂…R_iS₁S₂…S_jT₁T₂…T_k ，其中 i, j, k=1, 2, 3,……

形成字典序后，由上文所得等式(*)，((S₁∩S₂)…∩S_N+1)=(S₁∩S_N+1) 包含于 (S_N∩S_N+1)
因此 (((R₁∩R₂)…∩R_i)∩S₁)=(R₁∩S₁) 包含于 (R_i∩S₁) ，同理 (((S₁∩S₂)…∩S_j)T₁)=(S₁∩T₁) 包含于 (S_j∩T₁)
又因为部分索引顺序 R₁R₂…R_iS₁ 所对应的序列R和S的任意公共子串为 (R_x∩S₁)=(((R_x∩R_x+1)∩R_x+2)…∩S₁) 包含于 (R_i∩S₁)   --->结论(1)
同理 S_jT₁T₂…T_k 的任意公共子串 (S_j∩T_z) 包含于 (S_j∩T₁)   --->结论(2)
设上述索引顺序 R₁R₂…R_iS₁S₂…S_jT₁T₂…T_k 包含的任意公共子串为 ((R_x∩S_y)∩T_z)，则
((R_x∩S_y)∩T_z)
=(R_x∩T_z) <---注意由字典序性质可得
=(((R_x∩T_z)∩S₁)∩S_j)    <---将集合扩展
=((R_x∩S₁)∩(S_j∩T_z))
因为由结论(1)和结论(2)可知，上述索引顺序 R₁R₂…R_iS₁S₂…S_jT₁T₂…T_k 的任意公共子串 ((R_x∩S_y)∩T_z) 均包含于 ((R_i∩S₁)∩(S_j∩T₁))=(R_i∩T₁)
因此由上述分析，我们即可得到查找N个原序列的公共子串的策略为：
    在后缀数组中寻找一个最小区间[x, y]，使得该区间中至少存在一个任意原序列的索引，之后将Suffix[x]与Suffix[y]比较即可得到当前公共子串，依次顺序执行上述操作，并根据新的公共子串判断是否执行更新，最终即得到N个原序列中出现的最长公共子串。

由于在后缀数组中寻找最小区间的策略相对复杂，所以在下面的实现中使用中文注释详细讲述我的实现方法：

getLssAmongMultiSeq(s[0], s[1], … s[n-1])  /*get the longest common substring among s[0]~s[n-1]*/
/*first calculate the sum of each string's length*/
/*then we merge each string to generate a whole string named 'S'*/
/*create the suffix array named 'suffix' which contains elements as many as the sum*/
/*create the array named 'sign' containing n elements which is used to store the sign*/
/*initialize all the elements in the array 'sign' with zeor*/
/*create a boolean array containing n elements and all of them were initialized with false*/
/*the boolean array was named 'isFind'*/
/*create an array containing n elements which was named 'record'*/
/*首先初始化后缀数组suffix以及用于存放标记的数组sign*/
currlength←0
for i←0 to n-1
  do for j←0 to length[s[i]]-1
       do suffix[currlength+j]←currlength+j+i  /*the index of sign needn't to put into the 'suffix'*/
     currlength←currlength+length[s[i]]
     sign[i]←currlength
/*将后缀数组排序最终形成字典序*/
sort the array 'suffix' with dictionary order
ix←0  /*current element in 'suffix'*/
length←0
location←0
cnt←0
while TRUE
  do for i←ix to length[suffix]-1
       do k←-1  /*the kth sequence*/
          repeat k←k+1
            until suffix[ix]>sign[k]  /*根据标记判断当前后缀数组中的元素属于第几个序列*/
          if isFind[k]=false  /*说明当前查找过程中还未出现第k个序列的索引*/
            then isFind[k]←true
                 cnt←cnt+1  /*记录当前探查区间中对应的原序列的个数*/
          record[k]←suffix[ix]  /*不断更新某个原序列对应的索引，从而使得区间最小化*/
          if cnt=n  /*当探查区间中已有n个原序列时，则退出循环，完成区间的查找*/
            then break
     if i=length[suffix]  /*在最后一次查找区间时，有可能是因为所有元素均已查找完，因此我们直接退出外层循环*/
       then break
     /*排序是为了查找最小元素及最大元素，从而进行比较，另外重新查找新区间时也需要将最大元素的下一个元素赋值给ix*/
     sort the array 'record'
     k←-1
     repeat k←k+1
       until record[0]>sign[k]  /*找出最小元素所属序列*/
     isFind[k]←false  /*将该位清零，表示在新的区间中不存在该序列的索引*/
     cnt←cnt-1  /*isFind[k]清零的同时将表示序列个数的cnt减1*/
     temp←-1
     repeat temp←temp+1
       until S[record[0]+temp]≠NULL && S[record[0]+temp]=S[record[n-1]+temp]
     if temp>length
       then length←temp
            location←record[0]
     /*将最大元素的下一个元素赋值给区间起始点，以重新开始查找新的满足条件的最小区间*/
     ix←record[n-1]+1  /*find new window from next index*/
     record[0]←record[1]  /*the start of new window*/
return length & location

可以发现在上述实现中，需要执行两种不同的排序，一次是对字符串的排序操作，另一次则是对整型数组的排序操作，一种比较偷懒的做法是直接调用库函数，然而在之前的文章中可以发现，库函数的运行时间并不理想，因此我在本文中为之前所写的健壮的快速排序设计了与库函数相同的通用接口，使其能够运用于不同类型的排序，并且效率更高，以下是通用快速排序的C语言实现：

typedef int (_cdecl *comp)(const void *, const void *);
extern void QuickSort(void *, size_t, size_t, comp);
static void RandomDoubDireSort(void *, int, int, size_t, comp);
static size_t RandomDoubDirePart(void *, size_t, size_t, size_t, comp);
static void InsertionSort(void *, size_t, size_t, comp);
static size_t GeneFrom(size_t, size_t); /*return the value belongs to [left, right]*/
static void swap(char *, char *, size_t);

void QuickSort(void *src, size_t NumOfElem, size_t NumOfBytes, comp f)
{
    if(!src || NumOfBytes == 0)
        return ;
    RandomDoubDireSort(src, 0,NumOfElem-1, NumOfBytes, f);
    InsertionSort(src, NumOfElem, NumOfBytes, f);
}

static void RandomDoubDireSort(
    void *src, 
    int low,  /*pay attention the type must be integer*/
    int high,/*otherwise the arithmetic (high-low) may be overflow*/
    size_t NumOfBytes,
    comp f)
{
    while(high-low >= 49)
    {
        size_t m_Sepa = RandomDoubDirePart(src, low, high, NumOfBytes, f);
        RandomDoubDireSort(src, low, m_Sepa, NumOfBytes, f);
        low = m_Sepa+1;
    }
}

static size_t RandomDoubDirePart(
    void *src, 
    size_t low, 
    size_t high,
    size_t NumOfBytes,
    comp f)
{
    size_t temp = GeneFrom(low, high);
    char *pivot = (char *)src+temp*NumOfBytes;
    int left = low-1, right = high+1;
    while(true)
    {
        do left++;
        while(f((char *)src+left*NumOfBytes, pivot) < 0);
        do right--;
        while(f((char *)src+right*NumOfBytes, pivot) > 0);
        if(right > left)
            swap((char *)src+left*NumOfBytes, 
            (char *)src+right*NumOfBytes, NumOfBytes);
        else
            return right;
    }
}

static void InsertionSort(
    void *src, 
    size_t len,
    size_t NumOfBytes,
    comp f)
{
    char *buf = (char *)malloc(NumOfBytes);
    for(size_t index = 1; index != len; index++)
    {
        memcpy(buf, (char *)src+index*NumOfBytes, NumOfBytes);
        int nestix = index-1;
        for(; nestix != -1; nestix--)
        {
            if(f((char *)src+nestix*NumOfBytes, buf) > 0)
                memcpy((char *)src+(nestix+1)*NumOfBytes, 
                (char *)src+nestix*NumOfBytes, NumOfBytes);
            else
                break;
        }
        memcpy((char *)src+(nestix+1)*NumOfBytes,
            buf, NumOfBytes);
    }
    free(buf);
}

static void swap(char *str1, char *str2, size_t NumOfBytes)
{
    char temp;
    if(str1 == str2)
        return;
    while(NumOfBytes--) {
        temp = *str1;
        *str1++ = *str2;
        *str2++ = temp;
    }
}

static size_t GeneFrom(size_t left, size_t right)
{
    return (rand()%(right-left+1))+left;
}

因为考虑到查找多个序列的公共子串算法比较繁琐，因此直接给出使用C的可变形参函数实现的版本，以方便复用。由于在具体实现的时候，往往需要借助语言本身提供的特性，因此可能和上述伪码有些区别，不过思想是一样。

/*注意可变形参需要包含 <stdarg.h> 头文件*/
/*该结构体包含指向一个字符串的指针，及其长度*/
typedef struct {
    char *LocaOfStr;
    int LengthOfLss;
} logStr;

/*该结构包含拼接形成的整个序列ls，即后缀数组suffix，标记数组sign等*/
/*NumOfSeq表示原序列的个数*/
typedef struct StrucForMultiSeq{
    logStr ls;
    char **suffix;
    char **sign;
    size_t *record;
    bool *isFind;
    size_t NumOfSeq;
} *pStrucForMultiSeq;

/*获得公共子串的主过程，该过程返回一个存储在堆上的字符串*/
/*因此需要在主调函数中释放堆内存*/
extern char *getLssAmongMultiSeq(size_t, ... );
/*将多个原序列拼接为整个序列的子过程*/
static logStr MergeMultiSeq(va_list *, size_t);
/*在拼接以及形成字典序后查找多个序列中的公共子串*/
static logStr FindLssAfterMergeAndSort(pStrucForMultiSeq);
/*分配资源及初始化操作*/
static void AllocResAndInit(pStrucForMultiSeq);
/*释放资源*/
static void FreeRes(pStrucForMultiSeq);

/*CALLBACK function*/
int _cdecl s_comp(const void *a, const void *b)  /*sort for C-Style string*/
{ return strcmp(*(char **)a, *(char **)b); }
int _cdecl i_comp(const void *a, const void *b)  /*sort for integer type*/
{ return *(int *)a-*(int *)b; }

/*类似于库函数strdup，返回的指针指向一个堆内存，
因此注意需要手动释放内存资源以免发生泄漏*/
char *getLssAmongMultiSeq(size_t num, ...)
{
    StrucForMultiSeq  s_FMS;
    logStr ls;
    char *ret;
    s_FMS.NumOfSeq = num;
    va_list ap;
    va_start(ap, num);
    s_FMS.ls = MergeMultiSeq(&ap, num);
    AllocResAndInit(&s_FMS);
    QuickSort(s_FMS.suffix, s_FMS.ls.LengthOfLss-s_FMS.NumOfSeq,
        sizeof(s_FMS.suffix[0]), s_comp);
    ls = FindLssAfterMergeAndSort(&s_FMS);
    ret = (char *)malloc(ls.LengthOfLss+1);
    memcpy(ret, ls.LocaOfStr, ls.LengthOfLss);
    ret[ls.LengthOfLss] = NULL;
    va_end(ap);
    FreeRes(&s_FMS);
    return ret;
}

static void AllocResAndInit(pStrucForMultiSeq ps)
{
    ps->isFind = (bool *)malloc(ps->NumOfSeq*sizeof(bool));
    ps->record = (size_t *)malloc(ps->NumOfSeq*sizeof(size_t));
    ps->sign = (char **)malloc(ps->NumOfSeq*sizeof(char *));
    ps->suffix = (char **)malloc((ps->ls.LengthOfLss-ps->NumOfSeq)*sizeof(char *));
    if(!(ps->isFind && ps->record && ps->sign && ps->suffix)) {
        printf("Fail to Allocate Memory");
        system("pause"); exit(0);
    }
    size_t currlength = 0, LenOfStr;
    char *pStr = ps->ls.LocaOfStr;
    for(size_t index = 0; index != ps->NumOfSeq; index++)
    {
        LenOfStr = strlen(pStr);
        for(size_t strix = 0; strix != LenOfStr; strix++)
            ps->suffix[currlength+strix] = pStr++;
        currlength += LenOfStr;
        ps->sign[index] = pStr++;
    }
    memset(ps->isFind, 0, ps->NumOfSeq);
}

static void FreeRes(pStrucForMultiSeq ps)
{
    free(ps->isFind);
    free(ps->ls.LocaOfStr);
    free(ps->record);
    free(ps->sign);
    free(ps->suffix);
}

static logStr MergeMultiSeq(va_list *ap, size_t num)
{
    enum { INITSIZE = 1024, TIMES = 2 };  /*argument for dynamic array*/
    struct { size_t current, max; char *str; } da;
    char *pVar, *pda;
    size_t szLength = 0;
    logStr ls;
    da.current = 0, da.max = INITSIZE, da.str = (char *)malloc(da.max);
    for(size_t index = 0; index != num; index++)
    {
        pVar = va_arg(*ap, char*);
        szLength = strlen(pVar);
        while(szLength+1>da.max-da.current) {
            pda = (char *)realloc(da.str, da.max*TIMES);
            if(pda == NULL) {
                printf("Fail to Allocate Memory!");
                system("pause"); exit(0);
            }
            da.str = pda, da.max *= TIMES;
        }
        strcpy(&da.str[da.current], pVar);
        da.current += szLength+1;
        da.str[da.current-1] = NULL;
    }
    ls.LocaOfStr = da.str;
    ls.LengthOfLss = da.current;
    return ls;
}

static logStr FindLssAfterMergeAndSort(pStrucForMultiSeq ps)
{
    size_t ix, cnt = ix = 0, index;        /*ix is the index of current element in 'suffix'*/
    size_t len = ps->ls.LengthOfLss-ps->NumOfSeq;
    int k, temp;                /*the kth original sequence*/
    logStr ls = {0, 0};
    for( ; ; )
    {
        for(index = ix; index != len; index++)
        {
            k = -1;
            do k++;
            while(ps->suffix[index] > ps->sign[k]);
            if(ps->isFind[k] == false) {
                ps->isFind[k] = true; cnt++;
            }
            ps->record[k] = index;
            if(cnt == ps->NumOfSeq) break;
        }
        if(index == len) break;
        QuickSort(ps->record, ps->NumOfSeq, sizeof(ps->record[0]), i_comp);
        k = -1;
        do k++;
        while(ps->suffix[ps->record[0]] > ps->sign[k]);
        ps->isFind[k] = false;
        cnt--;
        temp = -1;
        do temp++;
        while(ps->suffix[ps->record[0]][temp] != NULL &&
            ps->suffix[ps->record[0]][temp] == ps->suffix[ps->record[ps->NumOfSeq-1]][temp]);
        if(temp > ls.LengthOfLss) {
            ls.LengthOfLss = temp; ls.LocaOfStr = ps->suffix[ps->record[0]];
        }
        ps->record[0] = ps->record[1];
        ix = ps->record[ps->NumOfSeq-1]+1;
    }
    return ls;
}

以上是在单个/多个序列中查找最长公共子串的算法，该算法提供了一种文本匹配的思想，具有借鉴意义。并且在有些情况下，排序是高效检索的前提，所以一个健壮而高效的通用排序算法也具有现实意义。

K阶马尔可夫随机文本
用于生成随机文本的一种比较简单的做法是根据所要生成的文本的长度，每次从26个英文字母及1个空格符中随机选出一个字符，最终构成一个序列。这样带来的问题是，文本的生成确实满足随机性的要求，然而似乎并不是给人看的，也就是说生成的文本不存在任何意义。那么如何生成有意义的随机文本呢？其实在回答这个问题之前，我们首先应该注意到的是现实世界中发生的任何事件所具有的意义均在于其上下文环境之中。比如甲问乙：你今年多大了？那么乙的回答应该是1~100之间的任意一个数，如果乙回答了1000，那么该答案显然不合上下文语境的要求。对于随机文本也是一样，某个字符的生成应当依赖于其上下文语境，这样至少是可以让人看懂的。最终的结果是，我们可以把随机文本的生成抽象成一个参数化的马尔可夫链，即
    在一个系统中，其在任意时刻均处于N个离散状态中的一个：S₁，S₂，…，S_N，时刻t的状态标记为Q_t，t=1，2，…
    设在某个离散时刻，根据之前的状态值转移到某个状态的概率为P(Q_t+1=S_j | Q_t=S_i，Q_t-1=S_k，…)
    对于一阶马尔可夫模型来说，系统在时刻t+1的状态仅仅依赖于其在时刻t的状态，因此转移概率为P(Q_t+1=S_j | Q_t=S_i，Q_t-1=S_k，…)=P(Q_t+1=S_j | Q_t=S_i)，下图是一阶马尔可夫链更形象的表示：

图5

其中a_ij=P(Q_t+1=S_j | Q_t=S_i)，即当前状态S_i以概率a_ij转移至状态S_j，这一概率在任何时刻t均相同且有a_ij≥0，∑a_ij=1(其中i=1，2，…)，π_i表示整个系统的首个状态为i的概率。

上述马尔可夫模型可以看成一个“随机自动机”（stochastic automation），选择一个初始状态并规定所要生成的序列长度，最终的结果则完全由系统自身决定，因此满足文本生成的随机性要求。并且由于系统的当前状态由前一时刻的状态值决定，因此我们可以保证该状态在局部上下文中存在意义，然而如果我们严格要求系统的某个状态由先前的所有状态决定，那么最终的生成结果将是完全预定义的，这种将”上下文”进行泛化所形成的模型我们称之为“K阶马尔可夫链”。有关该文本更详细的介绍参见信息论之父香农的著名论文[A Mathematical Theory of Communication]。

问题的关键在于：整个系统中的所有状态及其转移概率如何确定？这可以通过选择一个已有文本，根据选定文本的内在结构构建出每个状态的统计模型，举例来说：
设选定的文本为“an apple”，则有

原状态	转移状态
a	n、p
n	[blank]
[blank]	a
p	p、l
l	e
e	/

根据以上统计模型（其中[blank]表示空格），我们知道如果时刻t的状态为'a'，那么下一时刻的状态可以为'n'也可以为'p'，若所选择的状态为'p'，那么接下来可选择的的状态为'p'或'l'，依此类推，直至系统已经构造完所需长度的序列，或由某一状态无法再进行转移为止。一般来说，所选择的文本的内在结构越丰富，那么最终的随机文本也就更加有趣。另外要说明的一点是，关于统计模型的构造虽然是以字符为例，但实际生成英文随机文本时都将整个单词作为一个状态，然而如果要生成的是中文文本，那么只能以单个汉字作为一个状态来构造统计模型，注意这里不考虑对句子进行分词。好在无论以什么作为状态，最终生成随机文本的原理都是相同的。

由之前的分析可知一旦构造完整个统计模型，那么在生成随机文本的过程中将不再发生改变。并且在生成随机文本时将根据现有状态随机选择下一个状态，因此将产生大量的检索操作。所以可以考虑在读取完整个文本之后，建立后缀数组并对其进行排序，执行完这些预处理之后每次都通过二分搜索来定位，那么文本的生成效率更高。这里我们不考虑使用链接形式的散列结构，虽然它能提供高效的搜索操作，然而其致命的缺点是我们无法保证文本的任意两个状态之间都不会产生碰撞。在最后的统计模型中，每个原状态所对应的转移状态的个数均不相同，并且对文本的统计通过机器自动完成，我们不可能事先预知一个特定状态的所有转移状态，因此带来的问题是如何在未知个数的转移状态中提供无差别的随机选择？给出以下方法：

RandSelFrmUnkState
i←1
state←0
while CurrentState is still in the set of Translation States
  do if rand(0,i)=0  /*generate integer in the interval [0,i]*/
       then state←i
     i←i+1
return state

其基本思想是总是以1的概率选择第一个转移状态，以1/2的概率选择第二个，以1/3的概率选择第三个，依此类推。在这一过程结束之后，我们即能保证在未知个数的转移状态中选择任意一个状态的概率均相同。简单证明为，假设在n个状态中（n未知）选择第k个状态，因此选择该状态时的概率必为1/k，接下来必不选择第k+1个状态，此时的概率为k/(k+1)，类推直至最终的结果为P{state=k}=1/k×k/(k+1)×…×(n-1)/n=1/n，结论得证。

根据以上分析，简易的随机文本生成程序如下：

#include <stdio.h>
#include <string.h>
#include <time.h>
#include <stdarg.h>
#include <Windows.h>

/*binary search*/
extern int binasearch(const void *,const void *, size_t, size_t, comp);

#define BYTESOFCNCHAR 2
#define TEXTLENGTH 160
#define ORDER 3
static size_t k_order = 1;
static char *inputFile = ".\\origtext.txt";
static char *outputFile = ".\\markovtext.txt";
static char *getFileContext(const char *);
static void geneMarkovText();
static void getFirstState(char *,char *, size_t, size_t);
static int getNextState(char *, char **, size_t, size_t);
static int wordNComp(const void *, const void *);
static size_t RandSelFrmUnkState(char *, char **, size_t);

int main(int argc, char *argv[])
{
    srand((unsigned int)time(NULL));
    geneMarkovText();
    system("pause");
    return EXIT_SUCCESS;
}

static void geneMarkovText()
{
    char *textBuf = getFileContext(inputFile);
    size_t szBuf = strlen(textBuf)/2;
    char MarkovText[4096], *pMark = MarkovText;
    memset(MarkovText, 0, 4096);
    /*pay attention that a chinese character take over 2 bytes*/
    char **SuffixArray = (char **)malloc(sizeof(char *)*szBuf);
    FILE *fp = fopen(outputFile, "w");
    for(size_t index = 0; index != szBuf; index++)
        SuffixArray[index] = textBuf+index*BYTESOFCNCHAR;
    QuickSort(SuffixArray, szBuf, sizeof(SuffixArray[0]), s_comp);
    for(size_t k = 0; k != ORDER; k++)
    {
        sprintf(pMark, "\n%d-阶\n", k_order);
        pMark = MarkovText+strlen(MarkovText);
        getFirstState(textBuf, pMark, szBuf, k_order);
        for(size_t index = 0; index != TEXTLENGTH; index++, pMark += BYTESOFCNCHAR)
            if(!getNextState(pMark, SuffixArray, szBuf, k_order))
                break;
        k_order++;
    }
    printf("%s\n", MarkovText);
    fprintf(fp, MarkovText);
    free(textBuf);
    free(SuffixArray);
    fclose(fp);
}

static void getFirstState(char *src, char *dst, size_t len, size_t k)
{
    size_t Loca = GeneFrom(0, len-k);
    memcpy(dst, src+Loca*BYTESOFCNCHAR, k*BYTESOFCNCHAR);
}

static int getNextState(char *buf, char **SuffixArray, size_t len, size_t k)
{
    short *dst = (short *)buf, *src;
    size_t choice = 1;
    int location = binasearch(buf, SuffixArray, len, sizeof(char *), wordNComp);
    if(location == -1)
        return 0;
    while(location != len && !wordNComp(SuffixArray+location, buf))
    {
        if(GeneFrom(1, choice++) == 1)
            src = (short *)*(SuffixArray+location);
        location ++;
    }
    *(dst+k) = *(src+k);
    return 1;
}

static int wordNComp(const void *p, const void *q)
{
    char *p1 = *(char **)p, *p2 = (char *)q;
    size_t left = k_order*BYTESOFCNCHAR, index;
    for(index = 0; index != left; index++, p1++, p2++)
        if(*p1 != *p2) break;
    if(index == left) return 0;
    else return (unsigned char)*p1-(unsigned char)*p2;
}

static char *getFileContext(const char *fileName)
{
    FILE *fp = fopen(fileName, "r");
    size_t szFile, result;
    char *buf;
    if(fp == NULL) { printf("Wrong : file not exist!\n"); exit(1); }
    /*obtain file size*/
    fseek(fp, 0, SEEK_END);
    szFile = ftell(fp);
    rewind(fp);
    /*allocate memory and read file*/
    buf = (char *)malloc(sizeof(char)*szFile+1);
    if(buf == NULL) { printf("Fail to Allocate Memory!\n"); exit(1); }
    result = fread(buf, 1, szFile, fp);
    if(result != szFile) { printf("Reading error!\n"); exit(1); }
    buf[szFile] = NULL;
    fclose(fp);
    return buf;
}

extern int binasearch(
    const void *elem,
    const void *src,
    size_t len,
    size_t NumOfBytes, 
    comp f)
{
    int low = -1, high = len-1, mid;
    while(low+1 != high)
    {
        mid = (low+high)/2;
        if(f((char *)src+NumOfBytes*mid, elem) < 0)
            low = mid;
        else
            high = mid;
    }
    if(f((char *)src+NumOfBytes*high, elem) != 0)
        high = -1;
    return high;
}

下面我选了06年春晚的一段绕口令来测试上述程序，原文本如下：

且南来了个哑巴，腰里别着个喇叭；且北来了个喇嘛，手里提溜着五斤塌目；别喇叭的哑巴要拿喇叭去换提溜塌目喇嘛的塌目；提溜塌目喇嘛就不拿塌目换别喇叭哑巴的喇叭；别喇叭哑巴要拿喇叭打提溜塌目喇嘛一喇叭；提溜塌目喇嘛要拿塌目打别喇叭哑巴一塌目；不知道是提溜塌目喇嘛拿塌目打了别喇叭哑巴一塌目；还是别喇叭哑巴拿喇叭打了提溜塌目喇嘛一喇叭；哑巴吹喇叭，那个喇嘛回家，他就炖塌目。

以下是生成的随机文本

1-阶
提溜塌目；不拿喇嘛的喇嘛，那个喇叭哑巴一塌目喇嘛的哑巴一喇叭哑巴一塌目；且南来了个喇叭；提溜塌目；提溜塌目打了别着五斤塌目喇嘛一塌目；不拿塌目；提溜塌目；别喇嘛，手里提溜塌目；哑巴的喇嘛一塌目；还是提溜塌目打了提溜塌目打别着五斤塌目喇叭哑巴要拿喇嘛拿喇叭哑巴拿喇嘛就不拿喇嘛一塌目打别喇叭哑巴吹喇叭；提溜塌目喇嘛回家，手
2-阶
嘛就不拿塌目；别喇叭；提溜塌目；不知道是提溜塌目喇嘛回家，他就炖塌目喇嘛的塌目换别喇叭哑巴，腰里别着个喇嘛一喇叭的哑巴要拿喇叭哑巴的喇叭，那个喇嘛要拿塌目打别喇叭；哑巴一塌目；还是别喇叭，那个喇叭；提溜塌目；别喇叭，那个喇叭；提溜塌目打了别喇叭；哑巴要拿喇叭；别喇叭哑巴吹喇叭哑巴一塌目。
3-阶
溜塌目喇嘛就不拿塌目打了别喇叭哑巴的喇叭；且北来了个喇嘛，手里提溜着五斤塌目；提溜塌目喇嘛拿塌目换别喇叭哑巴拿喇叭打了提溜塌目喇嘛一喇叭；且北来了个哑巴，腰里别着个喇叭；哑巴吹喇叭，那个喇嘛，手里提溜着五斤塌目；提溜塌目喇嘛的塌目；不知道是提溜塌目喇嘛要拿塌目打了别喇叭哑巴的喇叭；提溜塌目喇嘛拿塌目换别喇叭哑巴一塌目；别喇叭

按理说，如果下一个状态所依赖的当前状态越短，则可供选择的状态越多，因而就更加缺乏内聚力，然而如果所依赖的状态越长，则趋向于原始输入的逐字拷贝。一般来讲，用连续的两个字符（即2-阶）选择下一个状态是一个较好的选择，因为此时既能重现原文本的味道，又伴随随机性带来的“恶搞”趣味。