首页 > 其他分享 >FFT

FFT

时间：2023-02-19 08:22:04浏览次数：29  
标签：frac int FFT len Complex const omega

概念

FFT全称Fast Fourier Transformation， 即快速傅里叶变换，可在 \(O(n \log n)\) 的复杂度计算多项式乘法

一般的多项式乘法是这样的：

\[\begin{aligned} & \ (x^2 + 2x + 1)(3x ^ 2 - x + 7) \\ =& \ x ^ 2(3x ^ 2 - x + 7) + 2x(3x ^ 2 - x + 7) + (3x ^ 2 - x + 7) \\ =& \ 3x ^ 4 - x ^ 3 + 7x ^ 2 + 6x ^ 3 - 2 x ^ 2 + 14x + 3x ^ 2 - x + 7 \\ =& \ 3x ^ 4 + 5x ^ 3 + 8x ^ 2 + 13x + 7 \end{aligned} \]

注意到这样是 \(O(n ^ 2)\) 的，慢的一批

FFT 基本思想

  • 多项式的表达方法
    • 系数表示法：\(F(x)\) 表示一个形如 \(a_{n} x^{n} + a_{n - 1} x^{n - 1} + \cdots + a_{1} x + a_{0}\) 的多项式，其中 \(a\) 是系数，且通常把 \(a_{i}\) 写作 \(F_i\) 。上式也可以写作 \(\sum_{i = 0}^{n} F_i x ^i\)
    • 点值表示法：注意到平面上的 \(n + 1\) 个点可以唯一确定一个 \(n\) 次多项式，所以我们可以用这 \(n + 1\) 个点表示一个多项式
  • 卷积
    • 注意到令 \(C = A \times B\) （其中 \(A, B, C\) 均为多项式），则 \(C_k = \sum_{i = 0}^{k} A_{i} \times B_{k - i} = \sum\limits_{i + j = k} A_{i} \times B_{j}\)
    • 定义形如 \(C_k = \sum\limits_{i \oplus j = k} A_{i} \times B_{j}\) 的式子称为卷积（其中 \(\oplus\) 为运算符号）
    • 显然，多项式乘法就是加法卷积
  • 点值表示法的乘法运算
    • 点值表示法下，多项式的乘法即为这两个多项式对应点的 \(y\) 坐标的乘积
    • 例如我们可以用 \((1, 4), (2, 9), (3, 16)\) 来表示 \(x ^ 2 + 2x + 1\) ，用 \((1, 9), (2,17), (3, 31)\) 表示 \(3x ^ 2 - x + 7\) ，那么这两个多项式的结果 \(3x ^ 4 + 5x ^ 3 + 8x ^ 2 + 13x + 7\) 均过点 \((1, 36), (2, 153), (3, 496)\) 。但是这 \(3\) 个点无法表示结果，我们再加入一些点就行了。
    • 注意到点值表示法下，多项式乘法的时间复杂度是线性的
    • 那么我们可以利用这个性质，将原式用点值表示法表示后相乘，在转化为系数表达法，这就是 FFT 的算法流程。其中系数表示法转点值表示法称为 DFT ，点值表示法转系数表示法称为 IDFT（DFT的逆运算）

单位根及其性质

前置芝士：复数（普通高中教科书·数学（A版）必修 第二册 第七章）

附上复数结构体的代码

struct Complex {
	double x, y;
	
	inline Complex(const double xx = 0, const double yy = 0) {
		x = xx, y = yy;
	}
	
	inline Complex operator + (const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	
	inline Complex operator - (const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	
	inline Complex operator * (const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
	
	inline Complex operator / (const Complex &b) const {
		double t = b.x * b.x + b.y * b.y;
		return Complex((x * b.x + y * b.y) / t, (y * b.x - x * b.y) / t);
	}
};

定义 \(n\) 次单位根 （ \(n \in Z\) ）是方程 \(x^n = 1\) 的复数解

亦可将 \(x\) 表示为 \(r^n(\cos{n \theta} + i\sin{n \theta})\)

不难发现方程的每个根都是单位圆上的 \(n\) 等分点

我们从 \(1\) 开始，依次将这些点表示为 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\)

一些性质：

  • 对于任意的 \(n\) ，都有 \(\omega_n^0 = 1\)
  • \(\omega_n^k = (\omega_n^1)^k\)
  • \(\omega_n^j \times \omega_n^k = \omega_n^{j + k}\)
  • \(\omega_{pn}^{pk} = \omega_n^k\)
  • \(\forall \ 2 | n \ ,\omega_n^{k + \frac{n}{2}} = -\omega_n^k\)
  • \(\omega_n^k = \omega_n^{k \bmod n}\)

DFT的基本思想

下面保证 \(n = 2^k (k \in N^{\ast})\)

若不满足，补上剩下的项即可（补的项系数均为 \(0\) ）

我们先将 \(n - 1\) 次多项式的系数按照奇偶拆分为这样

\[\begin{aligned} F(x) &= F_0 + F_1 x +\cdots + F_n x^n \\ &= (F_0 + F_2 x ^ 2 + \cdots + F_{n - 2}) + (F_1 x + F_3 x ^ 3 + \cdots + F_{n - 1} x^{n - 1}) \end{aligned} \]

令

\[A = F_0 + F_2 x + \cdots + F_{n - 2} x ^{\frac{x}{2} - 1} \\ B = F_1 + F_3 x + \cdots + F_{n - 1} x ^{\frac{x}{2} - 1} \]

则

\[F(x) = A(x ^ 2) + xB(x ^ 2) \]

我们将 \(\omega_n^k , (k < \dfrac{n}{2})\) 带入，得

\[F(\omega_n^k) = A((\omega_n^k) ^ 2) + \omega_n^k B((\omega_n^k) ^ 2) \]

由 \((\omega_n^k)^2 = \omega_n^{2k} = \omega_{\frac{n}{2}}^k\) ，得

\[F(\omega_n^k) = A(\omega_{\frac{n}{2}}^k) + \omega_n^k B(\omega_{\frac{n}{2}}^k) \]

我们将 \(\omega_n^{k + \frac{n}{2}} , (k < \dfrac{n}{2})\) 带入 \(F(x) = A(x ^ 2) + xB(x ^ 2)\) ，得

\[F(\omega_n^{k + \frac{n}{2}}) = A((\omega_n^{k + \frac{n}{2}}) ^ 2) + \omega_n^{k + \frac{n}{2}} B((\omega_n^{k + \frac{n}{2}}) ^ 2) \]

由 \((\omega_n^j)^k = \omega_n^{jk}\)， 得

\[F(\omega_n^{k + \frac{n}{2}}) = A(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}} B(\omega_n^{2k + n}) \]

由 \(\omega_n^k = \omega_n^{k \bmod n}\) ，得

\[F(\omega_n^{k + \frac{n}{2}}) = A(\omega_n^{2k}) + \omega_n^{k + \frac{n}{2}} B(\omega_n^{2k}) \]

由 \(\omega_n^{2k} = \omega_{\frac{n}{2}}^k\) ，得

\[F(\omega_n^{k + \frac{n}{2}}) = A(\omega_{\frac{n}{2}}^k) + \omega_n^{k + \frac{n}{2}} B(\omega_{\frac{n}{2}}^k) \]

由 \(\omega_n^{k + \frac{n}{2}} = -\omega_n^k\) ，得

\[F(\omega_n^{k + \frac{n}{2}}) = A(\omega_{\frac{n}{2}}^k) - \omega_n^kB(\omega_{\frac{n}{2}}^k) \]

对比上一个式子 \(F(\omega_n^k) = A(\omega_{\frac{n}{2}}^k) + \omega_n^k B(\omega_{\frac{n}{2}}^k)\) ，发现他们的区别仅仅只有 \(B\) 的正负

不难发现，若我们已知 \(A, B\) 在 \(\omega_{\frac{n}{2}}^0, \omega_{\frac{n}{2}}^1, \cdots,\omega_{\frac{n}{2}}^{\frac{n}{2} - 1}\) 的 点值表示，就可以在线性时间复杂度得出 \(\omega_{n}^0, \omega_{n}^1, \cdots,\omega_{n}^{n - 1}\) ，即 \(F(x)\) 的点值表示

那么我们就可以不断分治求解DFT

代码：

void DFT(Complex *f, int len) {
    if (len == 1)
        return ;

    Complex *A = f, *B = f + (len >> 1);

    for (int i = 0; i < len; ++i)
        tmp[i] = f[i];

    for (int i = 0; i < (len >> 1); ++i) {
        A[i] = tmp[i << 1];
        B[i] = tmp[i << 1 | 1];
    }

    DFT(A, len >> 1), DFT(B, len >> 1);
    
    Complex tG(cos(2 * Pi / len), sin(2 * Pi / len)), buf(1, 0);

    for (int i = 0; i < (len >> 1); ++i) {
        Complex tt = buf * B[i];
        tmp[i] = A[i] + tt;
        tmp[i + (len >> 1)] = A[i] - tt;
        buf = buf * tG;
    }

    for (int i = 0; i < len; ++i)
        f[i] = tmp[i];
}

IDFT的基本思想

我们令用单位根点值表达的多项式 \(F(x)\) 的系数表达式为 \(G(x)\) ，则有 \(G_k = \sum_{i = 0}^{n - 1} (\omega_n^k)^i F_i\)

结论： \(n F_k = \sum_{i = 0}^{n - 1} (\omega_n^{-k})^iG_i\)

证明：

\[\begin{aligned} 右边 &= \sum_{i = 0}^{n - 1} (\omega_n^{-k})^iG_i \\ &= \sum_{i = 0}^{n - 1} \omega_n^{-ki}(\sum_{j = 0}^{n - 1} \omega_n^{ij} F_j) \\ &= \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} \omega_n^{-ki} \omega_n^{ij} F_j \\ &= \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} \omega_n^{i(j - k)} F_j \end{aligned} \]

当 \(j = k\) 时，原式 \(=\sum_{i = 0}^{n - 1}\omega_n^0F_k = nF_k\)

当 \(j \not = k\) 时，设 \(t = j - k\) ，则原式 \(= \sum_{i = 0}^{n - 1}\omega_n^{ip}F_{k + p} = \omega_n^p(\sum_{i = 0}^{n - 1}\omega_n^i)F_{k + p}\)

注意到 \(\sum_{i = 0}^{n - 1}\omega_n^i\) 实际就是等比数列求和，利用等比数列求和公式得到

\[\sum_{i = 0}^{n - 1}\omega_n^i = \dfrac{\omega_n^0 - \omega_n^n}{1 - \omega_n^1} = \dfrac{\omega_n^0 - \omega_n^0}{1 - \omega_n^1} = 0 \]

所以此时原式的值为 \(0\)

综上所述，\(n F_k = \sum_{i = 0}^{n - 1} (\omega_n^{-k})^iG_i\)

再看这个结论，实际就是把 \(G\) 当作系数，带入 \(\{\omega_n^0, \omega_n^{-1}, \cdots, \omega_n^{-n + 1}\}\) 求值

因为 \(\omega_n^{-n + 1} = \cos{\dfrac{2 \pi}{n}} + i \sin{\dfrac{2 \pi}{n}}\) ，我们直接利用 \((\omega_n^{-1})^k = \omega_n^{-k}\) 求值即可

代码和 DFT区别不大，只要把 tG 的 sin 前面加上负号即可

代码：

#include <bits/stdc++.h>
using namespace std;
const double Pi = acos(-1);
const int N = 2e6 + 7;

struct Complex {
	double x, y;
	
	inline Complex(const double xx = 0, const double yy = 0) {
		x = xx, y = yy;
	}
	
	inline Complex operator + (const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	
	inline Complex operator - (const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	
	inline Complex operator * (const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
	
	inline Complex operator / (const Complex &b) const {
		double t = b.x * b.x + b.y * b.y;
		return Complex((x * b.x + y * b.y) / t, (y * b.x - x * b.y) / t);
	}
} f[N << 1], p[N << 1], tmp[N << 1];

int n, m;

void FFT(Complex *f, int len, int sign) {
    if (len == 1)
        return ;

    Complex *A = f, *B = f + (len >> 1);

    for (int i = 0; i < len; ++i)
        tmp[i] = f[i];

    for (int i = 0; i < (len >> 1); ++i) {
        A[i] = tmp[i << 1];
        B[i] = tmp[i << 1 | 1];
    }

    FFT(A, len >> 1, sign), FFT(B, len >> 1, sign);
    
    Complex tG(cos(2 * Pi / len), sin(2 * Pi / len) * sign), buf(1, 0);

    for (int i = 0; i < (len >> 1); ++i) {
        Complex tt = buf * B[i];
        tmp[i] = A[i] + tt;
        tmp[i + (len >> 1)] = A[i] - tt;
        buf = buf * tG;
    }

    for (int i = 0; i < len; ++i)
        f[i] = tmp[i];
}

signed main() {
	scanf("%d%d", &n, &m);
	
	for (int i = 0; i <= n; ++i)
		scanf("%lf", &f[i].x);
	
	for (int i = 0; i <= m; ++i)
		scanf("%lf", &p[i].x);
	
	for (m += n, n = 1; n <= m; n <<= 1);
	
	FFT(f, n, 1), FFT(p, n, 1);
	
	for (int i = 0; i < n; ++i)
		f[i] = f[i] * p[i];
	
	FFT(f, n, -1);
	
	for (int i = 0; i <= m; ++i)
		printf("%d ", (int) (f[i].x / n + 0.5));
	
	return 0;
}

优化

可以发现这种写法常数非常大

我们尝试观察一下这个递归分组的过程

0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7

注意到递归后 \(6\) 排在了 \(3\) 的位置，且 \(6 = (110)_2, 3 = (011)_2\)

观察其他的数字后，我们发现递归的本质就是二进制翻转

怎么求呢？我们递推实现

我们设 \(i\) 二进制翻转后的数为 \(tr_i\)

假设我们现在要求 \(i\) 二进制翻转后的数，注意到 \(i = (\cdots x)_2 (x = i \bmod 2)\)

那么有 \(tr_i = (x\cdots)_2\) ，其中 \(\cdots\) 的部分我们直接用之前求出的 \(tr\) 值求解即可

代码：

for (int i = 0; i < n; ++i)
	tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);

为了减小常数，我们甚至可以化递归为迭代，也就是这样

#include <bits/stdc++.h>
using namespace std;
const double Pi = acos(-1);
const int N = 5e6 + 7;

struct Complex {
	double x, y;
	
	inline Complex(const double xx = 0, const double yy = 0) {
		x = xx, y = yy;
	}
	
	inline Complex operator + (const Complex &b) const {
		return Complex(x + b.x, y + b.y);
	}
	
	inline Complex operator - (const Complex &b) const {
		return Complex(x - b.x, y - b.y);
	}
	
	inline Complex operator * (const Complex &b) const {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
	
	inline Complex operator / (const Complex &b) const {
		double t = b.x * b.x + b.y * b.y;
		return Complex((x * b.x + y * b.y) / t, (y * b.x - x * b.y) / t);
	}
} f[N], p[N];

int tr[N];

int n, m, limit, highest;

inline void FFT(Complex *f, int sign) {
	for (int i = 0; i < limit; ++i)
		if (i < tr[i])
			swap(f[i], f[tr[i]]);
	
	for (int p = 2; p <= limit; p <<= 1) {
		int len = p >> 1;
		Complex tG(cos(2 * Pi / p), sin(2 * Pi / p) * sign);
		
		for (int k = 0; k < limit; k += p) {
			Complex buf(1, 0);
			
			for (int l = k; l < k + len; ++l) {
				Complex tt = buf * f[len + l];
				f[len + l] = f[l] - tt;
				f[l] = f[l] + tt;
				buf = buf * tG;
			}
		}
	}
}

signed main() {
	scanf("%d%d", &n, &m);
	
	for (int i = 0; i <= n; ++i)
		scanf("%lf", &f[i].x);
	
	for (int i = 0; i <= m; ++i)
		scanf("%lf", &p[i].x);
	
	limit = 1, highest = n + m;
	
	while (limit <= highest)
		limit <<= 1;
	
	for (int i = 0; i < limit; ++i)
		tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? (limit >> 1) : 0);
	
	FFT(f, 1), FFT(p, 1);
	
	for (int i = 0; i < limit; ++i)
		f[i] = f[i] * p[i];
	
	FFT(f, -1);
	
	for (int i = 0; i <= highest; ++i)
		printf("%d ", (int) (f[i].x / limit + 0.5));
	
	return 0;
}

标签：frac,int,FFT,len,Complex,const,omega
From： https://www.cnblogs.com/wshcl/p/17134184.html

相关文章

  • 用于ARM上的FFT与IFFT源代码-C语言
    /*********************************************************************************程序名称：快速傅里叶变换(FFT)**程序描述：本程序实现快速傅里叶变换**程序作者：宋......
  • 深入理解 FFT
    理论前置知道啥是多项式（即\(f(x)=\displaystyle\sum_{i=0}^{n-1}f_ix^i\)这一类东西）。知道啥是多项式的卷积（即\((f\timesg)(x)=h(x)\)，其中\(h_i=\displaystyle\sum_......
  • FFT&NTT
    FFT快速傅里叶变换<NTTFFT和NTT是\(O(nlogn)\)处理两个多项式相乘的算法（FFT<NTT）前置知识复数一个复数可以表示为\[z=a+ib~~a,b\inR\]我们把他看做平面上的一个点，......
  • 算法学习笔记(17): 快速傅里叶变换（FFT）
    快速傅里叶变换（FFT）有趣啊，都已经到NOI的难度了，救命首先，我们先讲述一下前置知识。已经明白的读者请移步后文虚数定义：\(z=a+bi\)，其中\(a,b\inR\\i=\sqrt{-1......
  • 基2和基4FFT
    1.FFT的必要索引变换基2算法需要位顺序的反转位逆序，而基4算法需要首先构成一个2位的数字，再反转这些数字，称为数字逆序。1.1位逆序和数字逆序2.FFT的复数乘法转实数乘法......
  • FFT快速傅里叶变换
    FFT快速傅里叶变换DFT：离散傅里叶变换—>\(O(n^2)\)计算多项式乘法FFT：快速傅里叶变换—>\(O(n\logn)\)计算多项式乘法FNTT/NTT：快速傅里叶变换的优化版—>优化常数及误差......
  • 快速傅里叶逆变换（IFFT）
    本文作者为JustinRochester。目录地址上一篇下一篇快速傅里叶逆变换（IFFT）......
  • 快速傅里叶变换（FFT）的分治实现
    本文作者为JustinRochester。目录地址上一篇下一篇......
  • 手撕fft系列之频移fftshift源码解析
    壹：fft在数字信号处理领域是一个神一样的存在。要好好熟悉一下。这里给出频移的算法源码解析。所谓的频移，就是把数字信号的频频顺序打乱，移动一些。这个在防止啸叫和......
  • 傅里叶级数_傅里叶变换_离散傅里叶变换(DFT)_快速傅里叶变换(FFT)
    一、傅里叶级数 核心思想：周期函数\(f(t)\)可以看成是一系列频率(周期)不同的周期函数\({f_k}(t)\)的叠加，即：\[\begin{array}{c}f(t)={c_1}{f_1}(t)+{c_2}{f_2}......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99