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

poj 2778link

题面

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

分析

  这题刚开始以为有关动态规划哈希什么的,看了题解才发现这是AC自动机的构造神题。AC自动机的基本介绍网上有很多,这里就不赘述。
  首先对于这 m m m 个字符串,可以建一个 AC 自动机,就如普通的AC自动机一样,维护 fail 数组, t r [ i ] [ j ] tr[i][j] tr[i][j] 表示第 i i i 个状态读入字符 j j j 之后变成的状态。然后对于每一个状态,需要记录这个状态是否包含一个完整的基因片段,其中 e [ i ] e[i] e[i] 表示从根到状态 i i i 的字符串是否为基因片段,然后用 l s t [ i ] lst[i] lst[i] 表示从根到状态 i i i 的字符串的某一个后缀是否为基因片段。
  这两个数组可以在建立自动机的时候维护,我们可以知道,若是到了某个状态, e [ i ] ! = 0   ∣ ∣   l s t [ i ] ! = 0 e[i] != 0 \ || \ lst[i] != 0 e[i]!=0  lst[i]!=0,这个状态就包含了一个基因片段,所以不可进入这个状态(不可达)。
  
  之后,我们建立一个邻接矩阵A,其中 A [ i ] [ j ] A[i][j] A[i][j] 表示从 i i i j j j 的转移方式(每个转移方式对应一个字符),但是要注意不可达状态的行和列要置为0。这样 A [ i ] [ j ] A[i][j] A[i][j] 就表示 从 状态 i i i 到 状态 j j j 长度为 1 的合法字符串个数, A n [ i ] [ j ] A^{n}[i][j] An[i][j] 就表示 从 状态 i i i 到 状态 j j j 长度为 n 的合法字符串个数,我们所要求的就是 ∑ A n [ 0 ] [ j ] \sum A^{n}[0][j] An[0][j]
  对于 A n [ i ] [ j ] A^{n}[i][j] An[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 L L 的字符串,其中这个字符串至少包括一个词根。
  要求答案对 2 64 2^{64} 264 取模,即自然溢出。

分析

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

  不超过 L L L 的字符串,数量总共有 2 6 1 + 2 6 2 + . . . + 2 6 L 26^1 + 26^2 +... + 26^L 261+262+...+26L 个,若是用倍增公式的话,有一个除数,取模的数不是素数,并且很大,不太好取逆元,快速幂是一种求和的方法:
  令 a k = 2 6 k ,   s k = ∑ i = 1 k a k a_k = 26 ^ k, \ s_k = \sum_{i = 1}^{k}a_k ak=26k, sk=i=1kak,则 s k = s i − k + a k s_k = s_{i-k} + a_k sk=sik+ak,写成矩阵的形式:
[ 26 0 1 1 ] [ a i s i − 1 ] = [ a i + 1 s i ] \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]\\ [26101][aisi1]=[ai+1si]
于是就有:
[ a n + 1 s n ] = [ 26 0 1 1 ] n [ a 1 s 0 ] = [ 26 0 1 1 ] n [ 26 0 ] \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] [an+1sn]=[26101]n[a1s0]=[26101]n[260]

  不超过 L L L 且不包含词根的字符串,若用poj那题 A [ i ] [ j ] A[i][j] A[i][j] 的定义,数量总共有 A 1 + A 2 + . . . + A L A^1 + A^2 +... + A^L A1+A2+...+AL 个。
  令 S k = ∑ i = 1 k A k S_k = \sum_{i = 1}^{k}A_k Sk=i=1kAk,则 S k = S i − k + A k S_k = S_{i-k} + A_k Sk=Sik+Ak,写成矩阵的形式:
[ A 0 1 1 ] [ A i S i − 1 ] = [ A i + 1 S i ] \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]\\ [A101][AiSi1]=[Ai+1Si]
于是就有:
[ A n + 1 S n ] = [ A 0 1 1 ] n [ A 1 S 0 ] = [ A 0 1 1 ] n [ A 0 ] \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] [An+1Sn]=[A101]n[A1S0]=[A101]n[A0]
值得注意的是,这里的 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();
    }
}

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值