FFT

https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

 

 

 

#include <cstdio>
#include <cmath>
#define Maxn 1350000
using namespace std;
const double Pi=acos(-1);
int n,m;
struct CP
{
  CP (double xx=0,double yy=0){x=xx,y=yy;}
  double x,y;
  CP operator + (CP const &B) const
  {return CP(x+B.x,y+B.y);}
  CP operator - (CP const &B) const
  {return CP(x-B.x,y-B.y);}
  CP operator * (CP const &B) const
  {return CP(x*B.x-y*B.y,x*B.y+y*B.x);}
  //除法没用
}f[Maxn<<1],sav[Maxn<<1];
void dft(CP *f,int len)
{
  if (len==1)return ;//边界
  //指针的使用比较巧妙 
  CP *fl=f,*fr=f+len/2;
  for (int k=0;k<len;k++)
       sav[k]=f[k];
  for (int k=0;k<len/2;k++)//分奇偶打乱
    {
	    fl[k]=sav[k<<1];
		fr[k]=sav[k<<1|1];
	}
  dft(fl,len/2);
  dft(fr,len/2);//处理子问题
  //由于每次使用的单位根次数不同(len次单位根),所以要重新求。
  CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0);
  for (int k=0;k<len/2;k++)
  {
    //这里buf = (len次单位根的第k个) 
    sav[k]=fl[k]+buf*fr[k];//(1)
    sav[k+len/2]=fl[k]-buf*fr[k];//(2)
    //这两条语句具体见Tips的式子
    buf=buf*tG;//得到下一个单位根。
  }
  for (int k=0;k<len;k++)
       f[k]=sav[k];
}
int main()
{
  scanf("%d",&n);
  for (int i=0;i<n;i++)
      scanf("%lf",&f[i].x);
  //一开始都是实数,虚部为0
  for(m=1;m<n;m<<=1);
  //把长度补到2的幂,不必担心高次项的系数,因为默认为0
  dft(f,m);
  for(int i=0;i<m;++i)
    printf("(%.4f,%.4f)\n",f[i].x,f[i].y);
  return 0;
}

　　

#include <cstdio>
#include <cmath>
#define Maxn 1350000
using namespace std;
const double Pi=acos(-1);
int n,m;
struct CP
{
  CP (double xx=0,double yy=0){x=xx,y=yy;}
  double x,y;
  CP operator + (CP const &B) const
  {return CP(x+B.x,y+B.y);}
  CP operator - (CP const &B) const
  {return CP(x-B.x,y-B.y);}
  CP operator * (CP const &B) const
  {return CP(x*B.x-y*B.y,x*B.y+y*B.x);}
  //除法没用
}f[Maxn<<1],p[Maxn<<1],sav[Maxn<<1];
void dft(CP *f,int len)
{
  if (len==1)return ;//边界
  //指针的使用比较巧妙 
  CP *fl=f,*fr=f+len/2;
  for (int k=0;k<len;k++)sav[k]=f[k];
  for (int k=0;k<len/2;k++)//分奇偶打乱
    {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];}
  dft(fl,len/2);
  dft(fr,len/2);//处理子问题
  //由于每次使用的单位根次数不同(len次单位根),所以要重新求。
  CP tG(cos(2*Pi/len),sin(2*Pi/len)),buf(1,0);
  for (int k=0;k<len/2;k++){
    //这里buf = (len次单位根的第k个) 
    sav[k]=fl[k]+buf*fr[k];//(1)
    sav[k+len/2]=fl[k]-buf*fr[k];//(2)
    //这两条语句具体见Tips的式子
    buf=buf*tG;//得到下一个单位根。
  }for (int k=0;k<len;k++)f[k]=sav[k];
}
void idft(CP *f,int len)
{
  if (len==1)return ;//边界
  //指针的使用比较巧妙 
  CP *fl=f,*fr=f+len/2;
  for (int k=0;k<len;k++)sav[k]=f[k];
  for (int k=0;k<len/2;k++)//分奇偶打乱
    {fl[k]=sav[k<<1];fr[k]=sav[k<<1|1];}
  idft(fl,len/2);
  idft(fr,len/2);//处理子问题
  CP tG(cos(2*Pi/len),  -  sin(2*Pi/len)),buf(1,0);
               //注意这 ↑ 个负号! 
  for (int k=0;k<len/2;k++)
  {
    //这里buf = (len次反单位根的第k个) 
    sav[k]=fl[k]+buf*fr[k];
    sav[k+len/2]=fl[k]-buf*fr[k];
    buf=buf*tG;//得到下一个反单位根。
  }
  for (int k=0;k<len;k++)
       f[k]=sav[k];
}
int main()
{
  scanf("%d%d",&n,&m);
  for (int i=0;i<=n;i++)
       scanf("%lf",&f[i].x);
  for (int i=0;i<=m;i++)
       scanf("%lf",&p[i].x);
  for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂
  dft(f,n);
  dft(p,n);//DFT
  for(int i=0;i<n;++i)
      f[i]=f[i]*p[i];//点值直接乘
  idft(f,n);//IDFT
  for(int i=0;i<=m;++i)
      printf("%d ",(int)(f[i].x/n+0.49));
  //注意结果除以n
  return 0;
}

　　

 

 

三次变两次"优化
根据 (a+bi)*(c+di)=a(c+di)*bi(c+di)=ac-bd+adi+bci
假设我们需要求 F(x)*G(x)
设复多项式 P(x)=F(x)+G(x)i
也就是实部为 F(x),虚部为 G(x).
则 P(x)^2=(F(x)+G(x)i)^2=F(x)^2-G(x)^2+2F(x)G(x)i

其实部为F(x)^2-G(x)^2,
其虚部为2F(x)G(x)

也就是说求出 P(x)^2之后,把它的虚部除以2即可.

#include<algorithm>
#include<cstdio>
#include<cmath>
#define Maxn 1350000
using namespace std;
const double Pi=acos(-1);
inline int read()
{
  register char ch=0;
  while(ch<48||ch>57)ch=getchar();
  return ch-'0';
}
int n,m;
struct CP
{
  CP (double xx=0,double yy=0){x=xx,y=yy;}
  double x,y;
  CP operator + (CP const &B) const
  {return CP(x+B.x,y+B.y);}
  CP operator - (CP const &B) const
  {return CP(x-B.x,y-B.y);}
  CP operator * (CP const &B) const
  {return CP(x*B.x-y*B.y,x*B.y+y*B.x);}
}f[Maxn<<1];//只用了一个复数数组 
int tr[Maxn<<1];
void fft(CP *f,bool flag)
{
  for (int i=0;i<n;i++)
    if (i<tr[i])swap(f[i],f[tr[i]]);
  for(int p=2;p<=n;p<<=1){
    int len=p>>1;
    CP tG(cos(2*Pi/p),sin(2*Pi/p));
    if(!flag)tG.y*=-1;
    for(int k=0;k<n;k+=p){
      CP buf(1,0);
      for(int l=k;l<k+len;l++){
        CP tt=buf*f[len+l];
        f[len+l]=f[l]-tt;
        f[l]=f[l]+tt;
        buf=buf*tG;
      }
    }
  }
}
int main()
{
  scanf("%d%d",&n,&m);
  for (int i=0;i<=n;i++)
      f[i].x=read(); //实部  
  for (int i=0;i<=m;i++)
      f[i].y=read(); //虚部 
  for(m+=n,n=1;n<=m;n<<=1);
  for(int i=0;i<n;i++)
       tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
  fft(f,1); //求出离散的值 
  for(int i=0;i<n;++i) //进行平方 
     f[i]=f[i]*f[i];
  fft(f,0);  //反演出系数 
  for(int i=0;i<=m;++i)  //虚部 /2，即为所要求的系数 
    printf("%d ",(int)(f[i].y/n/2+0.49));
  return 0;
}