原题题面
加里敦大学的生物研究所，发现了决定人喜不喜欢吃藕的基因序列
S
S
S,有这个序列的碱基序列就会表现出喜欢吃藕的性状，但是研究人员发现对碱基序列
S
S
S，任意修改其中不超过 3 个碱基，依然能够表现出吃藕的性状。现在研究人员想知道这个基因在 DNA 链
S
0
S_0
S0上的位置。所以你需要统计在一个表现出吃藕性状的人的 DNA 序列
S
0
S_0
S0上，有多少个连续子串可能是该基因，即有多少个
S
0
S_0
S0的连续子串修改小于等于三个字母能够变成
S
S
S。
输入格式
第一行有一个整数
T
T
T，表示有几组数据。 每组数据第一行一个长度不超过
1
0
5
10^5
105的碱基序列
S
0
S_0
S0。 每组数据第二行一个长度不超过
1
0
5
10^5
105的吃藕基因序列
S
S
S。
输出格式
共
T
T
T行，第
i
i
i行表示第
i
i
i组数据中，在
S
0
S_0
S0中有多少个与
S
S
S等长的连续子串可能是表现吃藕性状的碱基序列。
输入样例
1 ATCGCCCTA CTTCA
输出样例
2
题面分析
这本是一道后缀自动机的题目，但也有题解用了FFT，本篇博客采用了FFT的方法。 我们记第一个串为
s
1
s1
s1，第二个串为
s
2
s2
s2。 对于字符集
{
A
,
C
,
G
,
T
}
\{A,C,G,T\}
{A,C,G,T}，我们对于每个字符依次处理，并记匹配到了该字符时为1，未匹配到该字符时为0，得到两个数组
A
,
B
A,B
A,B。 那么对于
s
1
s1
s1串中的
[
k
,
k
+
l
e
n
(
s
2
)
]
[k,k+len(s2)]
[k,k+len(s2)]和
s
2
s2
s2的
[
0
,
l
e
n
(
s
2
)
−
1
]
[0,len(s2)-1]
[0,len(s2)−1]可以成功匹配的位数为
∑
i
=
0
l
e
n
(
s
2
)
−
1
A
[
k
+
i
]
B
[
i
]
\sum_{i=0}^{len(s2)-1}{A[k+i]B[i]}
i=0∑len(s2)−1A[k+i]B[i] 这个格式长得就很像卷积，我们设
C
[
i
]
=
B
[
l
e
n
(
s
2
)
−
i
−
1
]
C[i]=B[len(s2)-i-1]
C[i]=B[len(s2)−i−1],得到原式等价于
∑
i
=
0
l
e
n
(
s
2
)
−
1
A
[
k
+
i
]
C
[
l
e
n
(
s
2
)
−
i
−
1
]
\sum_{i=0}^{len(s2)-1}{A[k+i]C[len(s2)-i-1]}
i=0∑len(s2)−1A[k+i]C[len(s2)−i−1] 即
∑
i
+
j
=
l
e
n
(
s
2
)
+
k
−
1
A
[
i
]
C
[
j
]
\sum_{i+j=len(s2)+k-1}{A[i]C[j]}
i+j=len(s2)+k−1∑A[i]C[j] 用FFT算出四个字符的结果后相加，再统计
[
l
e
n
(
s
2
)
−
1
,
l
e
n
(
s
1
)
]
[len(s2)-1,len(s1)]
[len(s2)−1,len(s1)]的位置的结果就是符合题意的串个数即可。 假设我们记FFt的结果中某个符合条件的为
a
n
s
[
i
]
ans[i]
ans[i],那么
a
n
s
[
i
]
+
3
≥
l
e
n
(
s
2
)
ans[i]+3\geq len(s2)
ans[i]+3≥len(s2)。
AC代码(900ms)
#include <bits/stdc++.h>
using namespace std
;
#define ll long long
#define cp complex<double>
const int maxn
=4e5;
cp a
[maxn
+10];
cp b
[maxn
+10];
char s1
[maxn
+10];
char s2
[maxn
+10];
ll ans
[maxn
+10];
char charset
[]={'A', 'C', 'G', 'T'};
namespace FFT
{
const double pi
=acos(-1.0);
ll rev
[maxn
+10];
int getBit(int k
)
{
int bit
=0;
while((1<<bit
)<k
)
bit
++;
return bit
;
}
void solve(cp
*a
, int n
, int inv
)
{
int bit
=getBit(n
);
for(int i
=0; i
<n
; ++i
)
{
rev
[i
]=(rev
[i
>>1]>>1)|((i
&1)<<(bit
-1));
if (i
<rev
[i
]) swap(a
[i
], a
[rev
[i
]]);
}
for(int mid
=1; mid
<n
; mid
*=2)
{
cp
w(cos(pi
*1.0/mid
), inv
*sin(pi
*1.0/mid
));
for(int i
=0; i
<n
; i
+=mid
*2)
{
cp
omega(1, 0);
for(int j
=0; j
<mid
; j
++, omega
*=w
)
{
cp x
=a
[i
+j
];
cp y
=omega
*a
[i
+j
+mid
];
a
[i
+j
]=x
+y
;
a
[i
+j
+mid
]=x
-y
;
}
}
}
if (inv
==-1)
{
for(int i
=0; i
<n
; i
++)
{
a
[i
].real(a
[i
].real()/n
);
}
}
}
}
void solve()
{
int t
;
scanf("%d", &t
);
while(t
--)
{
scanf("%s%s", s1
, s2
);
int len1
=strlen(s1
);
int len2
=strlen(s2
);
int len
=1<<FFT
::getBit(len1
+len2
);
for(int i
=0; i
<len
; ++i
)
ans
[i
]=0;
for(int z
=0; z
<4; ++z
)
{
for(int i
=0; i
<len1
; ++i
)
{
a
[i
].real(s1
[i
]==charset
[z
] ? 1 : 0);
a
[i
].imag(0);
}
for(int i
=len1
;i
<len
;++i
)
{
a
[i
]=cp(0,0);
}
for(int i
=0; i
<len2
; ++i
)
{
b
[len2
-i
-1]=(s2
[i
]==charset
[z
] ? 1 : 0);
b
[i
].imag(0);
}
for(int i
=len2
;i
<len
;++i
)
{
b
[i
]=cp(0,0);
}
FFT
::solve(a
, len
, 1);
FFT
::solve(b
, len
, 1);
for(int i
=0; i
<len
; ++i
)
{
a
[i
]*=b
[i
];
}
FFT
::solve(a
, len
, -1);
for(int i
=0; i
<len
; ++i
)
{
ans
[i
]+=(ll
) (a
[i
].real()+0.5);
}
}
ll cnt
=0;
for(int i
=len2
-1; i
<len1
; ++i
)
{
if (ans
[i
]+3>=len2
)
cnt
++;
}
printf("%lld\n", cnt
);
}
}
signed main()
{
#ifdef ACM_LOCAL
freopen("in.txt", "r", stdin);
freopen("out.txt", "w", stdout);
long long test_index_for_debug
=1;
char acm_local_for_debug
;
while(cin
>>acm_local_for_debug
)
{
cin
.putback(acm_local_for_debug
);
if (test_index_for_debug
>100)
{
throw runtime_error("Check the stdin!!!");
}
auto start_clock_for_debug
=clock();
solve();
auto end_clock_for_debug
=clock();
cout
<<"\nTest "<<test_index_for_debug
<<" successful"<<endl
;
cerr
<<"Test "<<test_index_for_debug
++<<" Run Time: "
<<double(end_clock_for_debug
-start_clock_for_debug
)/CLOCKS_PER_SEC
<<"s"<<endl
;
cout
<<"--------------------------------------------------"<<endl
;
}
#else
solve();
#endif
return 0;
}
后记
FFT居然可以处理部分字符串问题，绝了… DrGilbert 2020.10.5
