FFT 学习记录

最新推荐文章于 2024-01-21 23:09:14 发布

取竹

最新推荐文章于 2024-01-21 23:09:14 发布

阅读量376

点赞数

分类专栏： FFT

FFT 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

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个取值点：

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

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

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

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

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

\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}})}}.

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

范例：

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

$f(4)=10$
$f(5)=5.25$
$f(6)=1$

要求 $f(18)$ 的值。

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

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

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

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

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

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

.\,\,\,\,\,\,\,\,\,\,=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)}}

.\,\,\,\,\,\,\,\,\,\,={\frac {1}{4}}(x^{2}-28x+136)

此时代入数值 $\ 18$ 就可以求出所需之值： $\ 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;  
}