POJ2778 DNA Sequence 题解（AC自动机+矩阵快速幂）

（题目描述略）

本题对于时间的要求比较严格，这意味着朴素的搜索算法是很难用优化手段通过的。为了解决这个问题，我们需要用一种不同于枚举的算法。

几乎任何一个问题都有其对应的图论模型，这个问题也不例外。我们可以将此问题转化为在其对应的图上求可接受的长度确定的路径条数。由于本题只有 AGCT 四个选择，所以其规模是我们可以接受的。

为了将这个问题转化为图论，我们首先要将其等价的图建立起来。根据定义，一个图由点集和边集组成。对于本题，我们若已知序列 DNA[1 .. n]，欲推导出 DNA[1 .. n+1]，我们可以将 DNA[1 .. n] 和 DNA[1 .. n+1] 作为点，即状态，将 DNA[n+1] 作为边，即状态的转移，记权值为 1 。当一个状态序列中包含禁阻子串时，我们将这个状态称为禁阻状态，即不合法的状态，而所有从其他状态转移到此状态（此点的入边）和从此状态转移到其他状态（此点的出边）都称为禁阻转移，即不合法的转移。我们将所有禁阻状态和禁阻转移从我们所构建的图中移去，则剩余图的所有状态必全为可接受的状态，即有效状态。我们再在此图中统计给定长度的路径条数，即为所求。

但是，以上做法不仅没有解决时间的问题，还增加了空间的问题。为了解决空间不足的问题，我们可以去除不需要的状态。设存在一个由有效状态 DNA[1 .. n] 转移到另一个有效状态 DNA[1 .. n+1]，且 n 大于最长的禁阻子串的长度，则在有效状态 DNA[1 .. n+1] 中 DNA[1] 为多余状态，因为其既不参与之后的转移，又不可能与之后的字符构成禁阻子串。推广开来，我们可以预处理出最长的禁阻子串的长度 length，则状态序列的最长长度即为 length - 1，当由有效状态 DNA[n-length+2 .. n] 转移到状态 DNA[n-length+3 .. n+1] 时，我们只需判断状态序列 DNA[n-length+3 .. n+1] 是否包含禁阻子串，若否则其必为有效状态。

尽管有上文的空间优化，其空间复杂度仍为指数级别。为了进一步优化空间，我们需要尽可能地压缩状态。由于禁阻子串的个数很小，我们可以发现在状态序列的某几位上某几个字符是几乎等价的。我们将这些状态合并，便可以使状态总数维持在可接受的范围之内。但是，这一步说起来容易，做起来却是比较复杂的。为了更高效地完成这一操作，我们将使用一种线性时间复杂度和线性空间复杂度的数据结构，即本题的正解 AC自动机 。

AC自动机是一种可以接受给定字符串集合内所有字符串的确定有穷状态自动机，可以在线性时间内构造。其具体实现方法在这里不做过多介绍。根据禁阻序列所构造的AC自动机其实已经完成了状态的压缩，我们只需删去所有禁阻状态和禁阻转移便可实现空间的优化。关于以上结论的解释如下，因为AC自动机中，我们每插入一个节点，便对其父节点进行了一次状态的划分，即将父节点所对应的剩余有效状态划分成通向新禁阻状态的转移和通向新有效状态的转移，而在当前状态下所有通向有效状态的转移是等价的。即使一个节点的状态转移被划分多次，剩余有效状态转移仍然等价。例如对于禁阻序列 {“AT”,”AC”,”AG”,”AA”}，当插入 “AT” 和 “AC” 完成后，状态 “A” 被划分成 ‘T’ 转移和 ‘C’ 转移，而此时状态 “A” 的 ‘G’ 转移和 ‘A’ 转移是等价的。基于以上原理，一个和问题等价的图论模型就建立起来了。

接下来我们要解决的是时间的问题。图建立好之后，我们要使用一种高效的手段计算以源点为路径起点的长度为 n 的不同路径的总数。在这里，我们仍然用状态转移的思想进行分析。我们记路径长度为 i，起始点 j，终止点 k 的不同路径总数为 f[i, j, k]， 则 f[i, j, k] = sum(f[i-1, j, l] * f[1, l, k]) （1 ≤ l ≤ m，m 为图中点的总数）（很像 Floyd）。翻译成中文，意思即为将长度为 i 的路径划分为长度为 i - 1 的路径和长度为 1 的路径， 枚举划分点求和。当然，你也可以将其划分为长度 i - 2 的路径和长度 2 的路径，或长度为 i - 3 的路径和长度 3 的路径，等等。 然而为了更好地探究其性质，划分成长度 i - 1 和长度 1 的路径方便一些。在列出上文的动态转移方程之后，我们惊奇地发现这与一个数据结构的运算非常相似，这个数据结构就是 矩阵，而这种运算就是 矩阵乘法。

关于矩阵及矩阵乘法，我在这里也不再做过多的介绍，只说明一点，那就是矩阵乘法的运算法则：

$$A_{p*q} * B_{q*r} = C_{p*r}\\C(i,j) = \sum_{k=1}^q A(i,k)*B(k,j)$$ 同时，矩阵乘法不满足交换律，满足结合律，这使得矩阵乘法可以用 矩阵快速幂 进行运算。我们将图用邻接矩阵的形式储存，记为 A，即 A[i, j] 表示从 i 到 j 的路径条数，初始路径长度为 1，将 A 看做矩阵， $\sum_{i=1}^m A^n(1,i)$即为所求。这样，就使得时间复杂度也降低到了我们可以接受的范围之内。

最后还有几个小问题。其一是结构体赋值的问题。对于矩阵乘法的运算，由于在矩阵快速幂中不止出现一次，最好写子程序或重载乘法运算符。这样不可避免地出现用结构体定义矩阵的情况，而结构体赋值便成了一个纠结的问题。其实，在 C/CPP 中，结构体赋值是允许的，其原理应该是逐位赋值。于是，重载赋值运算符便成了不必要的了。

其二是整型溢出的问题。虽然题目描述中提到最后结果对 100,000 取模，但事实上矩阵乘法运算过程中仍有整型溢出的可能。为解决这个问题，long long 超长整型需要在运算过程中使用。

其三是矩阵 1 的问题。对于矩阵乘法，我们需要定义一个初始矩阵 C，使得任意矩阵 A * C = A，即类似于数 1 的效果。这样的矩阵是存在的，其为:

$$C(i,j) = \begin{cases}1 & \text{i = j}\\0 & \text{i ≠ j}\end{cases}$$

代码如下：

#include"stdio.h"
#include"string.h"
bool end[105];
char s[15];
int f[105],fail[105],m,next[105][4],q[105],top=1;
struct Matrix
{
    long long data[105][105];
}mat;
int reflect(char c)
{
    switch(c)
    {
        case 'A':return 0;
        case 'C':return 1;
        case 'G':return 2;
        case 'T':return 3;
    }
}
void Insert(char str[])
{
    int now=0;
    for(int i=0;str[i];now=next[now][reflect(str[i++])])
        if(next[now][reflect(str[i])]==-1)
            next[now][reflect(str[i])]=top++;
    end[now]=false;
}
void Build()
{
    top=1;
    for(int i=0;i<4;i++)
        if(next[0][i]>-1)
            fail[next[0][i]]=0,q[top++]=next[0][i];
        else
            next[0][i]=0;
    for(int i=1;i<top;i++)
        for(int j=0;j<4;j++)
            if(next[q[i]][j]>-1)
                fail[next[q[i]][j]]=next[fail[q[i]]][j],q[top++]=next[q[i]][j];
            else
                next[q[i]][j]=next[fail[q[i]]][j];
}
Matrix Multiply(Matrix A,Matrix B)
{
    Matrix C;
    for(int i=0;i<m;i++)
        for(int j=0;j<m;j++)
        {
            C.data[i][j]=0;
            for(int k=0;k<m;k++)
                C.data[i][j]=(A.data[i][k]*B.data[k][j]%100000+C.data[i][j])%100000;
        }
    return C;
}
Matrix pow(Matrix A,int n)
{
    Matrix C;
    for(int i=0;i<m;i++)
        for(int j=0;j<m;j++)
            C.data[i][j]=(i==j)?1:0;
    for(Matrix B=A;n>0;n>>=1,B=Multiply(B,B))
        if(n&1)
            C=Multiply(B,C);
    return C;
}
int add(int n)
{
    int sum=0;
    mat=pow(mat,n);
    for(int i=0;i<m;i++)
        sum=(mat.data[0][i]+sum)%100000;
    return sum;
}
int main()
{
    int n;
    scanf("%d %d\n",&m,&n);
    memset(end,true,sizeof(end));
    memset(next,-1,sizeof(next));
    while(m--)
        gets(s),Insert(s);
    Build();
    m=0;
    memset(f,-1,sizeof(f));
    for(int i=0;i<top;i++)
    {
        for(int temp=i;temp>0&&end[i];temp=fail[temp])
            end[i]&=end[temp];
        if(end[i])
            f[i]=m++;
    }
    memset(mat.data,0,sizeof(mat.data));
    for(int i=0;i<top;i++)
        if(end[i])
            for(int j=0;j<4;j++)
                if(f[next[i][j]]>-1)
                    mat.data[f[i]][f[next[i][j]]]++;
    printf("%d",add(n));
    return 0;
}

对于样例 {“AT”,”AC”,”AG”,”AA”}，所建矩阵

$$\begin{pmatrix}3 & 1\\0 & 0\\\end{pmatrix}$$ 矩阵快速幂后结果 $$\begin{pmatrix}27 & 9\\0 & 0\\\end{pmatrix}$$

本题非常经典，建议编写。由于本题样例过弱，不建议用样例检验程序的正确性。