转载请注明出处,谢谢http://blog.csdn.net/bigtiao097?viewmode=contents
参考文章http://blog.csdn.net/gatevin/article/details/46892383
题意:
给两个字符串a、b以及数k,求字符串a中与b Hamming距离小于等于k的不同子串的个数
(Hamming距离的定义为两个相同长度的字符串对应位置字符不同的位置数量,如“abbab”和”bbabb”的Hamming距离是3因为有3个字符不一样)
思路:
大概做过几道FFT的题,但是没想出来这个怎么用FFT,最后看了别人的博客学会的,确实用的巧妙
我们构造多项式,字符‘a’代表系数为0,字符‘b’代表系数为1,然后将字符串‘b’反着写,将两个多项式相乘,那么只有当两个字符串都是‘b’(因为我们记‘b’为1)时多项式多对应那一项结果的系数会+1,我们举一个例子:
字符串A:aabbab
字符串B:ab
多项式A:001101(左边代表低位,右边代表高位)
多项式B:010000—-000010(反着写)
两者相乘可得:000000011010
第七位代表A的子串‘aa’与B相同的字母个数,第八位代表A的子串‘ab’与B相同的字母个数,依此类推……
然后我们反着来一遍(‘a’代表系数为1,‘b’代表系数为0)
多项式A:110010(左边代表低位,右边代表高位)
多项式B:100000—-000001(反着写)
两者相乘可得:000000110010
含义同上
感觉非常巧妙,有点那种n+1 = 1+n=2+(n-1) + 3+(n-2)=……的感觉
这样我们枚举A的每一个子串,可以O(1)的求出该子串与B中不同的字母个数(子串长度-相同‘a’的个数-相同‘b’的个数),就知道了Hamming距离是不是小于K,对于不同子串这个要求,Hash一下利用set统计一下就可以了
Result: Accepted Time: 2155ms
具体代码如下:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const double PI = acos(-1);
const ll base = 163;
const int maxn = 1e5+5;//maxn = max(len1,len2)
const int maxm = 4*maxn;//k是>=maxn的最小的2的幂,maxm=2*k
int xx[maxn], yy[maxn], len1, len2,len;//len1,len2分别为两个多项式的最高次数+1
int num[maxm];
int temp[maxn];
int suma[maxm];
int sumb[maxm];
char a[maxn];
char b[maxn];
int lena,lenb;
int k;
int t;
ll ha[maxn];
ll p[maxn];
set<ll> s;
void init()
{
ha[0] = 0;
int n = strlen(a+1);
for(int i = 1; i <= n; i++)
ha[i] = ha[i-1]*base+(a[i]-'a'+1);
}
ll get(int l, int r)
{
return ha[r]-ha[l-1]*p[r-l+1];
}
//复数结构体
struct Complex
{
double x, y;//实部为x,虚部为y
Complex(double x=0, double y=0):x(x),y(y){}
Complex operator+(const Complex &rhs) const { return Complex(x+rhs.x, y+rhs.y);}
Complex operator-(const Complex &rhs) const { return Complex(x-rhs.x, y-rhs.y);}
Complex operator*(const Complex &rhs) const { return Complex(x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x);}
};
Complex x1[maxm],x2[maxm];
/*
进行FFT和IFFT前的反转变换。
将位置i和(i二进制反转后位置)互换。
len必须取2的幂
*/
void change(Complex *x, int len)
{
Complex t;
for(int i = 1, j = len/2; i < len-1; i++)
{
//交换下标互反的元素,i<j保证交换一次
//i做正常的+1,j做反转的+1,始终保持i和j是反转的
if(i < j)
{
t = x[i];
x[i] = x[j];
x[j] = t;
}
int k = len / 2;
while(j >= k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}
/*
做FFT
len必须为2的幂
on==1时是DFT,on==-1时是IDFT
*/
void fft(Complex *x, int len, int on)
{
change(x, len);
for(int h = 2; h <= len; h <<= 1)
{
Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
for(int j = 0; j < len; j += h)
{
Complex w(1, 0);
for(int k = j; k < j+h/2; k++)
{
Complex u = x[k];
Complex t = w*x[k+h/2];
x[k] = u+t;
x[k+h/2] = u-t;
w = w*wn;
}
}
}
if(on == -1)
for(int i = 0; i < len; i++) x[i].x /= len;
}
void workFFT()
{
len = 1;
memset(x1,0,sizeof x1);
memset(x2,0,sizeof x2);
memset(num,0,sizeof num);
while(len < len1*2 || len < len2*2) len <<=1;
for(int i = 0; i < len; i++)
x1[i] = x2[i] = Complex(0,0);
for(int i = 0; i < len1; i++) x1[i] = Complex(xx[i],0);
for(int i = 0; i < len2; i++) x2[i] = Complex(yy[i],0);
fft(x1,len,1);
fft(x2,len,1);
for(int i = 0; i < len; i++) x1[i] = x1[i]*x2[i];
fft(x1,len,-1);
for(int i = 0; i < len; i++) num[i] = (int)(x1[i].x+0.5);
}
int main()
{
p[0] = 1;
for(int i = 1; i <maxn; i++)
p[i] =p[i-1] * base;
while(scanf("%d",&k)&&k!=-1)
{
scanf("%s",a+1);
scanf("%s",b+1);
init();
s.clear();
memset(suma,0,sizeof suma);
memset(sumb,0,sizeof sumb);
memset(temp,0,sizeof temp);
memset(xx,0,sizeof xx);
memset(yy,0,sizeof yy);
lena = strlen(a+1);
lenb = strlen(b+1);
len1 = lena+1;
len2 = len1;
for(int i=1;i<=lena;i++)
xx[i] = a[i]-'a';
for(int i=1;i<=lenb;i++)
temp[i] = b[i]-'a';
for(int i=1;i<len2;i++)
yy[i] = temp[len1-i];
workFFT();
for(int i=0;i<len;i++)
sumb[i] = num[i];
for(int i=1;i<=lena;i++)
xx[i] = 1-xx[i];
for(int i=1;i<=lenb;i++)
temp[i] = 1-temp[i];
for(int i=1;i<len2;i++)
yy[i] = temp[len1-i];
workFFT();
for(int i=0;i<len;i++)
suma[i] = num[i];
for(int i=lenb;i<=lena;i++)
if(lenb-suma[lena+i-lenb+1]-sumb[lena+i-lenb+1]<=k)
s.insert(get(i-lenb+1,i));
printf("Case %d: %d\n",++t,s.size());
}
}