首页 > 其他 > 详细

芝士:FFT

时间:2019-12-18 21:13:52      阅读:77      评论:0      收藏:0      [点我收藏+]

FFT

背景

对于两个多项式乘法\(O(n^2)\)的时间复杂度难以令人满意

所以将其优化,得到了FFT算法,其时间复杂度为优秀的\(O(n*log_n)\)

前置芝士

系数表示法

对于一个\(n+1\)项的\(n\)次多 ♂ 项式\(f(x)=\sum_{i=0}^{n} a_i*x^i\)

也可如此表示\(f(x)=\{a_0,a_1,a_2,·····,a_{n}\}\)

点值表示法

对于一个函数,很明显我们可以将其与平面直角坐标系相联系

我们回忆一下

对于一次函数,有了两个点就能确定唯一一个的一次函数

对于二次函数,有了三个点就能确定唯一一条二次函数

那么对于一个\(n\)次函数,有了\(n+1\)个点,就能确定一条唯一的\(n\)次函数

所以\(f(x)=\{\{x_1,y_1\},\{x_2,y_2\},·····,\{x_{n+1},y_{n+1}\}\}\)

复数

首先我们知道\(i=\sqrt{-1}\)

形如\(a+b*i\)的形式成为虚数

其中\(a\)被称为实部,\(b*i\)成为虚部

那么运算呢?

\(x=a_1+b_1*i,y=a_2+b_2*i\)

\(\begin{cases}x+y=(a_1+a_2)+(b_1+b_2)*i\\x*y=(a_1*a_2-b_1*b_2)+(a_1*b_2+a_2*b_1)*i\end{cases}\)

进入正题

以下所有讨论的N为\(2^k\)

首先我们要明白我们要的是什么,

对于一个函数\(c(x)=a(x)*b(x)\)

我们先假设\(n为a的项数,m为b的项数\)

我们要求的就是c的表达式

结合上面所讲,如果要将c的表达式求出,

我们必须要求出\(n+m+1\)\(c(x)\)的对应值

并且带入的x必须互不相同。

但是,如果我们知道\(a(x)\)所对应的\((n+m+1)\)个值和\(b(x)\)所对应的\((n+m+1)\)的值,

我们就可以\(O(n+m+1)\)的时间复杂度求出\(c\)的表达式

时间复杂度的瓶颈就在于转换上面

FFT就是解决这个问题

首先讨论\(x\)

我们首先需要\(x\)满足以下性质

1.这n个数不同

2.\(x_n^k=x_{2*n}^{2*k}\)。n表示将一个单位圆分成n份,取其中的k份

3.\(x_n^k=-x_{n}^{k+\frac{n}{2}}\)

4.\(x_n^a*x_n^b=x_n^{a+n}\)

我们将以上满足性质的\(x\)称为\(\omega\)也就是单位复根

我们将\(a(x)\)分成两个函数\(a_0(x)\)\(a_1(x)\)

我们将\(a(x)\)函数中的偶数项的\(a_i\)提取出来,
成为\(a_0(x)\)的每一项,但是函数中的x依然是每次递增1

那么\(a(x)\)可以怎么表示呢?

很明显

\(a(x)=a_0(x^2)+x*a_1(x^2)\)

你是否明白了些什么?

也就是我们如果明白可以一个\(O(1)\)算的函数

我们就可以将这个函数往上递推,

以一个\(O(n)\)的复杂度递推

这就是时间复杂度\(O(n*log_n)\)的来源,也是FFT的过程

常数优化(真的是常数)

非递归

真的如标题所示,将递归的写成非递归,就是将数组进行划分就行了

蝴蝶优化

\(\omega\)写成一个常量,同时新定义一个变量作为中介变量就行了

板子

#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
#define maxn 10000005
#define x first
#define y second
const double pi=acos(-1.0);
int n,m;
int limit=1;
pair<double,double> a[maxn];
pair<double,double> b[maxn];
int l;
int r[maxn];
pair<double,double> operator + (pair<double,double> a,pair<double,double> b)
{
    return make_pair(a.x+b.x,a.y+b.y);
}
pair<double,double> operator - (pair<double,double> a,pair<double,double> b)
{
    return make_pair(a.x-b.x,a.y-b.y);
}
pair<double,double> operator * (pair<double,double> a,pair<double,double> b)
{
    return make_pair(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
#undef x
#undef y
void fft(pair<double,double> *a,int t)
{
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(a[i],a[r[i]]);
    for(int mid=1;mid<limit;mid<<=1)
    {
        pair<double,double> wn=make_pair(cos(pi/mid),t*sin(pi/mid));
        for(int r=mid<<1,j=0;j<limit;j+=r)
        {
            pair<double,double> w=make_pair(1,0);
            for(int k=0;k<mid;k++,w=w*wn)
            {
                pair<double,double> x=a[j+k],y=w*a[j+mid+k];
                a[j+k]=x+y;
                a[j+mid+k]=x-y;
            }
        }
    }
}
void work()
{
    fft(a,1);
    fft(b,1);
    for(int i=0;i<=limit;i++)
        a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=n+m;i++)
        cout<<(int)(a[i].first/limit+0.5)<<' ';
}
void prepare()
{
    while(limit<=n+m)
    {
        limit<<=1;
        l++;
    }
    for(int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
int main()
{
    ios::sync_with_stdio(false);
    cin>>n>>m;
    for(int i=0;i<=n;i++)
        cin>>a[i].first;
    for(int i=0;i<=m;i++)
        cin>>b[i].first;
    prepare();
    work();
    return 0;
}

NTT

背景

因为FFT的精度损失会影响答案,

导致笔者心态@#!@#^&$

相信这也是大多数人的想法

所以NTT横空出世

模仿

我们思考

芝士:FFT

原文:https://www.cnblogs.com/loney-s/p/12063408.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!