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=si−k+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][aisi−1]=[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=Si−k+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][AiSi−1]=[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();
}
}