poj2778 DNA Sequence （AC自动机+矩阵乘法）

poj2778 DNA Sequence

题目大意：一段DNA序列里面如果带有某些特征序列，那么这段DNA序列就有问题。问给出若干段DNA序列，再给定一个长度，问这个长度的所有可能序列中，有多少是没问题的。
嗯，总之先建机行事。
先建好Trie树，升级成AC自动机。由于m很小，每个有病序列的长度也很小，但是n很大，所以想到用快速幂优化，邻接矩阵表示各个状态之间的关系。
假设矩阵$A_{n×n}$中的第i行，第j列元素$A_{ij}$表示从i状态走1步到达j状态有几种方法，$B_{n×n}$矩阵中的第j行，第k列元素$\ B_{kj}$表示从k状态走1步到达j状态有几种方法，那么矩阵$\ C_{n×n} = A × B$中的$\begin{equation*}\ C_{ij} = \sum_{k=1}^n(A_{ik} × B_{kj}) \end{equation*}$ 即表示从i走到j走两步有几种方法。

通过建好的自动机，可以找到每个状态走一步可以到达的下一个状态，由状态之间的关系可以建立一个邻接矩阵，其中第i行第j列表示从状态i走一步到状态j有几种方法，然后这个矩阵的n次方就是每个状态在经过n步之后的结果，其中第一行的总和就是从树根到达其它状态的n步之后一共有几种方式。在跑自动机的时候，需要筛掉那些到达有病节点的状态。

随手写个样例：
3 4
ACC
AGT
GT
邻接矩阵应该长这样（可能行列位置会有变化）：
2 1 0 0 0 0 1 0
1 1 1 0 1 0 0 0
1 1 0 0 0 0 1 0
0 0 0 0 0 0 0 0
1 1 0 0 0 0 1 0
0 0 0 0 0 0 0 0
1 1 0 0 0 0 1 0
0 0 0 0 0 0 0 0
结果应该是201

另外（为了偷懒？？），用了一点trivial的trick，就是ACTG的ASCII码的二进制表示的倒数第第二个和第三个二进制数刚好是从0~3的（A00 C01 T10 G11），所以可以直接将字符右移一位然后&3就可以完成从A到G的字符分选了。
代码比较冗余，也没有怎么优化过，各种手动开空间，一位数组当二维数组做，乘法找地址，跑起来速度一般（172ms）Orz

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<queue>
#include<cstdlib>
using namespace std;
struct node
{
    node *fail;
    node *son[4];
    int index;
    bool vis1, vis2, isdest;
};
char s[20];
node *root;
int index = 0;
long long *matrix, Size = 1, bytesize = 0, qwordsize = 0;
void multiply(long long *a, long long *b, long long *save)
{
    memset(save, 0, bytesize);
    for(int i = 0; i < Size; i++)
        for (int j = 0; j < Size; j++)
        {
            long long *aa = a + i*Size, *bb = b + j, *sv = save + i*Size + j;
            for (int k = 0; k < Size; k++)
                *sv += (*(aa + k) )*( *(bb + k*Size));
            *sv %= 100000;
        }
    memcpy(a, save, bytesize);
}
void f_pow(long long a[], int b)
{
    long long ans = 0;
    long long *base = new long long[qwordsize];
    long long *save = new long long[qwordsize];
    memset(base, 0, bytesize);
    for (int i = 0; i < Size; i++)
        base[i*Size + i] = 1;
    while (b > 0)
    {
        if (b & 1)
            multiply(base, a, save);
        multiply(a, a, save);
        b >>= 1;
    }
    for (int i = 0; i < Size; i++)
        ans += base[i];
    printf("%d\n", ans % 100000);
}
void buildtrie(char *s, node *root)
{
    int len = strlen(s), i;
    node *temp = root, **son;
    for (i = 0; i < len - 1; i++)
    {

        int ch = (s[i] >> 1) & 0x3;
        son = &temp->son[ch];
        if (*son == NULL)
        {
            *son = new node;
            memset(*son, 0, sizeof(node));
            Size++;
            index++;
            (*son)->index = index;
        }
        temp = *son;
    }
    int ch = (s[i] >> 1) & 0x3;
    son = &temp->son[ch];
    if (*son == NULL)
    {
        *son = new node;
        memset(*son, 0, sizeof(node));
        Size++;
        index++;
        (*son)->index = index;
    }
    (*son)->isdest = true;
}
void buildfail(node *root)
{
    queue<node*> que;
    for (int i = 0; i < 4; i++)
    {
        node *temp = root->son[i];
        if (temp != NULL)
        {
            que.push(temp);
            temp->fail = root;
        }
    }
    while (!que.empty())
    {
        node *temp = que.front();
        node **son = temp->son;
        que.pop();
        for (int i = 0; i < 4; i++)
        {
            if (son[i] != NULL)
            {
                que.push(son[i]);
                node *t = temp;
                do
                {
                    t = t->fail;
                    if (t->son[i] != NULL)
                        break;
                } while (t != root);
                if (t->son[i] != NULL)
                    son[i]->fail = t->son[i];
                else
                    son[i]->fail = root;
            }
        }
    }
}
bool finddest(node *temp)
{
    if (temp->vis1)
        return temp->isdest;
    if (temp->isdest)
        return true;
    if (temp != root)
        temp->isdest = finddest(temp->fail);
    temp->vis1 = true;
    return temp->isdest;
}
void ac_automation(node *root)
{
    node *temp = root;
    queue<node*> que;
    que.push(root);
    root->vis2 = true;
    while (!que.empty())
    {
        node *temp = que.front();
        for (int i = 0; i < 4; i++)
        {
            node *t = temp;
            while (t != root)
            {
                if (t->son[i] != NULL)
                    break;
                t = t->fail;
            }
            if (t->son[i] != NULL)
            {
                if (!finddest(t->son[i]))
                {
                    if (!t->son[i]->vis2)
                    {
                        que.push(t->son[i]);
                        t->son[i]->vis2 = true;
                    }
                    matrix[temp->index*Size+t->son[i]->index]++;
                }
            }
            else
            {
                if(root->son[i]!=NULL)
                    matrix[temp->index*Size+root->son[i]->index]++;
                else
                    matrix[temp->index*Size]++;
            }
        }
        que.pop();
    }
}
int main()
{
    int m, n;
    scanf("%d%d", &m, &n);
    root = new node;
    memset(root, 0, sizeof(node));
    while (m--)
    {
        scanf("%s", &s);
        buildtrie(s, root);
    }
    matrix = new long long[Size*Size];
    qwordsize = Size*Size;
    bytesize = qwordsize << 3;
    memset(matrix, 0, bytesize);
    buildfail(root);
    ac_automation(root);
    f_pow(matrix, n);
}