从N方到NlogN的转变——FFT

1.Why

为什么我们信息学竞赛要用到FFT
因为我们要优化卷积啊
将n边为log是一个非常大的优化啊

What

什么是FFT呢
wiki上是这么说的

快速傅里叶变换(英语:Fast Fourier Transform, FFT),是计算序列的离散傅里叶变换(DFT)或其逆变换的一种算法。傅里叶分析将信号从原始域(通常是时间或空间)转换到频域的表示或者逆过来转换。FFT会通过把DFT矩阵分解为稀疏(大多为零)因子之积来快速计算此类变换。[1] 因此,它能够将计算DFT的复杂度从只用DFT定义计算需要的 $O(n^{2})$,降低到 $O(n\log n)$,其中 n 为数据大小。

啥你说啥我听不懂(黑人问号?)

如果简单的说就是求两个多项式的乘积

How

其实关于FFT的具体理论很复杂,我一篇博客肯定不能接受清楚
所以请有兴趣的人移步Google 或 baidu

如果简单来说FFT用的分制的方法
所以才能有log啊

将一个多项式按奇偶分开
但如果我们直接找值去计算时间复杂度还是$O(n^2)$的

所以出现了一个特殊的东西他叫做n次单位根

然后利用他的性质就可以搞了

先贴两个板子

UOJ 34 多项式乘法

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int maxn = 1 << 18 | 5;
namespace FFT
{
const double pi = acos(-1.0);
struct Complex
{
    double x, y;
    Complex operator+(const Complex &a)
    {
        return (Complex){x + a.x, y + a.y};
    }
    Complex operator-(const Complex &a)
    {
        return (Complex){x - a.x, y - a.y};
    }
    Complex operator*(const Complex &a)
    {
        return (Complex){x * a.x - y * a.y, x * a.y + y * a.x};
    }
    Complex operator*(const double &a)
    {
        return (Complex){x * a, y * a};
    }
    Complex Get()
    {
        return (Complex){x, -y};
    }
} c[maxn], c1[maxn], V = (Complex){0.5, 0} * (Complex){0, -0.5};
int rev[maxn], N, FU;
double Inv;
void FFt(Complex *a, int op)
{
    Complex w, wn, t;
    for (int i = 0; i < N; i++)
        if (i < rev[i])
            swap(a[i], a[rev[i]]);
    for (int k = 2; k <= N; k <<= 1)
    {
        wn = (Complex){cos(op * 2 * pi / k), sin(op * 2 * pi / k)};
        for (int j = 0; j < N; j += k)
        {
            w = (Complex){1, 0};
            for (int i = 0; i < (k >> 1); i++, w = w * wn)
            {
                t = a[i + j + (k >> 1)] * w;
                a[i + j + (k >> 1)] = a[i + j] - t;
                a[i + j] = a[i + j] + t;
            }
        }
    }
    if (op == -1)
        for (int i = 0; i < N; i++)
            a[i] = a[i] * Inv;
}
void FFt(int *a, int *b, int *res, int n)
{
    int j;
    N = 1;
    while (N < n)
        N <<= 1;
    FU = N - 1;
    Inv = 1. / N;
    for (int i = 0; i < N; i++)
        if (i & 1)
            rev[i] = (rev[i >> 1] >> 1) | (N >> 1);
        else
            rev[i] = (rev[i >> 1] >> 1);
    for (int i = 0; i < N; i++)
        c[i].x = a[i], c[i].y = b[i];
    FFt(c, 1);
    for (int i = 0; i < N; i++)
        j = (N - i) & FU, c1[i] = (c[i] + c[j].Get()) * (c[i] - c[j].Get()) * V;
    FFt(c1, -1);
    for (int i = 0; i < N; i++)
        res[i] = round(c1[i].x);
}
} // FFT

int main(int argc, char const *argv[])
{
    static int a[maxn],b[maxn];
    int n1, n2, m;
    scanf("%d%d", &n1, &n2);
    m = n1 + n2 + 1;
    for (int i = 0; i <= n1; i++)
        scanf("%d", &a[i]);
    for (int i = 0; i <= n2; i++)
        scanf("%d", &b[i]);
    FFT::FFt(a, b, a, m);
    for (int i = 0; i < m; i++)
        printf("%d ", a[i]);
    return 0;
}

COGS 1473. 很强的乘法问题

#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
const double pi = acos(-1.);
struct Complex
{
    double x, y;
    Complex() { ; }
    Complex(double a, double b) : x(a), y(b) {}
    Complex operator+(const Complex &a) { return Complex(x + a.x, y + a.y); }
    Complex operator-(const Complex &a) { return Complex(x - a.x, y - a.y); }
    Complex operator*(const Complex &a) { return Complex(x * a.x - y * a.y, x * a.y + y * a.x); }
    Complex operator*(const double a) { return Complex(x * a, y * a); }
    Complex Get() { return Complex(x, -y); }
} A[150005 << 3], B[150005 << 3];
int rev[150005 << 3], N, FU;
double INV;
void FFt(Complex *a, int op)
{
    Complex w, wn, t;
    for (int i = 0; i < N; i++)
        if (i < rev[i])
            std::swap(a[i], a[rev[i]]);
    for (int k = 2; k <= N; k <<= 1)
    {
        wn = Complex(cos(op * 2 * pi / k), sin(op * 2 * pi / k));
        for (int j = 0; j < N; j += k)
        {
            w = Complex(1, 0);
            for (int i = 0; i < (k >> 1); i++, w = w * wn)
            {
                t = a[i + j + (k >> 1)] * w;
                a[i + j + (k >> 1)] = a[i + j] - t;
                a[i + j] = a[i + j] + t;
            }
        }
    }
    if (op == -1)
        for (int i = 0; i < N; i++)
            a[i] = a[i] * INV;
}
void FFt(const int *a, const int *b, int *res, int n)
{
    N = 1;
    while (N < n)
        N <<= 1;
    INV = 1. / N;
    for (int i = 0; i < N; i++)
        if (i & 1)
            rev[i] = (rev[i >> 1] >> 1) | (N >> 1);
        else
            rev[i] = (rev[i >> 1] >> 1);
    for (int i = 0; i < N; i++)
        A[i].x = a[i], B[i].x = b[i];
    FFt(A, 1), FFt(B, 1);
    for (int i = 0; i < N; i++)
        A[i] = A[i] * B[i];
    FFt(A, -1);
    for (int i = 0; i < N; i++)
        res[i] = round(A[i].x);
}
char s[150005 << 2];
struct BigNum
{
    int len;
    int a[1000005];
    void read()
    {
        scanf("%s", s);
        len = strlen(s);
        for (int i = len - 1, j = 0; i >= 0; i--, j++)
            a[j] = s[i] - '0';
    }
    BigNum operator*(const BigNum &b)
    {
        BigNum ans;
        ans.len = len + b.len - 1;
        FFt(a, b.a, ans.a, ans.len + 1);
        for (int i = 0; i <= ans.len + 2; i++)
        {
            if (ans.a[i] > 9)
            {
                ans.a[i + 1] += ans.a[i] / 10;
                ans.a[i] %= 10;
            }
        }
        while (ans.a[ans.len])
            ans.len++;
        return ans;
    }
    void print()
    {
        for (int i = len - 1; i >= 0; i--)
            printf("%d", a[i]);
        printf("\n");
    }
} a, b;
 
int main()
{
    freopen("bettermul.in","r",stdin);
    freopen("bettermul.out","w",stdout);
    a.read();
    b.read();
    (a * b).print();
    //while (1)
        ;
}
本文作者 : NekoMio
知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
本文链接 : https://www.nekomio.com/2017/08/14/96/
上一篇
下一篇