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}}对应着函数在这个位置的取值。
假设任意两个不同的xj都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：
{\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; }
