FFT

今天看了算法导论，对FFT感受颇深。

感觉我就在抄算法导论。

回归正题。

多项式的表示

系数表示法

系数表示法其实非常常见，其实就是：

$$A(x)=sumlimits_{j=0}^{n-1}a_jx^j$$

这是一个n次多项式，每个项的次数为0,1,...,n-1

如果用系数表示法来做多项式加法，那么时间复杂度是$O(N)$的

但是用系数表示法来做多项式乘法，那么时间复杂度是$O(N^2)$的，这样非常慢，所以这里介绍的FFT可以把多项式乘法的时间复杂度优化到$O(Nlog_2N)$

点值表示法

 一个n次多项式$A(x)$的点值表示就是n个点值对所形成的集合：

$${(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}$$

其中所有的$x_k$各不相同，并且当k=0,1,...,n-1时有

$$y_k=A(x_k)$$

一个多项式可以有很多种点值表示，只要选取的n个$x_k$互异即可。

对于许多关于多项式的操作来说，点值表示法是很方便的。

如果用点值表示法来做多项式加法，那么时间复杂度是$O(N)$的：

如果$C(x)=A(x)+B(x)$，则对于任意的$x_k$，$C(x_k)=A(x_k)+B(x_k)$。准确地说，如果已知$A(x)$的点集表示：

$${(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}$$

和$B(x)$的点集表示

$${(x_0,y^{'}_0),(x_1,y^{'}_1),...,(x_{n-1},y^{'}_{n-1})}$$

（注意，A和B对相同的n个点求值）则$C(x)$的点集表示为：

$${(x_0,y_0+y^{'}_0),(x_1,y_1+y^{'}_1),...,(x_{n-1},y_{n-1}+y^{'}_{n-1})}$$

如果点值表示法来做多项式乘法，那么时间复杂度也是$O(N)$的：

如果$C(x)=A(x)B(x)$，则对于任意的$x_k$，$C(x_k)=A(x_k)B(x_k)$。准确地说，如果已知$A(x)$的点集表示：

$${(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}$$

和$B(x)$的点集表示

$${(x_0,y^{'}_0),(x_1,y^{'}_1),...,(x_{n-1},y^{'}_{n-1})}$$

（注意，A和B对相同的n个点求值）则$C(x)$的点集表示为：

$${(x_0,y_0y^{'}_0),(x_1,y_1y^{'}_1),...,(x_{n-1},y_{n-1}y^{'}_{n-1})}$$

系数表示法转化成点值表示法称为求值；

点值表示法转化成系数表示法称为插值。

我们发现，对于多项式乘法，点值表示法比系数表示法有优势。

如果我们任选n个点求值（把系数表示法转化成点值表示法），那么时间复杂度是$O(N^2)$，但是稍后我们可以看到，如果我们巧妙地选取$x_k$，那么时间复杂度可以做到$O(Nlog_2N)$，并且插值（把点值表示法转化成系数表示法）也可以做到$O(Nlog_2N)$

于是我们得到一个思路：

 

单位复根

这里为稍后我们讨论奠定了一下基础。

n次单位复根是指满足$w^n=1$的复数$w$。n次单位复根有n个，分别是$e^{2pi ik/n},k=0,1,..,n-1$

我们可以从欧拉公式理解：

$$e^{ix}=cosx+isinx$$

如图我们可以看到，n个单位复根平均分布在以原点为圆心的复平面上，其中

$$w_n=e^{2pi i/n}$$

称为主n次单位复根，任何一个n次单位复根都是主n次单位复根$w_n$的幂

所以n个n次单位复根就是$w_n^0,w_n^1,...,w_n^{n-1}$

所以其实$w_n^k=e^{2pi ik/n},k=0,1,..,n-1$

然后然后n次单位复根有一些性质：

引理一（相消引理）：对于任意整数$ngeq 0$，$kgeq 0$，$d>0$，有$w_{dn}^{dk}=w_n^k$

 证明：$w_{dn}^{dk}=e^{2pi idk/dn}=e^{2pi ik/n}=w_n^k$

引理二：对于偶数$n>0$，有$w_{n}^{n/2}=w_2=-1$

引理三（折半引理）：对于偶数$n>0$，n个n次单位复根的平方等于n/2个n/2次单位复根。

引理四（求和引理）：对于任意整数$ngeq 1$和不能被n整除的非负整数k，有：

$$sumlimits_{j=0}^{n-1}(w_{n}^{k})^j=0$$

证明：等比数列的求和公式也可以用于复数：

$sumlimits_{j=0}^{n-1}(w_{n}^{k})^j=frac{(w_{n}^{k})^{n}-1}{w_{n}^{k}-1}=frac{(w_{n}^{n})^{k}-1}{w_{n}^{k}-1}=frac{(1)^{k}-1}{w_{n}^{k}-1}=0$

 从系数表示法到点值表示法

不失一般性，我们设n是2的幂，如果n不是2的幂，可以适当增大n。

前面已经提到，其实我们是对n个n次单位复根$w_n^0,w_n^1,...,w_n^{n-1}$进行求值，即：

如果

$$A(x)=sumlimits_{j=0}^{n-1}a_jx^j$$

那么

$$y_k=A(w_n^{k})=sumlimits_{j=0}^{n-1}a_jw_n^{kj},其中k=0,1,...,n-1$$

我们运用分治策略，用$A(x)$中偶数下标的系数和奇数下标的系数，分别定义了两个新的n/2次多项式$A^{[0]}(x)$和$A^{[1]}(x)$：

$A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}$

$A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1}$

那么有如下式：

$$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$$

我们发现，$A^{[0]}(x)$和$A^{[1]}(x)$是一个子问题，我们递归处理。

现在我们已经知道了$A^{[0]}(x)$在$w_{n/2}^0,w_{n/2}^1,...,w_{n/2}^{n/2-1}$处的取值：

$${(w_{n/2}^0,y^{'}_{0}),(w_{n/2}^1,y^{'}_{1}),...,(w_{n/2-1}^{n/2-1},y^{'}_{n/2-1})}$$

$$其中y^{'}_{k}=A^{[0]}(w_{n/2}^k)=A^{[0]}(w_{n}^{2k})，k=0,1,...,n/2-1$$

和$A^{[1]}(x)$在$w_{n/2}^0,w_{n/2}^1,...,w_{n/2}^{n/2-1}$处的取值：

$${(w_{n/2}^0,y^{''}_{0}),(w_{n/2}^1,y^{''}_{1}),...,(w_{n/2-1}^{n/2-1},y^{''}_{n/2-1})}$$

$$其中y^{''}_{k}=A^{[1]}(w_{n/2}^k)=A^{[1]}(w_{n}^{2k})，k=0,1,...,n/2-1$$

 我们想知道$A(x)$在$w_n^0,w_n^1,...,w_n^{n-1}$处的取值的求值：

$${(w_{n/2}^0,y_{0}),(w_{n/2}^1,y_{1}),...,(w_{n-1}^{n-1},y_{n-1})}$$

运算一下：

$若k=0,1,...,n/2-1$

$y_{k}$
$=A(w_n^{k})$
$=A^{[0]}((w_n^{k})^2)+w_n^{k}A^{[1]}((w_n^{k})^2)$
$=A^{[0]}(w_n^{2k})+w_n^{k}A^{[1]}(w_n^{2k})$
$=y^{'}_{k}+w_n^{k}y^{''}_{k}$

$y_{k+n/2}$
$=A(w_n^{k+n/2})$
$=A^{[0]}((w_n^{k+n/2})^2)+w_n^{k+n/2}A^{[1]}((w_n^{k+n/2})^2)$
$=A^{[0]}(w_n^{2k+n})+w_n^{k+n/2}A^{[1]}(w_n^{2k+n})$
$=A^{[0]}(w_n^{2k})-w_n^{k}A^{[1]}(w_n^{2k})$

$=y^{'}_{k}-w_n^{k}y^{''}_{k}$

（其实就是上几行的等式，不懂可以往回看几行）

我们发现我们已经知道了$y^{'}_{k}$和$y^{''}_{k}$了！！！！！

所以我们完全可以算出$y_{0},y_{1},...y_{n-1}$了。

好了，到现在为止可以在$O(Nlog_2N)$的时间内完成从系数表示法到点值表示法的转化。

伪代码是这样的：

 从点值表示法到系数表示法

我们通过矩阵来理解：

从左到右3个矩阵分别记为$Y$，$V$和$A$，其中$Y_{j,1}=y_j(j=0,1,...,n-1)$，$V_{i,j}=w_n^{ij}(i=0,1,...,n-1，j=0,1,...,n-1)$

我们可以表示为矩阵积：$Y=VA$

现在我们已经知道了$Y$和$V$，想求$A$

如果我们知道$V$的逆$V^{-1}$，那么我们可以求$A$：

$Y=VA$

$V^{-1}Y=V^{-1}VA$

$V^{-1}Y=A$

我们先给出一个结论：$V^{-1}_{i,j}=w_n^{-ij}/n(i=0,1,...,n-1，j=0,1,...,n-1)$

证明：

记$X=V^{-1}V$

$X_{i,j}=sumlimits_{k=0}^{n-1}V^{-1}_{i,k}V_{k,j}=sumlimits_{k=0}^{n-1}frac{w_n^{-ik}w_n^{kj}}{n}=sumlimits_{k=0}^{n-1}frac{(w_n^{j-i})^k}{n}$

若$i=j$，那么$w_n^{j-i}=1$，所以$X_{i,j}=1$

若$i eq j$，根据上面单位复根的求和引理，可知$X_{i,j}=0$

$X$是一个单位矩阵

所以$V^{-1}$是$V$的逆。

根据$V^{-1}Y=A$我们可以得到：

$a_k=sumlimits_{j=0}^{n-1}V^{-1}_{k,j}y_j=sumlimits_{j=0}^{n-1}w_n^{-kj}y_j/n=frac{1}{n}sumlimits_{j=0}^{n-1}y_jw_n^{-kj}$

即：

$$a_k=frac{1}{n}sumlimits_{j=0}^{n-1}y_jw_n^{-kj},其中k=0,1,...,n-1$$

我们来看一下”从系数表示法到点值表示法“时我们是怎么求值的：

$$y_k=sumlimits_{j=0}^{n-1}a_jw_n^{kj},其中k=0,1,...,n-1$$

我们发现这两条式子惊人地相似，就是$y_k$换成了$a_k$，$a_j$换成了$y_j$，$w_n^{kj}$换成了$w_n^{-kj}$，还多了了一个$frac{1}{n}$，不过这个很好处理。

我们完全可以沿用求值的函数就可以，将($y_0,y_1,...,y_{n-1})$传递进函数，然后$w_n^{kj}$变成$w_n^{-kj}$即可；最后还要除以$n$

多项式乘法的程序是这样的：

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<fstream>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<utility>
#include<set>
#include<bitset>
#include<vector>
#include<functional>
#include<deque>
#include<cctype>
#include<climits>
#include<complex>

using namespace std;

typedef long long LL;
typedef double DB;
typedef pair<int,int> PII;
typedef complex<DB> CP;

#define mmst(a,v) memset(a,v,sizeof(a))
#define mmcy(a,b) memcpy(a,b,sizeof(a))
#define re(i,a,b)  for(i=a;i<=b;i++)
#define red(i,a,b) for(i=a;i>=b;i--)
#define fi first
#define se second
#define m_p(a,b) make_pair(a,b)
#define SF scanf
#define PF printf
#define two(k) (1<<(k))

template<class T>inline T sqr(T x){return x*x;}
template<class T>inline void upmin(T &t,T tmp){if(t>tmp)t=tmp;}
template<class T>inline void upmax(T &t,T tmp){if(t<tmp)t=tmp;}

const DB EPS=1e-9;
inline int sgn(DB x){if(abs(x)<EPS)return 0;return(x>0)?1:-1;}
const DB Pi=acos(-1.0);

inline int gint()
  {
        int res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }
inline LL gll()
  {
      LL res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }

/*
注意事项
1.数组要开大 保证高位有足够的0
2.n必须是2的幂次，补0即可
3.最后所以的数要除n
4.n次多项式与m次多项式的乘积为n+m次多项式
  长度为n的多项式与长度为m的多项式的乘积为长度为n+m-1的多项式 
*/
 
const int maxN=100000*4;

int alen,blen,clen,N;
CP a[maxN+1000],b[maxN+1000],c[maxN+1000];

CP t[maxN+1000];
inline void FFT(CP *a,int N,int l,int r,int flag)
  {
      if(l==r) return;
      int i,len=r-l+1,mid=(l+r)/2,ln=l,rn=mid+1;
      for(i=l;i<=r;i+=2)t[ln++]=a[i],t[rn++]=a[i+1];
      re(i,l,r)a[i]=t[i];
      FFT(a,N,l,mid,flag);FFT(a,N,mid+1,r,flag);
      CP wn(cos(2*Pi/DB(len)),sin(DB(flag)*2*Pi/DB(len))),w(1.0,0);//注意没有=号，常错 
      re(i,l,mid){CP x=a[i],t=w*a[i+len/2];a[i]=x+y;a[i+len/2]=x-y;w=w*wn;}
      if(flag==-1 && len==N)re(i,l,r)a[i]/=DB(N);
  }

int main()
  {
      freopen("UOJ#34.in","r",stdin);
      freopen("UOJ#34.out","w",stdout);
      int i;
      alen=gint()+1;blen=gint()+1;
      re(i,0,alen-1)a[i]=DB(gint());
      re(i,0,blen-1)b[i]=DB(gint());
      clen=alen+blen-1;
      for(N=1;N<clen;N<<=1);
      FFT(a,N,0,N-1,1);
      FFT(b,N,0,N-1,1);
      re(i,0,N-1)c[i]=a[i]*b[i];
      FFT(c,N,0,N-1,-1);
      re(i,0,clen-1)printf("%d ",int(c[i].real()+0.5));printf("
");
      return 0;
  }

最后我们得到一个定理：

模运算下的FFT

上面我们谈论FFT都是用实数的。

但有时候我们要将系数$a_i$模一个数$P$，准确地说：

$已知A(x)和B(x):$

$$A(x)=sumlimits_{j=0}^{n-1}a_jx^j，其中0leq a_j<P$$

$$B(x)=sumlimits_{j=0}^{n-1}b_jx^j，其中0leq b_j<P$$

$A(x)和B(x)的积为C(x)=A(x)B(x),那么$

$$c_j=sumlimits_{k=0}^{j}a_kb_{j-k}\%P，其中0leq j leq2n-2$$

那么怎么做呢？

用原根（如果不知道原根，可以看我另一篇博客）

假设$P$是质数，那么$P$一定有原根（不妨设为$G$）。

设$P=2^{a} imes b+1,其中b为奇数$

我们知道：

$G^{varphi (P)}equiv 1(\% P)$

由费马小定理得：

$G^{P-1}equiv 1(\% P)$

又因为$P=2^{a} imes b+1$，所以$P-1=2^{a} imes b$，即：

$G^{2^{a} imes b}equiv 1(\% P)$

$(G^{b})^{2^{a}}equiv 1(\% P)$

类似于复数的做法，我们重新定义$w_n=(G^{b})^{2^{a}/n}\%P，其中n是2的幂，且1 leq nleq 2^{a}$

那么$w_n^{k}=w_n=(G^{b})^{k*2^{a}/n}\%P，其中k=0,1,...,n-1$

我们发现此时任然满足n次单位复根的性质：

引理一（相消引理）：对于任意整数$ngeq 0$，$kgeq 0$，$d>0$，有$w_{dn}^{dk}=w_n^k$

引理二：对于偶数$n>0$，有$w_{n}^{n/2}=w_2=-1$

证明：

我们先考虑一个问题：

$求解x^2equiv 1(\%P)，其中P是质数，且0leq x< P$

变形得：

$(x^2-1)\%P=0$

$(x+1)(x-1)\%P=0$

因为$P$是质数

所以$x-1=0$或$x+1=P$，即$x=1或P-1$

所以$w_2=G^{b*2^{a-1}}$要么是$1$，要么是$P-1$

又因为$G$是原根

所以$G^0,G^1,...,G^{P-2}$和$1,2,...,P-1$是一一对应的。

显然$G^0$对应了$1$

所以$G^{b*2^{a-1}}$不能对应$1$，只能对应$P-1$

所以$w_2=G^{b*2^{a-1}}=P-1$

又因为在模$P$意义下，$-1$和$P-1$是一样的

所以$w_{n}^{n/2}=w_2=-1$

引理三（折半引理）：对于偶数$n>0$，$n$个$w_n^{k}(k=0,1,...,n-1)$的平方等于$n/2$个$w_{n/2}^{k}(k=0,1,...,n/2-1)$

引理四（求和引理）：对于任意整数$ngeq 1$和不能被n整除的非负整数k，有：

$$sumlimits_{j=0}^{n-1}(w_{n}^{k})^j=0$$

我们发现居然有和n次单位复根的性质一模一样的性质。

我们只需要改动一下即可。

因此，若满足$2^{k}|varphi(p)且2^{k}geq n$,我们可以用$g^{ frac{varphi(p)}{2^k}}$作为单位根，其中$g$是$p$的原根

如果$p$没有原根呢？我们可以将$p$分解质因数，取若干个便于FFT的大质数进行取模，然后再用中国剩余定理还原系数即可。

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<fstream>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<utility>
#include<set>
#include<bitset>
#include<vector>
#include<functional>
#include<deque>
#include<cctype>
#include<climits>
#include<complex>
//#include<bits/stdc++.h>适用于CF,UOJ，但不适用于poj
 
using namespace std;

typedef long long LL;
typedef double DB;
typedef pair<int,int> PII;
typedef complex<DB> CP;

#define mmst(a,v) memset(a,v,sizeof(a))
#define mmcy(a,b) memcpy(a,b,sizeof(a))
#define fill(a,l,r,v) fill(a+l,a+r+1,v)
#define re(i,a,b)  for(i=(a);i<=(b);i++)
#define red(i,a,b) for(i=(a);i>=(b);i--)
#define ire(i,x) for(typedef(x.begin()) i=x.begin();i!=x.end();i++)
#define fi first
#define se second
#define m_p(a,b) make_pair(a,b)
#define p_b(a) push_back(a)
#define SF scanf
#define PF printf
#define two(k) (1<<(k))

template<class T> T sqr(T x){return x*x;}
template<class T> void upmin(T &t,T tmp){if(t>tmp)t=tmp;}
template<class T> void upmax(T &t,T tmp){if(t<tmp)t=tmp;}

const DB EPS=1e-9;
int sgn(DB x){if(abs(x)<EPS)return 0;return(x>0)?1:-1;}
const DB Pi=acos(-1.0);

int gint()
  {
        int res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }
LL gll()
  {
      LL res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }

const int maxm=8000;
const LL MOD=1004535809;
const LL G=3;

int n,m,X,cnt,d;

LL power(LL a,LL k,LL Mod){LL x=1;while(k){if(k&1)x=x*a%Mod;a=a*a%Mod;k>>=1;}return x;}

struct GF
  {
      LL a[4*maxm+200];
      LL& operator [](int i){return a[i];}
      GF(){mmst(a,0);}
      void FFT(LL *a,int n,int type)
        {
            if(n==1)return;
            static LL t[4*maxm+200];
            int i,ln=0,rn=n>>1;
            for(i=0;i<=n-1;i+=2)t[ln++]=a[i],t[rn++]=a[i+1];
            re(i,0,n-1)a[i]=t[i];
            LL *l=a,*r=a+(n>>1);
            FFT(l,n>>1,type);
            FFT(r,n>>1,type);
            LL wn=power(G,(LL)type*(MOD-1)/n%(MOD-1),MOD),w=1;
            for(i=0;i<n>>1;i++,w=w*wn%MOD)t[i]=l[i]+w*r[i],t[i+(n>>1)]=l[i]-w*r[i];
            re(i,0,n-1)a[i]=(t[i]%MOD+MOD)%MOD;
        }
      GF& operator *=(GF &f)
        {
            static LL b[4*maxm+200];
            int i;
            re(i,0,d-1)b[i]=f[i];
            FFT(a,d,1);
            FFT(b,d,1);
            re(i,0,d-1)a[i]=a[i]*b[i]%MOD;
            FFT(a,d,MOD-2);
                re(i,m-1,(m-2)<<1)(a[i-(m-1)]+=a[i])%=MOD,a[i]=0;
            LL inv=power(d,MOD-2,MOD);
            re(i,0,m-2)a[i]=a[i]*inv%MOD;
            return *this;
        }
  };