计蒜客 疑似病毒 (AC自动机 + 可达矩阵)

---

**链接 : ** Here!

背景 : ** 开始我同学是用 AC自动机 + DP 的方法来做这道题, 这道题的标签是 AC自动机, 动态规划, 矩阵, 按道理来说 AC自动机 + DP 应该是能过的, 但是他不幸的 $T$ 了, $QAQ$, 后来神犇Hug**给我们提供了一个思路!!!

思路 : ** 题目要求是 "如果其中两个或者两个以上**的 $DNA$ 序列是一个 $DNA$ 序列 $A$ 的子串，那么 $DNA$ 序列 $A$ 是疑似病毒的 $DNA$ 序列", 那么也就是说字符串 $A$, 它要是能在自动机中走到两次终态, 那么它就是疑似病毒的 $DNA$ 序列. 因为AC自动机可以等价转化为一个图, 所以说这个问题就变成了, 从 $root$ 节点出发, 在一个图中 $l$ 步通过终态点两次或者两次以上的方法数....................!!! $GG$

神犇告诉我们, 可以一个自动机在逻辑层面抽象成三层自动机, 然后用搜索来初始化可达矩阵, 这个可达矩阵$Matrix[i][j] = cnt$ 代表从 $i$ 点经过 $1$ 步到达 $j$ 点的方法数为 $cnt$. 然后将初始的可达矩阵求 $l$ 次方, 就可以得到一个增量矩阵 $Matrix'$ , 最后只需要用初始的矩阵 $[1, 0, 0, ....] * Matrix'$ 就可以得到答案了.

**细节 : **

1. 在初始化可达矩阵的时候, 对于任意一个节点, 它都能接受 $A, T, C, G$ 四个字符!!!, (因为这个地方wa了无数次).
2. 搜索时从当前层究竟能调到第几层的条件是自动机的节点的cal_value的所决定的.
3. 注意建立自动机的时候需要注意, 靠下的终态点需要加上它所有 $fail$ 指针的base_value.

补充 : ** 这道题的确非常有价值, 能够加深对AC自动机的理解, 而且还有将字符串匹配问题转为图的问题, 进而用可达矩阵**来表示连通性, 然后通过矩阵快速幂来求出 $l$ 长度的 $DNA$ 序列种类数.

**最后 : ** 待补充吧... 看起来理解还是不太够, 虽然听着思路能A, 但是找到方向才是更重要的事情!!!

---

/*************************************************************************
	> File Name: ac_machine.cpp
	> Author: 
	> Mail: 
	> Created Time: 2017年11月25日 星期六 11时01分11秒
 ************************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;

const int MAX_N = 10000;
const int MAT_X = 200;
const int SIZE = 4;
const int BASE = '0'; 
const int LAYER = 3;
const int MOD = 10007;
int node_cnt = 1;  // 统计节点数量用于分配index, 在insert中统计, 在automaton中分配
int mat_size;

// Matrix 结构声明
struct Matrix {
    Matrix() {
        memset(this->m, 0, sizeof(this->m));
    }
    int m[MAT_X][MAT_X];
    void show() {
        for (int i = 0 ; i < mat_size ; ++i) {
            for (int j = 0 ; j < mat_size ; ++j) {
                printf("%4d ", this->m[i][j]);
            }
            printf("
");
        }
    }
};
Matrix unit_mat;
Matrix CGMatrix; // 可达矩阵

void init_unint_mat() {
    for (int i = 0 ; i < MAT_X ; ++i) {
        unit_mat.m[i][i] = 1;
    }
}

Matrix operator* (const Matrix &a, const Matrix &b) {
    Matrix c;
    for (int i = 0 ; i < mat_size ; ++i) {
        for (int j = 0 ; j < mat_size ; ++j) {
            int sum = 0;
            for (int k = 0 ; k < mat_size ; ++k) {
                sum += (a.m[i][k] * b.m[k][j]);
                sum %= MOD;
            }
            c.m[i][j] = sum % MOD;
        }
    }
    return c;
}

// 快速幂 : 计算矩阵a的x次方
Matrix quick_mat_pow(Matrix a, int x) {
    Matrix ret = unit_mat;
    while (x) {
        // TODO
        if (x & 1) ret = ret * a;
        a = a * a;
        x >>= 1;
    }
    return ret;
}

// Trie 结构声明
// base_val : 基础值
// cal_val  : 计算值
// index    : 下标
typedef struct Trie {
    int base_val;
    int cal_val[LAYER];
    int index[LAYER];
    struct Trie *fail;
    struct Trie **childs;
} Node, *Tree;

Trie* new_node() {
    Trie *p = (Trie *)malloc(sizeof(Trie));
    p->childs = (Trie **)malloc(sizeof(Trie *) * SIZE);
    memset(p->childs, 0, sizeof(Trie *) * SIZE);
    p->fail = NULL;
    p->base_val = 0;
    memset(p->index, 0, sizeof(int) * LAYER);
    memset(p->cal_val, 0, sizeof(int) * LAYER);
    return p;
}

void clear(Trie *node) {
    if (node == NULL) return;
    for (int i = 0 ; i < SIZE ; ++i) {
        clear(node->childs[i]);
    }
    free(node->childs);
    free(node);
}

void insert(Trie *node, char *str) {
    Trie *p = node;
    for (int i = 0 ; str[i] ; ++i) {
        if (p->childs[str[i] - BASE] == NULL) {
            p->childs[str[i] - BASE] = new_node();
            // 更新节点数量
            ++node_cnt;
        }
        p = p->childs[str[i] - BASE];
    }
    p->base_val = 1;
}

// 建立自动机
void build_automaton(Trie *root) {
    root->fail = NULL;
    Trie *que[MAX_N];
    int l = 0, r = 0, k = 0; // k用于计算index
    que[r++] = root;
    while (l < r) {
        Trie *now = que[l++];
        
        // 更新三层下标
        now->index[0] = k;
        now->index[1] = k + node_cnt;
        now->index[2] = k + node_cnt * 2;
        ++k;
        
        for (int i = 0 ; i < SIZE ; ++i) {
            if (now->childs[i] == NULL) continue;
            Trie *child = now->fail;
            while (child && child->childs[i] == NULL) {
                child = child->fail;
            }
            if (child == NULL) {
                child = root;
            } else {
                child = child->childs[i];
            }
            now->childs[i]->fail = child;
            now->childs[i]->base_val += now->childs[i]->fail->base_val;
            que[r++] = now->childs[i];
        }
    }
}

// 判断在第几个自动机
int inLayer(int x) {
    return (x <= 2 ? x : 2);
}

// 得到孩子的下标
int getChildIndex(Trie *now, Trie *child, int now_ind) {
    return child->index[inLayer(now->cal_val[now_ind / node_cnt] + child->base_val)];
}

// 更新计算权值
int updataCalVal(Trie *now, Trie *child, int now_ind) {
    return now->cal_val[now_ind / node_cnt] + child->base_val;
}

// BFS初始化可达矩阵
void BFS(Trie *root) {
    Trie *que[MAX_N * 3];
    int ind[MAX_N * 3];
    int vis[MAX_N * 3] = {0};
    int que_l = 0, que_r = 0;
    int ind_l = 0, ind_r = 0;
    que[que_r++] = root;
    ind[ind_r++] = 0; 
    while (ind_l < ind_r) {
        Trie *now = que[que_l++];
        int now_ind = ind[ind_l++];
        vis[now_ind] = 1;
        for (int i = 0 ; i < SIZE ; ++i) { 
            Trie *child;
            if (!now->childs[i]) {
                // 寻找失败指针中是否出现childs[i], 如果没出现过, 那么就会走回root节点
                Trie *p_fail = now->fail;
                while (p_fail != NULL && p_fail->childs[i] == NULL) {
                    p_fail = p_fail->fail;
                }
                // 如果p_fail == NULL 那么这个一定为root
                if (p_fail == NULL) {
                    child = root;
                } else {
                    child = p_fail->childs[i];
                }
            } else {
                child = now->childs[i];
            }
            int child_ind = getChildIndex(now, child, now_ind);
            child->cal_val[child_ind / node_cnt] = updataCalVal(now, child, now_ind);
            CGMatrix.m[now_ind][child_ind]++;
            if (vis[child_ind] == 0) {
                ind[ind_r++] = child_ind;
                que[que_r++] = child;
                vis[child_ind] = 1;
            }
        }
    }
}

// 转换函数 : 将ATCG转换为0123
void transStr(char *str) {
    for (int j = 0 ; str[j] ; ++j) {
        switch(str[j]) {
            case 'A' :
                str[j] = '0';
                break;
            case 'T' :
                str[j] = '1';
                break;
            case 'C' :
                str[j] = '2';
                break;
            case 'G' :
                str[j] = '3';
                break;
        }
    }
}

int main() {
    int n, L;
    char str[200];
    while (scanf("%d%d", &n, &L) != EOF) {
        Trie *root = new_node();
        node_cnt = 1;
        for (int i = 0 ; i < n ; ++i) {
            getchar();
            scanf("%s", str);
            transStr(str);
            insert(root, str);
        }
        build_automaton(root);
    
        // 设置矩阵大小
        mat_size = node_cnt * 3;
        init_unint_mat();

        BFS(root);
        CGMatrix = quick_mat_pow(CGMatrix, L);

        int ans = 0;
        for (int j = 2 * node_cnt ; j < mat_size ; ++j) {
            ans += CGMatrix.m[0][j];
            ans %= MOD;
        }
        printf("%d
", ans);
        memset(CGMatrix.m, 0, sizeof(CGMatrix.m));
        clear(root);
    }
    return 0;
}