• 51nod1565 FFT


    思路:

    显然拆位FFT 不解释

    //By SiriusRen
    #include <bits/stdc++.h>
    using namespace std;
    const int N=550000;
    const double pi=acos(-1);
    int n=1,L,S,T,k,sa[N],st[N],sc[N],sg[N],R[N],csa,cst,csc,csg,ans;
    char s[N],t[N],A[N];
    struct Complex{double x,y;Complex(){}Complex(double X,double Y){x=X,y=Y;}}Sa[N],St[N],Sc[N],Sg[N],Ta[N],Tt[N],Tc[N],Tg[N];
    Complex operator+(Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
    Complex operator-(Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
    Complex operator*(Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
    Complex operator/(Complex a,int b){return Complex(a.x/b,a.y/b);}
    void FFT(Complex *a,int f){
        for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);
        for(int l=1;l<n;l<<=1){
            Complex wn=Complex(cos(pi/l),f*sin(pi/l));
            for(int j=0;j<n;j+=(l<<1)){
                Complex w=Complex(1,0);
                for(int k=0;k<l;k++,w=w*wn){
                    Complex X=a[j+k],Y=w*a[j+k+l];
                    a[j+k]=X+Y,a[j+k+l]=X-Y;
                }
            }
        }if(f==-1)for(int i=0;i<n;i++)a[i]=a[i]/n;
    }
    int main(){
        scanf("%d%d%d%s%s",&S,&T,&k,s+1,t+1);
        for(;n<=S+T;n<<=1)L++;
        for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
        for(int i=1;i<=S;i++){
            if(s[i]=='A')sa[max(0,i-k)]++,sa[i+k+1]--;
            else if(s[i]=='T')st[max(0,i-k)]++,st[i+k+1]--;
            else if(s[i]=='C')sc[max(0,i-k)]++,sc[i+k+1]--;
            else if(s[i]=='G')sg[max(0,i-k)]++,sg[i+k+1]--;
        }
        for(int i=1;i<=S;i++)sa[i]+=sa[i-1],st[i]+=st[i-1],sc[i]+=sc[i-1],sg[i]+=sg[i-1];
        for(int i=1;i<=T;i++){
            if(t[i]=='A')Ta[T-i]=Complex(1,0),csa++;
            else if(t[i]=='T')Tt[T-i]=Complex(1,0),cst++;
            else if(t[i]=='C')Tc[T-i]=Complex(1,0),csc++;
            else if(t[i]=='G')Tg[T-i]=Complex(1,0),csg++;
        }
        for(int i=1;i<=S;i++){
            if(sa[i])Sa[i-1]=Complex(1,0);
            if(st[i])St[i-1]=Complex(1,0);
            if(sc[i])Sc[i-1]=Complex(1,0);
            if(sg[i])Sg[i-1]=Complex(1,0);
        }
        FFT(Sa,1),FFT(St,1),FFT(Sc,1),FFT(Sg,1),FFT(Ta,1),FFT(Tt,1),FFT(Tc,1),FFT(Tg,1);
        for(int i=0;i<n;i++)Sa[i]=Sa[i]*Ta[i];
        for(int i=0;i<n;i++)St[i]=St[i]*Tt[i];
        for(int i=0;i<n;i++)Sc[i]=Sc[i]*Tc[i];
        for(int i=0;i<n;i++)Sg[i]=Sg[i]*Tg[i];
        FFT(Sa,-1),FFT(St,-1),FFT(Sc,-1),FFT(Sg,-1);
        for(int i=0;i<S;i++)A[i]=1;
        for(int i=0;i<S;i++)
            if((int)(Sa[i+T-1].x+0.2)!=csa||(int)(St[i+T-1].x+0.2)!=cst||(int)(Sc[i+T-1].x+0.2)!=csc||(int)(Sg[i+T-1].x+0.2)!=csg)A[i]=0;
        for(int i=0;i<S;i++)if(A[i])ans++;
        printf("%d
    ",ans);
    }
  • 相关阅读:
    nginx centos 服务开机启动设置实例详解
    CentOS打开关闭永久防火墙指定端口
    使用 nginx 反向代理 sqlserver 访问 配置
    Springboot集成Mybatis
    linux中查看java进程
    mybatis关于jdbc连接报错,5.5.62MySQL连接,出现com.mysql.cj.jdbc.exceptions.CommunicationsException: Communications link failure等问题解决方法
    索引的分析
    MySQL慢查询分析工具
    MySQL锁
    nGrinder介绍、编写脚本与执行(完整版)
  • 原文地址:https://www.cnblogs.com/SiriusRen/p/7058043.html
Copyright © 2020-2023  润新知