以前刷过的FFT

Gym - 101667H 

2017-2018 ACM-ICPC, Asia Daejeon Regional Contest

#include<bits/stdc++.h>
using namespace std;
#define maxn 4000005
const double pi=acos(-1.0);
struct com
{
    double x,y;
    com(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
}a[maxn],b[maxn];
com operator + (com a,com b) {return com(a.x+b.x,a.y+b.y);}
com operator - (com a,com b) {return com(a.x-b.x,a.y-b.y);}
com operator * (com a,com b) {return com(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int S,T,n,m,L,R[maxn],ans[maxn];
long long F[maxn];
char s[maxn],t[maxn];
double f[maxn],g[maxn];
void FFT(com a[maxn],int opt)
{
    for (int i=0;i<n;++i)
        if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1)
    {
        com wn=com(cos(pi/k),opt*sin(pi/k));
        for (int i=0;i<n;i+=(k<<1))
        {
            com w=com(1,0);
            for (int j=0;j<k;++j,w=w*wn)
            {
                com x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}
void calc(int opt)
{
    FFT(a,1);FFT(b,1);
    for (int i=0;i<=n;++i) a[i]=a[i]*b[i];
    FFT(a,-1);
    for (int i=S-1;i<T;++i)
    {
         F[i]+=(long long)(a[i].x/n+0.5)*opt;//对于每种匹配位置累计赢的次数
    }
}
int main()
{
    scanf("%d%d",&T,&S);
    scanf("%s%s",t,s);//S短,T长
    for(int i=0;i<S/2;++i) swap(s[i],s[S-i-1]);//短串逆置
    T=T+S;
    for(int i=T-S;i<T;i++)t[i]='B';
    t[T]=0;
    m=S+T-2;
    for(int i=0;i<T;i++)
    {
        if(t[i]=='S')t[i]='R';
        else if(t[i]=='R')t[i]='P';
        else if(t[i]=='P')t[i]='S';
    }
    for (n=1;n<=m;n<<=1) ++L;
    for (int i=0;i<n;++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0;i<T;++i) f[i]=(t[i]=='S');
    for (int i=0;i<S;++i) g[i]=(s[i]=='S');
    for (int i=0;i<=n;++i) a[i]=com(0,0),b[i]=com(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    for (int i=0;i<T;++i) f[i]=(t[i]=='R');
    for (int i=0;i<S;++i) g[i]=(s[i]=='R');
    for (int i=0;i<=n;++i) a[i]=com(0,0),b[i]=com(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    for (int i=0;i<T;++i) f[i]=(t[i]=='P');
    for (int i=0;i<S;++i) g[i]=(s[i]=='P');
    for (int i=0;i<=n;++i) a[i]=com(0,0),b[i]=com(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    long long  ans=0;
    for (int i=S-1;i<T;++i)//这个范围自己考虑一下就好了
      ans=max(ans,F[i]);//所有位置取max
    printf("%lld\n",ans);
}

 模仿这个题目在TOJ出了一个题,还无人AC 烟村四五家

Fuzzy Search

 CodeForces - 528D 

Leonid works for a small and promising start-up that works on decoding the human genome. His duties include solving complex problems of finding certain patterns in long strings consisting of letters 'A', 'T', 'G' and 'C'.

Let's consider the following scenario. There is a fragment of a human DNA chain, recorded as a string S. To analyze the fragment, you need to find all occurrences of string T in a string S. However, the matter is complicated by the fact that the original chain fragment could contain minor mutations, which, however, complicate the task of finding a fragment. Leonid proposed the following approach to solve this problem.

Let's write down integer k ≥ 0 — the error threshold. We will say that string Toccurs in string S on position i (1 ≤ i ≤ |S| - |T| + 1), if after putting string Talong with this position, each character of string T corresponds to the some character of the same value in string S at the distance of at most k. More formally, for any j (1 ≤ j ≤ |T|) there must exist such p (1 ≤ p ≤ |S|), that |(i + j - 1) - p| ≤ k and S[p] = T[j].

For example, corresponding to the given definition, string "ACAT" occurs in string "AGCAATTCAT" in positions 2, 3 and 6.

Note that at k = 0 the given definition transforms to a simple definition of the occurrence of a string in a string.

Help Leonid by calculating in how many positions the given string T occurs in the given string S with the given error threshold.

Input

The first line contains three integers |S|, |T|, k (1 ≤ |T| ≤ |S| ≤ 200 000, 0 ≤ k ≤ 200 000) — the lengths of strings S and T and the error threshold.

The second line contains string S.

The third line contains string T.

Both strings consist only of uppercase letters 'A', 'T', 'G' and 'C'.

Output

Print a single number — the number of occurrences of T in S with the error threshold k by the given definition.

Examples

Input
10 4 1
AGCAATTCAT
ACAT
Output
3

Note

If you happen to know about the structure of the human genome a little more than the author of the problem, and you are not impressed with Leonid's original approach, do not take everything described above seriously.

 

#include<bits/stdc++.h>
using namespace std;
typedef double db;
const db PI=acos(-1.0);
struct Complex
{
    db x,y;
    Complex(db _x=0.0,db _y=0.0):x(_x),y(_y){}
    Complex operator + (const Complex &b)const
    {
        return Complex(x+b.x,y+b.y);
    }
    Complex operator - (const Complex &b)const
    {
        return Complex(x-b.x,y-b.y);
    }
    Complex operator * (const Complex &b)const
    {
        return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
    }
};
void change(Complex y[],int len)
{
    for(int i=1,j=len/2;i<len-1;i++)
    {
        if(i<j)swap(y[i],y[j]);
        int k=len/2;
        while(j>=k)
        {
            j-=k;
            k/=2;
        }
        if(j<k)j+=k;
    }
}
void fft(Complex y[],int len,int on)
{
    change(y,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=y[k];
                Complex t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)
        for(int i=0;i<len;i++)
            y[i].x/=len;
}
const char ch[4]={'A','C','G','T'};
int lens,lent,k;
char s[200005],t[200005];
int es[200005],ps[200005];
Complex ss[600005],tt[600005];
int cnt[4][200005];
void work(int id)
{
    memset(es,0,sizeof(es));
    memset(ps,0,sizeof(ps));
    memset(ss,0,sizeof(ss));
    memset(tt,0,sizeof(tt));
    for(int i=0;i<lens;i++)
        if(s[i]==ch[id])
        {
            es[max(i-k,0)]++;
            es[min(i+k+1,lens)]--;
        }
    ps[0]=es[0];
    for(int i=1;i<lens;i++)ps[i]=es[i]+ps[i-1];
    for(int i=0;i<lens;i++)ss[i].x=(ps[i]>0);
    for(int i=0;i<lent;i++)tt[i].x=(t[i]==ch[id]);
    reverse(tt,tt+lent);
    int len=1;
    while(len<2*lens || len<2*lent)len<<=1;
    fft(ss,len,1);
    fft(tt,len,1);
    for(int i=0;i<len;i++)ss[i]=ss[i]*tt[i];
    fft(ss,len,-1);
    for(int i=0;i<lens;i++)cnt[id][i]=(int)(ss[i+lent-1].x+0.5);
}
int main()
{
    scanf("%d%d%d",&lens,&lent,&k);
    scanf("%s%s",s,t);
    for(int i=0;i<4;i++)work(i);
    int ans=0;
    for(int i=0;i<lens;i++)
    {
        int tot=0;
        for(int j=0;j<4;j++)tot+=cnt[j][i];
        ans+=(tot==lent);
    }
    printf("%d\n",ans);
    return 0;
}

 

posted @ 2018-08-10 20:37  暴力都不会的蒟蒻  阅读(210)  评论(0编辑  收藏  举报