[AC自动机+矩阵快速幂] poj2778 与 hdu2243

W_Zifan

于 2021-04-01 15:31:18 发布

阅读量292

点赞数 2

CC 4.0 BY-SA版权

文章标签：算法数据结构字符串矩阵

本文链接：https://blog.youkuaiyun.com/ECNU_SF/article/details/115373898

本文介绍了如何使用AC自动机解决两类构造字符串的问题：一类是构造不含特定子串的字符串；另一类是构造包含特定词根的字符串。文章通过具体题目详细讲解了AC自动机的构建过程及矩阵快速幂的应用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

poj 2778link

题面

给定了 $m$ 个 A, C, G, T 构成的基因片段(字符串)，让你构造一个不含这 $m$ 个片段的，长度为 $n$ 的由A, C, G, T构成的字符串，问有几种方式。

分析

这题刚开始以为有关动态规划哈希什么的，看了题解才发现这是AC自动机的构造神题。AC自动机的基本介绍网上有很多，这里就不赘述。
首先对于这 $m$ 个字符串，可以建一个 AC 自动机，就如普通的AC自动机一样，维护 fail 数组，用 $t r [i] [j]$ 表示第 $i$ 个状态读入字符 $j$ 之后变成的状态。然后对于每一个状态，需要记录这个状态是否包含一个完整的基因片段，其中 $e [i]$ 表示从根到状态 $i$ 的字符串是否为基因片段，然后用 $l s t [i]$ 表示从根到状态 $i$ 的字符串的某一个后缀是否为基因片段。
这两个数组可以在建立自动机的时候维护，我们可以知道，若是到了某个状态， $\ || \ lst[i] != 0$ ，这个状态就包含了一个基因片段，所以不可进入这个状态(不可达)。

之后，我们建立一个邻接矩阵A，其中 $A [i] [j]$ 表示从 $i$ 到 $j$ 的转移方式(每个转移方式对应一个字符)，但是要注意不可达状态的行和列要置为0。这样 $A [i] [j]$ 就表示从状态 $i$ 到状态 $j$ 长度为 1 的合法字符串个数， $A^{n}[i][j]$ 就表示从状态 $i$ 到状态 $j$ 长度为 n 的合法字符串个数，我们所要求的就是 $\sum A^{n}[0][j]$ 。
对于 $A^{n}[i][j]$ 则需要用矩阵快速幂来计算，由于poj卡的比较紧，所以取mod的地方要在保证正确的时候尽可能减少。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<stack>
#include<queue>
#include<map>
#include<set>
#include<vector>
#include<cmath>
#include<string>
#include<sstream>
#include<bitset>
#define here printf("modassing [%s] in LINE %d\n", __FUNCTION__, __LINE__);
#define debug(x) cout << #x << ":\t" << (x) << endl;

using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> P;
const int maxn = 1000 + 10;
const int maxm = 2e6 + 10;
const int INF = 0x3f3f3f3f;
const ll mod = 1e9 + 7;
const double pi = acos(-1.0);
const double eps = 0.0000001;

int N;

struct mt                            //N表示方阵的维数
{
    ll a[102][102];
}mat;

mt mult(mt A, mt B, ll mod)          //计算方阵A * B
{
    mt res;
    for(int i = 0; i < N; i++)
    {
        for(int j = 0; j < N; j++)
        {
            res.a[i][j]=0;
            for(int k = 0; k < N; k++)
            {
                res.a[i][j] += A.a[i][k] * B.a[k][j];
            }
            res.a[i][j] %= mod;
        }
    }
    return res;
}
 
mt power(mt a,ll b, ll mod)         //计算方阵A^b
{
    mt res;
    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++)
            res.a[i][j] = 0;
    for(int i = 0; i < N; i++) res.a[i][i] = 1;   //res为单位矩阵
 
    while(b)
    {
        if(b & 1) res = mult(res, a, mod);
        b >>= 1;
        a = mult(a, a, mod);
    }
    return res;
}


//tot记录trie树结点个数，0是根结点
int tr[maxn][4], tot;
int e[maxn], fail[maxn], lst[maxn];          //e记录有几个这样的字符串

int getch(char ch)
{
    if(ch == 'A')
        return 0;
    if(ch == 'C')
        return 1;
    if(ch == 'G')
        return 2;
    return 3;
}

void insertString(char *s)         //字符串从1开始
{
    int u = 0;
    for (int i = 1; s[i]; i++)
    {
        if (!tr[u][getch(s[i])]) tr[u][getch(s[i])] = ++tot;
        u = tr[u][getch(s[i])];
    }
    e[u] = 1;
}

void build()
{
    queue<int> q;
    for (int i = 0; i < 4; i++)
        if (tr[0][i]) q.push(tr[0][i]);    //先插入深度为1的结点
    
    while(q.size())
    {
        int u = q.front();
        q.pop();
        for (int i = 0; i < 4; i++)
        {
            if(tr[u][i]) //有i+'a'的子结点，则让它的fail指针指向 它父亲的fail指针指向的那个结点 的i+'a'结点
            {
                fail[tr[u][i]] = tr[fail[u]][i];
                q.push(tr[u][i]);
                lst[tr[u][i]] = (e[fail[tr[u][i]]] ? fail[tr[u][i]] : lst[fail[tr[u][i]]]);	 //lst指向存储了完整串的fail节点
            }
            else tr[u][i] = tr[fail[u]][i];  //若是没有子结点，让它指向 它父亲的fail指针指向的那个结点 的i+'a'结点
        }
    }
}

int m;
ll n;
char s[110];

int main()
{
    scanf("%d %lld", &m, &n);
    for (int i = 1; i <= m; i++)
    {
        scanf("%s", s + 1);
        insertString(s);
    }
    build();

    for (int i = 0; i <= tot; i++)
    {
        if(e[i] || lst[i]) continue;
        for (int j = 0; j < 4; j++)
        {
            int tmp = tr[i][j];
            if(e[tmp] || lst[tmp]) continue;
            mat.a[i][tmp]++;
        }
    }
    N = tot + 1;
    mt ans = power(mat, n, 100000);
    ll res = 0;
    for (int i = 0; i <= tot; i++)
        res = (res + ans.a[0][i]) % 100000;
    printf("%lld\n", res);
}

hdu 2243link

题面

也是给你一些字符串表示词根，让你构造一个答案不超过 $L$ 的字符串，其中这个字符串至少包括一个词根。
要求答案对 $2^{64}$ 取模，即自然溢出。

分析

其实若是求恰好为 $L$ 的字符串，就是poj那题的补集。若是不超过，答案就是 (不超过 $L$ 的字符串 - 不超过 $L$ 且不包含词根的字符串)。

不超过 $L$ 的字符串，数量总共有 $26^1 + 26^2 +... + 26^L$ 个，若是用倍增公式的话，有一个除数，取模的数不是素数，并且很大，不太好取逆元，快速幂是一种求和的方法：
令 $a_k = 26 ^ k, \ s_k = \sum_{i = 1}^{k}a_k$ ，则 $s_k = s_{i-k} + a_k$ ，写成矩阵的形式：
$\left[ \begin{matrix} 26 & 0\\ 1 & 1\\ \end{matrix} \right] \left[ \begin{matrix} a_i\\ s_{i-1}\\ \end{matrix} \right] = \left[ \begin{matrix} a_{i+1}\\ s_{i}\\ \end{matrix} \right]\\$
于是就有：
$\left[ \begin{matrix} a_{n+1}\\ s_{n}\\ \end{matrix} \right] = \left[ \begin{matrix} 26 & 0\\ 1 & 1\\ \end{matrix} \right]^n \left[ \begin{matrix} a_1\\ s_{0}\\ \end{matrix} \right] = \left[ \begin{matrix} 26 & 0\\ 1 & 1\\ \end{matrix} \right]^n \left[ \begin{matrix} 26\\ 0\\ \end{matrix} \right]$

不超过 $L$ 且不包含词根的字符串，若用poj那题 $A [i] [j]$ 的定义，数量总共有 $A^1 + A^2 +... + A^L$ 个。
令 $S_k = \sum_{i = 1}^{k}A_k$ ，则 $S_k = S_{i-k} + A_k$ ，写成矩阵的形式：
$\left[ \begin{matrix} A & 0\\ 1 & 1\\ \end{matrix} \right] \left[ \begin{matrix} A_i\\ S_{i-1}\\ \end{matrix} \right] = \left[ \begin{matrix} A_{i+1}\\ S_{i}\\ \end{matrix} \right]\\$
于是就有：
$\left[ \begin{matrix} A_{n+1}\\ S_{n}\\ \end{matrix} \right] = \left[ \begin{matrix} A & 0\\ 1 & 1\\ \end{matrix} \right]^n \left[ \begin{matrix} A_1\\ S_{0}\\ \end{matrix} \right] = \left[ \begin{matrix} A & 0\\ 1 & 1\\ \end{matrix} \right]^n \left[ \begin{matrix} A\\ 0\\ \end{matrix} \right]$
值得注意的是，这里的 1 指的不是全 1 矩阵，而是单位矩阵。还有多组数据记得清空。

#include <bits/stdc++.h>
#define here printf("modassing [%s] in LINE %d\n", __FUNCTION__, __LINE__);
#define debug(x) cout << #x << ":\t" << (x) << endl;

using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> P;
const int maxn = 1000 + 10;
const int maxm = 2e6 + 10;
const int INF = 0x3f3f3f3f;
const ll mod = 1e9 + 7;
const double pi = acos(-1.0);
const double eps = 0.0000001;

ull qpow(ull a, ll b)
{
    ull t = 1;
    for(; b; b >>= 1, a = a * a)
        if(b & 1)  t = t * a;
    return t;
}

int N;

struct mt                            //N表示方阵的维数
{
    ull a[70][70];
}mat, M;

mt mult(mt A, mt B)          //计算方阵A * B
{
    mt res;
    for(int i = 0; i < N; i++)
    {
        for(int j = 0; j < N; j++)
        {
            res.a[i][j]=0;
            for(int k = 0; k < N; k++)
            {
                res.a[i][j] += A.a[i][k] * B.a[k][j];
            }
        }
    }
    return res;
}
 
mt power(mt a,ll b)         //计算方阵A^b
{
    mt res;
    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++)
            res.a[i][j] = 0;
    for(int i = 0; i < N; i++) res.a[i][i] = 1;   //res为单位矩阵
 
    while(b)
    {
        if(b & 1) res = mult(res, a);
        b >>= 1;
        a = mult(a, a);
    }
    return res;
}


//tot记录trie树结点个数，0是根结点
int tr[maxn][26], tot;
int e[maxn], fail[maxn], lst[maxn];          //e记录有几个这样的字符串

void insertString(char *s)         //字符串从1开始
{
    int u = 0;
    for (int i = 1; s[i]; i++)
    {
        if (!tr[u][s[i]-'a']) tr[u][s[i]-'a'] = ++tot;
        u = tr[u][s[i]-'a'];
    }
    e[u] = 1;
}

void build()
{
    queue<int> q;
    for (int i = 0; i < 26; i++)
        if (tr[0][i]) q.push(tr[0][i]);    //先插入深度为1的结点
    
    while(q.size())
    {
        int u = q.front();
        q.pop();
        for (int i = 0; i <26; i++)
        {
            if(tr[u][i]) //有i+'a'的子结点，则让它的fail指针指向 它父亲的fail指针指向的那个结点 的i+'a'结点
            {
                fail[tr[u][i]] = tr[fail[u]][i];
                q.push(tr[u][i]);
                lst[tr[u][i]] = (e[fail[tr[u][i]]] ? fail[tr[u][i]] : lst[fail[tr[u][i]]]);	 //lst指向存储了完整串的fail节点
            }
            else tr[u][i] = tr[fail[u]][i];  //若是没有子结点，让它指向 它父亲的fail指针指向的那个结点 的i+'a'结点
        }
    }
}

void clear()
{
    memset(mat.a, 0, sizeof(mat.a));
    memset(M.a, 0, sizeof(M.a));
    for (int i = 0; i <= tot; i++)
    {
        e[i] = fail[i] = lst[i] = 0;
        for (int j = 0; j < 26; j++)
            tr[i][j] = 0;
    }
    tot = 0;
}


int m;
ll n;
char s[110];

int main()
{
    while(~scanf("%d %lld", &m, &n))
    {
        N = 2;
        mt g1;
        g1.a[0][0] = 26;
        g1.a[0][1] = 0;
        g1.a[1][0] = g1.a[1][1] = 1;
        mt g2 = power(g1, n);
        ull res = g2.a[1][0] * 26;

        for (int i = 1; i <= m; i++)        
        {
            scanf("%s", s + 1);
            insertString(s);
        }
        build();                    //构造AC自动机

        for (int i = 0; i <= tot; i++)            //构造A
        {
            if(e[i] || lst[i]) continue;
            for (int j = 0; j < 26; j++)
            {
                int tmp = tr[i][j];
                if(e[tmp] || lst[tmp]) continue;
                mat.a[i][tmp]++;
            }
        }

        N = 2 * tot + 2; 
        for (int i = 0; i < N; i++)            //构造M
        {
            if(i <= tot)
                for (int j = 0; j <= tot; j++)
                    M.a[i][j] = mat.a[i][j];
            else
                for (int j = 0; j < N; j++)
                    if(j == i - (tot+1) || j == i) M.a[i][j] = 1;
        }

        mt ans = power(M, n);

        N = tot + 1;
        mt tmp;
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                tmp.a[i][j] = ans.a[i + N][j];
        ans = mult(mat, tmp);

        for (int i = 0; i < N; i++)
            res -=ans.a[0][i];
        printf("%I64u\n", res);
        clear();
    }
}