zoukankan      html  css  js  c++  java
  • FFT

      void FFT(complex a[],int n,int fl){
          for (int i=1,j=n/2;i<n;i++){
            if (i<j) {complex t=a[i];a[i]=a[j];a[j]=t;};
            int k;
          for (k=(n>>1);j&k;j^=k,k>>=1);
          j^=k;
        }
        
         for (int i=2;i<=n;i<<=1){
          complex w;w.r=cos(fl*2*pi/i);w.i=sin(fl*2*pi/i);
          for (int j=0;j<n;j+=i){
              complex wi;wi.r=1;wi.i=0;
              for (int k=j;k<j+i/2;k++){
                complex u=a[k],v=a[k+i/2]*wi;
              a[k]=u+v;a[k+i/2]=u-v;
              wi=wi*w;    
            }
          }    
        }
        
        if (fl==-1) for (int i=0;i<n;i++) a[i].r/=n;
      }

    fl==1求点值 fl==-1插值

    _________________________
    DFT

    void DFT(){
      for (int i=0;i<n;i++){
          cp t=(cp){cos(2*pi/n*i),sin(2*pi/n*i)},bas=(cp){1,0};
          for (int j=0;j<n;j++){
            tmp[i]=tmp[i]+bas*a[j];
          bas=bas*t;    
        }
      }
      for (int i=0;i<n;i++) a[i]=tmp[i];
    }
    
    void IDFT(){
      for (int i=0;i<n;i++) P[i].r=P[i].i=0;
      for (int i=0;i<n;i++){
          cp t=(cp){cos(-2*pi/n*i),sin(-2*pi/n*i)},bas=(cp){1,0};
        for (int j=0;j<n;j++){
          P[i]=P[i]+tmp[j]*bas;
          bas=bas*t;     
        }
      }    
      
      for (int i=0;i<n;i++) P[i].r=(int)(P[i].r/n+0.5);
    }

     -------------------------------------------------------------------------

    CODECHEF JUNE15 MOREFB

    系数为多项式的分治FFT

    #include <bits/stdc++.h>
    #define LL long long 
    #define LDB long double
    using namespace std;
    
      const LDB pi=acos(-1);
      const LL mo=99991;
    
      int a[1000001],n,k;
    
      struct data{
          LL a_1,a_n;
      };
      
      inline data operator*(data a,data b) {return((data){(a.a_n*b.a_n+a.a_1*b.a_1)%mo,(a.a_n*b.a_n+a.a_1*b.a_n+a.a_n*b.a_1)%mo});};
    
      struct cp{
          LDB r_1,r_n,r_n2,i_1,i_n,i_n2;
          
          void set0(){
            r_1=r_n=r_n2=i_1=i_n=i_n2=0;
        }
        
        void set(data a){
          set0();
          r_1=a.a_1;r_n=a.a_n;
        }
      }A[1000001],B[1000001];
      
      inline cp operator *(cp a,cp b){
          return((cp)
          {
              a.r_1*b.r_1-a.i_1*b.i_1,
              a.r_n*b.r_1+a.r_1*b.r_n-a.i_n*b.i_1-a.i_1*b.i_n,
              a.r_n2*b.r_1+a.r_n*b.r_n+a.r_1*b.r_n2-a.i_n2*b.i_1-a.i_n*b.i_n-a.i_1*b.i_n2,
              a.r_1*b.i_1+a.i_1*b.r_1,
              a.r_n*b.i_1+a.r_1*b.i_n+a.i_n*b.r_1+a.i_1*b.r_n,
              a.r_n2*b.i_1+a.r_n*b.i_n+a.r_1*b.i_n2+a.i_n2*b.r_1+a.i_n*b.r_n+a.i_1*b.r_n2
          });
      }
      
      inline cp operator +(cp a,cp b){
          return((cp)
          {
              a.r_1+b.r_1,
              a.r_n+b.r_n,
              a.r_n2+b.r_n2,
              a.i_1+b.i_1,
              a.i_n+b.i_n,
              a.i_n2+b.i_n2,
          });
      }
      
      inline cp operator -(cp a,cp b){
          return((cp)
          {
              a.r_1-b.r_1,
              a.r_n-b.r_n,
              a.r_n2-b.r_n2,
              a.i_1-b.i_1,
              a.i_n-b.i_n,
              a.i_n2-b.i_n2,
          });      
      }
    
      data getkth(int po){
          data ret;ret.a_1=1;ret.a_n=0;
          data bas;bas.a_1=0;bas.a_n=1;
          for (;po;bas=bas*bas){
            if (po&1) ret=ret*bas;
          po>>=1;    
        }
        return(ret);
      }
    
       void FFT(cp a[],int n,int fl){
        for (int i=1,j=n/2;i<n;i++){
          if (i<j) {cp t=a[i];a[i]=a[j];a[j]=t;};
          int k;
          for (k=(n>>1);j&k;j^=k,k>>=1);
          j^=k;
        }
        
         for (int i=2;i<=n;i<<=1){
          cp w;w.set0();
          w.r_1=cos(fl*2*pi/i);w.i_1=sin(fl*2*pi/i);
          for (int j=0;j<n;j+=i){
            cp wi;wi.set0();
            wi.r_1=1;wi.i_1=0;
            for (int k=j;k<j+i/2;k++){
              cp u=a[k],v=a[k+i/2]*wi;
              a[k]=u+v;a[k+i/2]=u-v;
              wi=wi*w;    
            }
          }    
        }
        
        if (fl==-1) for (int i=0;i<n;i++) a[i].r_n2/=n,a[i].r_n/=n,a[i].r_1/=n;
      }
      vector <data> MUL(vector <data> &a,vector <data> &b){
          int n=1;
          while (n<a.size()+b.size()) n<<=1;
          for (int i=0;i<n;i++){
            if (i<a.size()) A[i].set(a[i]);else A[i].set0();
          if (i<b.size()) B[i].set(b[i]);else B[i].set0();    
        }
        
        FFT(A,n,1);FFT(B,n,1);
        for (int i=0;i<n;i++) A[i]=A[i]*B[i];
        FFT(A,n,-1);
        
        vector <data> ret;ret.resize(a.size()+b.size());
        for (int i=0;i<a.size()+b.size();i++){
          LL tn2=(A[i].r_n2+0.5),tn=(A[i].r_n+0.5),t1=(A[i].r_1+0.5);
          ret[i].a_n=(tn2+tn)%mo;ret[i].a_1=(tn2+t1)%mo;
        }
        return(ret);
      }
    
      void show(vector <data> tmp1,vector <data> tmp2){
          for (int i=0;i<tmp1.size();i++) printf("%lld %lld|",tmp1[i].a_n,tmp1[i].a_1);printf("
    ");
          for (int i=0;i<tmp2.size();i++) printf("%lld %lld|",tmp2[i].a_n,tmp2[i].a_1);printf("
    ");
          vector <data> res=MUL(tmp1,tmp2);
          for (int i=0;i<res.size();i++) printf("%lld %lld|",res[i].a_n,res[i].a_1);printf("
    
    ");
      }
    
      vector <data> solve(int l,int r){
          if (l==r){
            vector <data> ret;ret.resize(2);
          ret[0].a_1=1;ret[0].a_n=0;
          ret[1]=getkth(a[l]);
          return(ret);    
        }
        
        int mid=(l+r)>>1;
        vector <data> tmp1=solve(l,mid);
        vector <data> tmp2=solve(mid+1,r);
    //    show(tmp1,tmp2);
        return(MUL(tmp1,tmp2));
      }
    
      int main(){
          scanf("%d%d",&n,&k);
          for (int i=1;i<=n;i++) scanf("%d",&a[i]);
          vector <data> ans=solve(1,n);
          printf("%lld
    ",ans[k].a_n%mo); 
      }
  • 相关阅读:
    MySQL学习记录
    Python3玩转儿 机器学习(4)
    Python3玩转儿 机器学习(3)
    C#-WebForm-文件上传-FileUpload控件
    C#将WebBowser控件替换为谷歌内核【转】
    各种【icon】矢量图
    WPF 获取鼠标全局坐标【精简】
    winfrom 的 各种效果【需要新浪帐号查看】
    JS 文字波纹效果【插件】
    C# 窗口抖动
  • 原文地址:https://www.cnblogs.com/zhujiangning/p/6158114.html
Copyright © 2011-2022 走看看