原题题面
加里敦大学的生物研究所,发现了决定人喜不喜欢吃藕的基因序列 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)
//inv是1时是系数转点值,-1是点值转系数
{
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);
// printf("%lld ", ans[i]);
}
// printf("\n");
}
ll cnt=0;
for(int i=len2-1; i<len1; ++i)
{
if (ans[i]+3>=len2)
cnt++;
}
printf("%lld\n", cnt);
}
}
signed main()
{
// ios_base::sync_with_stdio(false);
// cin.tie(0);
// cout.tie(0);
#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