【学习笔记】开发福利院保护（FFT）

2022/1/23 23:06:30

编程Tag： int 笔记 define fu Omega fft operatorname DFT 福利院

本文主要是介绍【学习笔记】开发福利院保护（FFT），对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

概述

FFT，即 快速傅里叶变换 ，是将多项式乘法从 \(O(n^2)\) 优化到 \(O(n\log n)\) 的算法。

本质上是优化卷积，卷积的一般形式：

\[C(i)=\sum\limits_{i\oplus j=k}A(i)B(i) \]

其中多项式乘法为加法卷积，即：

\[C(i)=\sum\limits_{i+j=k}A(i)B(i) \]

系数表示法：

我们可以用每个项的系数来表示这个多项式：

\[f(x)=a_nx^n+\dots +a_2x^2+a_1x+a_0\Rightarrow f(x) = \{a_0,a_1,\dots,a_n\} \]

点值表示法：

我们可以把多项式看成一个函数，那么从平面上取 \(n+1\) 个点，则可以确定一个 \(n+1\) 项的函数。

原因可以理解成 \(n+1\) 个多项式，也就是 \(n+1\) 个方程，可以解出 \(n+1\) 元的方程。

FFT 的原理就是将对于每一个多项式 \(A(x)\) 和 \(B(x)\) 由系数表示法变成点值表示法，\(A(x)\) 和 \(B(x)\) 相乘得 \(C(x)\) ，再变成系数表示法。将系数表示法变成点值表示法的过程称为 DFT，点值表示法变成系数表示法的过程称为 IDFT。

但如果暴力算的话，\(A(x)\) 和 \(B(x)\) 变成点值表示法是 \(O(n^2)\) 的，再用高斯消元求解是 \(O(n^3)\) ，所以 FFT 可以加速这个过程。

单位根

对于点值表示法的 \(n\) 个 \(A(x)\) 和 \(B(x)\) ，需要快速算出它们的值。

直接算 \(O(n^2)\) 的，发现对于 \(x^n=1\) 的 \(x\) ，可以 \(O(n)\) 计算出 \(A(x)\) 和 \(B(x)\) 。

在实数域中，只有 \(1\) 和 \(-1\) 满足 \(O(n)\) 算出多项式。

在复数域中，则有 \(i\) 和 \(-i\) 满足要求。

考虑复数乘法，表示为 模长相乘，辐角相加 ，则两个模长为 \(1\) 的向量相乘，得到的还是模长为 \(1\) 的向量。

可以定义 \(x^n=1\) 在复数意义下的解为 \(n\) 次复根，即 \(\omega_{n}\) 。这样的复根有 \(n\) 个，表示为 \(\omega_{n}^k\) ，其中 \(k=0,1,\dots,n-1\)

OI-WIKI上的图：

性质：

\[\omega_{n}^k=\omega_{2n}^{2k} \\ \omega_{2n}^{n+k}=-\omega_{2n}^k \]

快速傅里叶变换：

FFT 的基本思想是分治。

DFT

我们将多项式分为奇次项和偶次项处理。

对于一个 \(8\) 项多项式，按照次数分为两组：

\[f(x)=(a_0+a_2x^2+a_4x^4+a_6x^6)+(a_1x+a_3x^3+a_5x^5+a_7x^7)\\=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6) \]

\[G(x)=a_0+a_2x+a_4x^2+a_6x^3\\ H(x)=a_1+a_3x+a_5x^2+a_7x^3 \]

则有：

\[f(x)=G(x^2)+x\times H(x) \]

\[\begin{aligned}\operatorname{DFT}(f(\omega_n^k))&=\operatorname{DFT}(G((\omega_n^k)^2))+\omega_n^k\times \operatorname{DFT}(H((\omega_n^k)^2))\\ &=\operatorname{DFT}(G(\omega_n^{2k})) + \omega_n^k\times \operatorname{DFT}(H(\omega_n^{2k}))\\ &=\operatorname{DFT}(G(\omega_{n/2}^k)) + \omega_n^k \times \operatorname{DFT}(H(\omega_{n/2}^k))\end{aligned} \]

同理可得：

\[\begin{aligned}\operatorname{DFT}(f(\omega_n^{k+n/2}))&=\operatorname{DFT}(G(\omega_n^{2k+n}))+\omega_n^{k+n/2}\times \operatorname{DFT}(H(\omega_n^{2k+n}))\\&=\operatorname{DFT}(G(\omega_n^{2k}))-\omega_n^k\times \operatorname{DFT}(H(\omega_n^{2k}))\\&=\operatorname{DFT}(G(\omega_{n/2}^k))-\omega_n^k\times \operatorname{DFT}(H(\omega_{n/2}^k))\end{aligned} \]

所以我们可以对 \(G\) 和 \(H\) 分别递归求解，得出 \(\operatorname{DFT}(f(\omega_{n}^{k}))\) 和 \(\operatorname{DFT}(f(\omega_{n}^{k+n/2}))\)

位逆序置换

DFT 很明显可以用递归求解，但它还可以继续优化。

考虑到递归的过程是不断对若干个长为 \(2^m\) 的子段拆分成两个长为 \(2^{m-1}\) 的子段，并且每次递归 \(m-1\) 直到 \(m=1\) ，那么可以考虑从 \(m=1\) 的状态开始合并。

容易发现每个数在拆分到 \(m=1\) 时的位置和原始位置的二进制是翻转的，那么求出 rev[i] 表示 \(i\) 的二进制反转后的值即可。

这里给出代码：

Code-rev

  for (int i = 0; i < len; ++i)
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (num - 1));

IDFT

IDFT 相当于已知 \(y_i=f(\omega_{n}^i)\) ，\(i\in\{0,1,\dots,n-\}\) ，求 \(\{a_0,a_1,\dots,a_{n-1}\}\)

我们取单位根的倒数，跑一边 FFT，再将求得的 \(y\) 除以 \(n\) ，得到原来的 \(a\) 。

证明略。

Code-FFT

void FFT(Fu *y, int on) {
  for (int i = 0; i < len; ++i)
    if (i < rev[i]) swap(y[i], y[rev[i]]);
  for (int h = 2; h <= len; h <<= 1) {
    Fu wn = (Fu) {cos(2 * PI / h), sin(on * 2 * PI / h)};
    for (int j = 0; j < len; j += h) {
      Fu w = (Fu) {1, 0};
      for (int k = j; k < j + h / 2; ++k) {
				Fu 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 (int i = 0; i < len; ++i)
      y[i].x /= len;
}

例题

P3803

模板题，给出代码。

Code

#include <bits/stdc++.h>

#define LL long long
#define ULL unsigned long long
#define DB double
#define PR pair <int, int>
#define MK make_pair
#define pb push_back
#define fi first
#define se second
#define RI register int
#define Low(x) (x & (-x))

using namespace std;

const int kN = 4e6 + 5;
const DB PI = acos(-1.0);

int n, m, len = 1, num = 0, rev[kN];

struct Fu {
  DB x, y;

  Fu operator + (const Fu &K) const {return (Fu) {x + K.x, y + K.y};}
  Fu operator - (const Fu &K) const {return (Fu) {x - K.x, y - K.y};}
  Fu operator * (const Fu &K) const {return (Fu) {x * K.x - y * K.y, x * K.y + y * K.x};}
} f[kN], g[kN];

void FFT(Fu *y, int on) {
  for (int i = 0; i < len; ++i)
    if (i < rev[i]) swap(y[i], y[rev[i]]);
  for (int h = 2; h <= len; h <<= 1) {
    Fu wn = (Fu) {cos(2 * PI / h), sin(on * 2 * PI / h)};
    for (int j = 0; j < len; j += h) {
      Fu w = (Fu) {1, 0};
      for (int k = j; k < j + h / 2; ++k) {
				Fu 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 (int i = 0; i < len; ++i)
      y[i].x /= len;
}

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", &g[i].x);
  while (len <= n + m) len <<= 1, num++;
  for (int i = 0; i < len; ++i)
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (num - 1));
  FFT(f, 1), FFT(g, 1);
  for (int i = 0; i < len; ++i) f[i] = f[i] * g[i];
  FFT(f, -1);
  for (int i = 0; i <= n + m; ++i) printf("%d ", (int) (f[i].x + 0.5));
  return 0;
}

这篇关于【学习笔记】开发福利院保护（FFT）的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！