新冠病毒的检测需要采样咽拭子，然后将其中含有的碱基对解螺旋、经过聚合酶链式反应（PCR）扩增后与相应的荧光探针结合，然后进行检测。

而制造能与病毒碱基序列结合的荧光探针，需要知道病毒的碱基序列。

在这个过程中，有一个算法发挥了巨大的作用。

如今的这个算法已经变得十分基础，也有更多的算法超越了它。

但是，它诞生的那一刻，必定闪耀着人类智慧的光芒。

我第一次见到这个算法时，经历了拒绝、疑惑、不解。

直到后来，一切全都化为了惊叹。

这个算法叫KMP算法（Knuth、Morris和Pratt对该算法进行了改进，称为KMP算法）。由于三个字母的原因，有人戏称为“看毛片”算法。但它和“毛片”真的没有一点关系。KMP实际是合力研究出该算法的三位大神的姓名首字母。是的，这个算法是当时三个大神共同努力的结果。

KMP算法实际是一个字符串匹配算法。

子串的定位运算称为串的模式匹配或串匹配。

假设有两个串S，T，设S为主串，也称正文串，T为子串，也称为模式，在主串S中查找与模式T相匹配的子串，如果查找成功，返回匹配的子串第一个字符在主串中的位置。

那它如何为基因测序做出贡献，又为什么让我拒绝、疑惑、不解，最后惊叹。我们一起来重新感受下。

首先，我们先给出算法解决的问题，这是一个极为常见的字符串匹配问题。

存在一个字符串S，请找出其是否含有子串T，如果有，请给出T所在的位置。

假设字符串S的长度为m，字符串T的长度为n。

则见到这个问题后，有经验的程序员不出5秒就能知道其时间复杂度为O(m×n)。是的，这是很多人的第一反应，使用暴力枚举法。

字符串S我们常称为主串，而T则被称为模式串。在主串的任意一个位置，我们都要以此为起点和模式串依次对比。模式串长度为n，而主串中有m个位置，则复杂度为O（m×n）。

#include<iostream>#include<cstring>using namespace std;#define Maxsize 100typedef char SString[Maxsize+1];//0号单元存放串的长度bool StrAssign(SString &T,char *chars)//生成一个其值等于chars的串T{	int i;	if(strlen(chars)>Maxsize)		return false;	else    {		T[0]=strlen(chars);		for(i=1;i<=T[0];i++)        {            T[i]=*(chars+i-1);            cout<<T[i]<<"  ";        }        cout<<endl;		return true;	}}int Index_BF(SString S,SString T,int pos) // BF算法{ 	// 求T在主串S中第pos个字符之后第一次出现的位置	// 其中，T非空，1≤pos≤s[0],s[0]存放S串的长度	int i=pos,j=1,sum=0;	while(i<=S[0]&&j<=T[0])    {        sum++;        if(S[i]==T[j]) // 如果相等，则继续比较后面的字符		{			i++;			j++;		}		else		{			i=i-j+2; // i退回到上一轮开始比较的下一个字符			j=1;     // j退回到第1个字符		}    }	cout<<"一共比较了"<<sum<<"次"<<endl;	if(j>T[0])      // 匹配成功		return i-T[0];	else		return 0;}int main(){	SString S,T;	char str[100];    cout<<"串S："<<"  ";    cin>>str;        // abaabaabeca	StrAssign(S,str);// 可以修改程序，自己输入字符串     cout<<"串T："<<"  ";    cin>>str;        // abaabe	StrAssign(T,str);	cout<<endl;	cout<<"BF算法运行结果："<<endl;	cout<<"主串和子串在第"<<Index_BF(S,T,1)<<"个字符处首次匹配\n";	cout<<endl;	cin.get();cin.get();	return 0;}/* 输入abaabaabecaabaabe*//* 运行效果串S：  abaabaabecaa  b  a  a  b  a  a  b  e  c  a串T：  abaabea  b  a  a  b  eBF算法运行结果：一共比较了15次主串和子串在第4个字符处首次匹配*/

听说有个算法能做到复杂度O(m+n)的时候，听说那个算法叫KMP。

怎么可能呢，O(m+n)意味着比较指针不能回撤，直觉觉得是不可能的。

搜了一下各种资料。

竟然真的是……

如以下字符串匹配：

按暴力枚举法，当元素匹配不成功时，i、j分别要回退：

i=i-j+2; // i退回到上一轮开始比较的下一个字符j=1;     // j退回到第1个字符

是否需要如此回退呢？

看此时能匹配的S、T的公共子串G：abaab。（可以设S的G部分为GS，T的G部分为GT）

因为是依次后移的匹配，可以考虑在GS中的最大后缀（从后向前取若干字符）去匹配GT的前缀（从前向后取若干字符），如以上数据，能匹配的最大子串是ab，ab的长度为L。可以证明，i、j的回退，i不动，j回退到L+1再去比较，算法的正确性不受影响。

而在GS中的最大后缀去匹配GT的前缀，因为GS=GT，也就是在GT中的最大后缀去匹配GT的前缀。也是一个在自身字符串中子串的匹配问题，随着j的不同，如何去求L呢？可以是暴力枚举：

可以用next[j]表示j要回退的位置。

当j=1时，next[j]=0;

当T'的相等前缀和后缀的最大长度为Lmax时，next[j]=Lmax+1；

当没有相等的前缀和后缀时，next[j]=1。

根据以上规则，可以很容易求出T="abaabe“的next数组：

1 2 3 4 5 6 | j

a b a a b e | T

0 1 1 2 2 3 | next[j]

对于回退位置next[]的求解方法，KMP算法应用的却是动态规划的思想：

首先大胆假设next[j]=k，则

那么next[j+1]=?

也就是以上子串再分别多考虑一个随后的字符：

可以区分这两个字符在相等和不等的情况下分别考虑：

回退位置next[]的求解方法的代码实现：

void get_next(SString T,int next[]) //计算next函数值{	int j=1,k=0;	next[1]=0;	while(j<T[0])      // 动态规划	{	    if(k==0 || T[j]==T[k])            next[++j]=++k;		else			k=next[k]; // 下标k的回退	}    cout<<"-----next[]-------"<<endl;    for(j=1;j<=T[0];j++)        cout<<next[j]<<"  ";    cout<<endl;}

KMP算法的代码实现，与暴力枚举相比，仅仅在于i、j的回退位置不同（i不回退）以及求next[j]的动态规划思想：

int Index_KMP(SString S,SString T,int pos,int next[]) // KMP算法{ 	// 利用模式串T的next函数求T在主串S中第pos个字符之后的位置的KMP算法	// 其中，T非空，1≤pos≤StrLength(S)	int i=pos,j=1,sum=0;	while(i<=S[0]&&j<=T[0])    {        sum++;        if(j==0||S[i]==T[j]) // 继续比较后面的字符		{			i++;			j++;		}		else			j=next[j];       // 模式串向右移动    }	cout<<"一共比较了"<<sum<<"次"<<endl;	if(j>T[0])               // 匹配成功		return i-T[0];	else		return 0;}

下面给出一个KMP算法的实现动画，大家关注一下最上面S串的比较指针，S串真的没有任何回撤操作。

KMP算法动画演示：

视频加载中...

等你真正理解了它，发现它并没有太复杂，主要是利用了模式串中的相关性。甚至，会有更优的算法被提出。

void get_next2(SString T,int next[]) // 计算next函数值改进算法{	int j=1,k=0;	next[1]=0;	while(j<T[0])	{	    if(k==0||T[j]==T[k])        {            j++;            k++;            if(T[j]==T[k])                next[j]=next[k];            else                next[j]=k;        }		else			k=next[k];	}    cout<<"-----next[]-------"<<endl;    for(j=1;j<=T[0];j++)        cout<<next[j]<<"  ";    cout<<endl;}

但是KMP算法由D.E.Knuth、J.H.Morris、V.R.Pratt发表于1977。在那个年代，相信KMP算法诞生时，一定闪耀着光芒。

那KMP算法为什么会和基因扯上关系呢？

疫情爆发，专家发现了一种新型环状病毒，这种病毒的DNA序列是环状的，而人类的DNA序列是线性的。专家把人类和病毒的DNA表示为字母组成的字符串序列，如果在某个患者的DNA中发现这种环状病毒，说明该患者已被感染病毒，否则没有感染。

假设病毒的DNA为aabb，患者的DNA为eabbacab，说明该患者已被感染。

如何处理环状病毒呢？

1.环形处理

使用循环存储的方式，循环队列或循环链表。

2.线性处理

将病毒序列扩大两倍，依次从每个下标开始，取m个字符。

算法步骤:

（1） 首先对环状病毒进行处理（环形处理或线性处理）。

（2） 依次把每一个环状病毒变种作为子串，把患者DNA序列作为主串，进行模式匹配。若匹配成功，返回已感染病毒；

（3） 重复运行（2）步。

（4） 如果检测所有病毒变种都未匹配成功，返回未感染病毒。

例如：患者的DNA 序列eabbacab，病毒DNA序列：aabb，检测患者是否感染病毒。

代码实现：

/***病毒检测算法***/#include<cstring>#include<iostream>using namespace std;#define Maxsize 100typedef char SString[Maxsize+1];		//0号单元存放串的长度bool StrAssign(SString &T, char *chars)//生成一个其值等于chars的串T{	int i;	if(strlen(chars)>Maxsize)		return false;	else    {		T[0]=strlen(chars);//0号单元存放串的长度		for(i=1;i<=T[0];i++)        {            T[i]=*(chars+i-1);            cout<<T[i]<<"  ";        }        cout<<endl;		return true;	}}void get_next(SString T,int next[])//计算next函数值{	int j=1,k=0;	next[1]=0;	while(j<T[0])	{	    if(k==0||T[j]==T[k])            next[++j]=++k;		else			k=next[k];	}}int Index_KMP(SString S, SString T, int pos)//KMP算法{ 	// 利用模式串T的next函数求T在主串S中第pos个字符之后的位置的KMP算法	//其中，T非空，1≤pos≤StrLength(S)	int *next=new int[T[0]+1]; // 生成T的next数组	get_next(T,next);	int i=pos, j=1;	while(i<=S[0]&&j<=T[0])    {        if(j==0||S[i]==T[j]) // 继续比较后面的字符		{			i++;			j++;		}		else			j=next[j]; // 模式串向右移动    }	if(j>T[0]) // 匹配成功		return i-T[0];	else		return 0;}bool Virus_detection(SString S, SString T) // 病毒检测{    int i,j;    SString temp;                    // temp记录病毒变种    for(i=T[0]+1,j=1;j<=T[0];i++,j++)// 将T串扩大一倍，T[0]为病毒长度        T[i]=T[j];    for(i=0;i<T[0];i++)        // 依次检测T[0]个病毒变种    {        temp[0]=T[0];          // 病毒变种长度为T[0]        for(j=1;j<=T[0];j++)   // 取出一个病毒变种            temp[j]=T[i+j];        if(Index_KMP(S,temp,1))// 检测到病毒            return 1;    }    return 0;}int main(){	SString S,T;	char str[100];    cout<<"患者DNA序列S："<<"  ";    cin>>str; // eabbacab	StrAssign(S,str);    cout<<"病毒DNA序列T："<<"  ";    cin>>str; //aabb	setbuf(stdin,NULL);	StrAssign(T,str);	if(Virus_detection(S,T))        cout<<"该患者已感染病毒！"<<endl;    else        cout<<"该患者未感染病毒！"<<endl;	cin.get();	return 0;}/*输入eabbacabaabb*//*运行效果患者DNA序列S：  eabbacabe  a  b  b  a  c  a  b病毒DNA序列T：  aabba  a  b  b该患者已感染病毒！*/

平时，我们的文本，也就几百上万的字符。即使是做大数据分析，接触到的字符串也不过几十万、上百万个。当然，在这个尺度下，O（m×n）算法和O（m+n）的性能差别已经是极为巨大的。

然而，人类的基因序列，就是腺嘌呤（A）、鸟嘌呤（G）、胞嘧啶（C）、胸腺嘧啶（T）组成的那些ATCGCTACG……等朴实无华的序列串，有30亿个（即30亿对碱基对）。

基因片段的拼接需要不断在上亿的序列（碱基对）中找一个目标序列，O（m×n）算法和O（m+n）的性能差别有万倍、百万倍，甚至更多。

只用暴力匹配的方式，难以完成。

而KMP算法将O（m×n）降到了O（m+n）。只需要读取一遍目标串和模式串，就能给出结果！

在KMP算法诞生的1977年，那时人类才发现DNA的双螺旋结构不过20年，再过大约15年后人类基因组计划才会开始。

当然，在KMP算法诞生后人类基因组计划开始前的15年中，字符串匹配算法可能又会有进步，并且字符串匹配算法也并不最适合基因测序。是新的算法（例如OLC、de Bruijn Graph等算法）用在了基因组计划中。但KMP算法在匹配方面的研究已经打好了基础，它已经提前为人类基因组计划扫除了部分算法上的障碍。

之后，也正因为基因测序技术的发展，人类能迅速锁定出新冠病毒的序列，然后研发测试试剂。帮我们开展测试，推动人类抵御这场突如其来的病毒。

所以，在这场全世界抗击疫情的斗争中，有算法的身影。低调、朴实，但也无可替代。

如今的它，低调，朴实，甚至让人觉着习以为常。

但我知道，它曾经闪耀着人类智慧的光芒；而今，也依旧令人惊艳和伟大。

－End－