这个 FFT ，看得我都 FFT 了

ACM算法日常

发布于 2021-09-07 15:00:24

1.5K0

文章被收录于专栏：ACM算法日常ACM算法日常

FFT 即快速傅立叶变换。在很多计算机领域都用用处，例如数字图像处理、计算机网络。但他在算法竞赛中主要是用于多项式和生成函数相关的题目。

多项式

表达方式

简介

系数表达式，即

f(x)=a_1x^{n-1}+a_2x^{n-2}+a_3x^{n-3}+\cdots +a_{n}

。

坐标形式。每一个坐标用

(x_i,f(x_i))

表示。显然，为了能够表示一个确定的多项式，需要

个不同的坐标来表示。

比较

对于系数表示，多项式加法的时间复杂度是

O(n)

，多项式乘法的时间复杂度是

O(n^2)

。

对于点值表示，多项式加法的时间复杂度同样是

O(n)

，但是乘法的时间复杂度就是

O(2n)

（因为多项式乘法以后最高项次数为

，我们只需要

个坐标表示）。

思考

这样一来，我们就有一个想法，多项式乘法，是不是可以利用坐标表示做多项式乘法特别快这点来优化算法。

于是需要解决的最大的问题就是，多项式两种表示方法之间的互相转换。

求值朴素做法的时间复杂度是

O(n^2)

，即随便选取一个值带入，暴力计算。

差值朴素的做法时间复杂度是

O(n^3)

，即根据坐标进行高斯消元。

在介绍 FFT 之前，得先学习 DFT （离散傅里叶变换）算法。

DFT

由于对一个多项式的点值表达式的取是任意的，所以好的取法可能会使一个算法产生本质性的蜕变。

我们选取

次单位复根作为

来取值。

单位复根

z^n=1(n=1,2,3,\cdots)

，这个方程的复数根

为

次单位根。

单位的

个单位根分别为

x_k=e^{\frac{2\pi ki}{n}}(k=0,1,2,\cdots ,n-1)

。

个单位根在复平面的坐标表示为

(cos(\frac{2\pi k}{n}),sin(\frac{2\pi k}{n}))

，我们将这个记为

w_{n}^{k}

。在复平面上形象的表示的话，就是下图：

单位根在多项式的应用

我们将

个单位根带入多项式可以得到

个因变量结果，记为

y_k

。

我们将

个单位根的倒数

w_{n}^{-k}

作为自变量带入由

y_k

作为系数的多项式，可以得到

z_k=\sum _{i=0}^{n-1}y_i(w_{n}^{k-1})^i\\=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1} a_j(w_{n}^{i})^j)(w_{n}^{-k})^i\\=\sum _{j=0}^{n-1}a_j(\sum _{i=0}^{n-1}(w_n^{j-k})^i)

。

而

\sum _{i=0}^{n-1}(w_n^{j-k})^i

当

j-k=0

的时候，它为

，其余时候，它为

（通过等比数列求和）。于是有

a_i=\frac{z_i}{n}

。

于是可以发现，将

个单位根的倒数带入变换后的多项式，可以得到原来的多项式。

这样一来，我们利用

个单位根的性质，完成了多项式两种表示方式之间的转换。

DFT算法

有了

x_k

的取值，我们就可以得到

y_k

的取值了。

y_k=f(x_k)=\sum _{i=0}^{n-1}a_ix_k

。

直接暴力计算，两个方向转换的时间复杂均为

O(n^2)

。

FFT

那么 FFT 算法是如何优化计算这一过程的？利用分治。

我们把一个多项式的计算分为偶数项的计算和奇数项的计算：

f(x)=xf_{odd}(x^2)+f_{even}(x^2)

f_{odd}(x)=a_1+a_3x+\cdots +a_{n-1}x^{\frac{n}{2}-1}

f_{even}(x)=a_0+a_2x+\cdots +a_{n-2}x^{\frac{n}{2}-1}

也就是说， FFT 的思想就是不断得把一个多项式拆分成奇数多项式和偶数多项式，然后合并两个多项式的信息。

也就是说，如果我们已经得到了

f_{odd}

和

f_{even}

，我们只需要

O(n)

就可以得到

f(x)

了。

每次都能把多项式的长度减小一半，于是时间复杂度就是

O(nlogn)

。

模版

C++ 是自带了复数 stl 的，即：

#include <complex>

但是常数大，跑得慢，不如自己重载的好。

下面的模版必须要保证

是

的整数次幂。

typedef double LD;
const LD PI = acos(-1);
struct C {
    LD r, i;
    C(LD r = 0, LD i = 0): r(r), i(i) {}
};
C operator + (const C& a, const C& b) {
    return C(a.r + b.r, a.i + b.i);
}
C operator - (const C& a, const C& b) {
    return C(a.r - b.r, a.i - b.i);
}
C operator * (const C& a, const C& b) {
    return C(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}

void FFT(C x[], int n, int p) {
    for (int i = 0, t = 0; i < n; ++i) {
        if (i > t) swap(x[i], x[t]);
        for (int j = n >> 1; (t ^= j) < j; j >>= 1);
    }
    for (int h = 2; h <= n; h <<= 1) {
        C wn(cos(p * 2 * PI / h), sin(p * 2 * PI / h));
        for (int i = 0; i < n; i += h) {
            C w(1, 0), u;
            for (int j = i, k = h >> 1; j < i + k; ++j) {
                u = x[j + k] * w;
                x[j + k] = x[j] - u;
                x[j] = x[j] + u;
                w = w * wn;
            }
        }
    }
    if (p == -1)
        FOR (i, 0, n)
            x[i].r /= n;
}

void conv(C a[], C b[], int n) {
    FFT(a, n, 1);
    FFT(b, n, 1);
    FOR (i, 0, n)
        a[i] = a[i] * b[i];
    FFT(a, n, -1);
}

例题

**A * B II**

https://acm.ecnu.edu.cn/problem/3007/

大整数相乘。

把每一位都看成是多项式其中一项的系数，那么大数最后的结果也就是多项式乘法系数的结果。

要进位处理。

Hnoi2017 礼物

显然是要计算

\sum_{i=1}^{n}(x_i-y_{i+x}+y)^2

的最小值，其中$0≤x

展开这个式子，

\sum_{i=1}^{n}(x_i-y_{i+x}+y)^2

=\sum_{i=1}^{n}(x_i^2-2x_iy_{i+x}+2yx_i+y_{i+x}^2-2yy_{i+x}+y^2)

=\sum_{i=1}^{n}(x_i^2+y_{i}^2)-2*\sum_{i=1}^{n}x_iy_{i+x}+2y*\sum_{i=1}^{n}(x_i-y_i)+n*y^2

除了

2*\sum_{i=1}^{n}x_iy_{i+x}

，其他的和

x_i

与

y_i

相关的项都可以在

O(n)

的时间内算出了

那么

2y*\sum_{i=1}^{n}(x_i-y_i)+y^2

配个方，就可以求出最小值了，而

\sum_{i=1}^{n}(x_i^2-y_{i}^2)

是固定的

现在的问题就是求

2*\sum_{i=1}^{n}x_iy_{i+x}

，我们可以用FFT来解决

如果我们把多项式

倒置，我们就能发现

2*\sum_{i=1}^{n}x_iy_{i+x}

式子的

和

的下标和可以相同，我们可以利用多项式乘法同时算出卷积。

设

x'_i=x_{n-i+1}

，那么

2*\sum_{i=1}^{n}x_iy_{i+x}=2*\sum_{i=1}^{n}x'_{n-i+1}y_{i+x}

，这样就可以用FFT算出来了

总的时间复杂度

O(nlogn)

#include<bits/stdc++.h>
#define inf 0x3fffffff
using namespace std;

typedef double LD;
const LD PI=acos(-1);
struct C
{
    LD r,i;
    C(LD r=0,LD i=0):r(r),i(i){}
};
C operator + (const C& a, const C& b){
    return C(a.r+b.r,a.i+b.i);
}
C operator - (const C& a, const C& b){
    return C(a.r-b.r,a.i-b.i);
}
C operator * (const C& a, const C& b){
    return C(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
}
 
void FFT(C x[],int n,int p)
{
    for (int i=0,t=0;i<n;++i)
    {
        if (i>t) swap(x[i],x[t]);
        for (int j=n>>1;(t^=j)<j;j>>=1);
    }
    for (int h=2;h<=n;h<<=1)
    {
        C wn(cos(p*2*PI/h),sin(p*2*PI/h));
        for (int i=0;i<n;i+=h)
        {
            C w(1,0),u;
            for (int j=i,k=h>>1;j<i+k;++j)
            {
                u=x[j+k]*w;
                x[j+k]=x[j]-u;
                x[j]=x[j]+u;
                w=w*wn;
            }
        }
    }
    if (p==-1)
        for (int i=0;i<n;i++) x[i].r/=n;
}
 
void conv(C a[], C b[], int n) {
    FFT(a,n,1);
    FFT(b,n,1);
    for (int i=0;i<n;i++)
        a[i]=a[i]*b[i];
    FFT(a,n,-1);
}

int n,m;
int a[400010],b[400010];
C c[400010],d[400010];

int main()
{
    scanf("%d%d",&n,&m);
    int ans=0;int y=0,x;
    for (int i=1;i<=n;i++)
        scanf("%d",&a[i]),ans+=a[i]*a[i],y+=a[i];
    for (int i=1;i<=n;i++)
        scanf("%d",&b[i]),ans+=b[i]*b[i],y-=b[i];
    x=round(((double)-y)/((double)n));
    ans+=n*x*x+2*x*y;
    int len=1;
    while(len<2*n) len*=2;
    for (int i=1;i<=n;i++)
        c[i-1]=C(a[n-i+1],0);
    for (int i=n;i<=len;i++)
        c[i]=C(0,0);
    for (int i=1;i<=n;i++)
        d[i-1]=C(b[i],0);
    for (int i=1;i<=n;i++)
        d[n+i-1]=C(b[i],0);
    for (int i=2*n;i<=len;i++)
        d[i]=C(0,0);
    conv(c,d,len);
    int Max=0;
    for (int i=n-1;i<2*n-1;i++)
        Max=max(Max,(int)round(c[i].r));
    printf("%d\n",ans-2*Max);
    return 0;
}

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-08-29，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自 ACM算法日常微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度