FFT 学习记录

https://zhuanlan.zhihu.com/p/19763231?refer=wille

http://blog.jobbole.com/58246/

http://blog.csdn.net/yeeman/article/details/6325693

https://www.zhihu.com/question/22298352

http://blog.csdn.net/iamzky/article/details/22712347#comments

之前一直都不理解，今天在度教(tlzmybm)的教导下，应该对于大部分地方都懂了。

拉格朗日插值法:

对某个多项式函数，已知有给定的k + 1个取值点：

{\displaystyle (x_{0},y_{0}),\ldots ,(x_{k},y_{k})}

其中{\displaystyle x_{j}}对应着自变量的位置，而{\displaystyle y_{j}}对应着函数在这个位置的取值。

假设任意两个不同的x_j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：

{\displaystyle L(x):=\sum _{j=0}^{k}y_{j}\ell _{j}(x)}

其中每个{\displaystyle \ell _{j}(x)}为拉格朗日基本多项式（或称插值基函数），其表达式为：

{\displaystyle \ell _{j}(x):=\prod _{i=0,\,i\neq j}^{k}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{j-1})}{(x_{j}-x_{j-1})}}{\frac {(x-x_{j+1})}{(x_{j}-x_{j+1})}}\cdots {\frac {(x-x_{k})}{(x_{j}-x_{k})}}.}

拉格朗日基本多项式{\displaystyle \ell _{j}(x)}的特点是在{\displaystyle x_{j}}上取值为1，在其它的点{\displaystyle x_{i},\,i\neq j}上取值为0。

范例：

假设有某个二次多项式函数{\displaystyle f}，已知它在三个点上的取值为：

• {\displaystyle f(4)=10}
• {\displaystyle f(5)=5.25}
• {\displaystyle f(6)=1}

要求{\displaystyle f(18)}的值。

首先写出每个拉格朗日基本多项式：

{\displaystyle \ell _{0}(x)={\frac {(x-5)(x-6)}{(4-5)(4-6)}}}

{\displaystyle \ell _{1}(x)={\frac {(x-4)(x-6)}{(5-4)(5-6)}}}

{\displaystyle \ell _{2}(x)={\frac {(x-4)(x-5)}{(6-4)(6-5)}}}

然后应用拉格朗日插值法，就可以得到{\displaystyle p}的表达式（{\displaystyle p}为函数{\displaystyle f}的插值函数）：

{\displaystyle p(x)=f(4)\ell _{0}(x)+f(5)\ell _{1}(x)+f(6)\ell _{2}(x)}

{\displaystyle .\,\,\,\,\,\,\,\,\,\,=10\cdot {\frac {(x-5)(x-6)}{(4-5)(4-6)}}+5.25\cdot {\frac {(x-4)(x-6)}{(5-4)(5-6)}}+1\cdot {\frac {(x-4)(x-5)}{(6-4)(6-5)}}}
{\displaystyle .\,\,\,\,\,\,\,\,\,\,={\frac {1}{4}}(x^{2}-28x+136)}

此时代入数值{\displaystyle \ 18}就可以求出所需之值：{\displaystyle \ f(18)=p(18)=-11}。

/*     algorithm : High-Precision FFT  */  #include <cstdio>  #include <cstring>  #include <cmath>  #include <algorithm>  #define N 200005  #define pi acos(-1.0) // PI值  using namespace std;  struct complex  {      double r,i;      complex(double real=0.0,double image=0.0){          r=real; i=image;      }      // 以下为三种虚数运算的定义      complex operator + (const complex o){          return complex(r+o.r,i+o.i);      }      complex operator - (const complex o){          return complex(r-o.r,i-o.i);      }      complex operator * (const complex o){          return complex(r*o.r-i*o.i,r*o.i+i*o.r);      }  }x1[N],x2[N];  char a[N/2],b[N/2];  int sum[N]; // 结果存在sum里  void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)  {      register int i,j,k;      for(i=1,j=l/2;i<l-1;i++)      {          if(i<j)  swap(y[i],y[j]); // 交换互为下标反转的元素                                  // i<j保证只交换一次          k=l/2;          while(j>=k) // 由最高位检索，遇1变0，遇0变1，跳出          {              j-=k;              k/=2;          }          if(j<k)  j+=k;      }  }  void fft(complex *y,int l,double on) // FFT O(nlogn)                              // 其中on==1时为DFT，on==-1为IDFT  {      register int h,i,j,k;      complex u,t;       brc(y,l); // 调用反转置换      for(h=2;h<=l;h<<=1) // 控制层数      {          // 初始化单位复根          complex wn(cos(on*2*pi/h),sin(on*2*pi/h));          for(j=0;j<l;j+=h) // 控制起始下标          {              complex w(1,0); // 初始化螺旋因子              for(k=j;k<j+h/2;k++) // 配对              {                  u=y[k];                  t=w*y[k+h/2];                  y[k]=u+t;                  y[k+h/2]=u-t;                  w=w*wn; // 更新螺旋因子              } // 据说上面的操作叫蝴蝶操作…          }      }      if(on==-1)  for(i=0;i<l;i++) y[i].r/=l; // IDFT  }  int main(void)  {      int l1,l2,l;      register int i;      while(scanf("%s%s",a,b)!=EOF)      {          l1=strlen(a);          l2=strlen(b);          l=1;          while(l<l1*2 || l<l2*2)   l<<=1; // 将次数界变成2^n                                          // 配合二分与反转置换          for(i=0;i<l1;i++) // 倒置存入          {              x1[i].r=a[l1-i-1]-'0';              x1[i].i=0.0;          }          for(;i<l;i++)    x1[i].r=x1[i].i=0.0;          // 将多余次数界初始化为0          for(i=0;i<l2;i++)          {              x2[i].r=b[l2-i-1]-'0';              x2[i].i=0.0;          }          for(;i<l;i++)    x2[i].r=x2[i].i=0.0;          fft(x1,l,1); // DFT(a)          fft(x2,l,1); // DFT(b)          for(i=0;i<l;i++) x1[i]=x1[i]*x2[i]; // 点乘结果存入a          fft(x1,l,-1); // IDFT(a*b)          for(i=0;i<l;i++) sum[i]=x1[i].r+0.5; // 四舍五入          for(i=0;i<l;i++) // 进位          {              sum[i+1]+=sum[i]/10;              sum[i]%=10;          }          l=l1+l2-1;          while(sum[l]<=0 && l>0)   l--; // 检索最高位          for(i=l;i>=0;i--)    putchar(sum[i]+'0'); // 倒序输出          putchar('\n');      }      return 0;  }