多项式

一个\(n\)次多项式可以表示为\(A(x)=\sum_{i=0}^{n}a_i x^i\)，另一个\(n\)次多项式可以表示为\(B(x)=\sum_{i=0}^{n}b_i x^i\)。

• 多项式加法
  将\(A(x)\)和\(B(x)\)相加，得到多项式\(C(x)=\sum_{i=0}^{n} (a_i+b_i) x^i\)。
  复杂度是\(O(n)\)的。
• 多项式乘法
  将\(A(x)\)和\(B(x)\)相乘，得到多项式\(C(x)=\sum_{i=0}^{n}\sum_{j=0}^{n}a_ib_jx^{i+j}\)。
  复杂度是\(O(n^2)\)的，不是很优，于是我们决定想办法来解决它。

快速傅里叶变换 (FFT)

概述

快速傅里叶变换，英文名\(Fast\ Fourier\ Transform\)，它可以在\(O(nlog(n))\)的复杂度内完成多项式乘法。

前置技能

函数的表达方法

• 系数表达法
  通过\(a_0,a_1,...,a_n\)来表示函数\(f(x)=\sum_{i=0}^{n}a_ix^i\)
• 点值表达法
  选取平面上\(n+1\)个不同的点来表示函数\(f(x)=\{ (x_0,y_0), (x_1,y_1), ..., (x_n, y_n)\}\)

我们发现，通过点值表达法，我们可以在\(O(n)\)的复杂度内计算出\(C(x)\)，于是我们考虑将系数和点值结合起来。

单位复数根

在复平面上，单位圆有一些其妙的性质：
我们默认\(n=2^t,t\in N\)

• \(w_n^k=w_{2n}^{2k}\)
• \(w_n^{k+\frac{n}{2}}=-w_n^k\)
• \(w_n^k=w_n^{k+n}\)
• 对于任意\(k\)，均满足\(\sum_{j=0}^{n-1}(w_n^k)^j=0\)

推导

默认\(n=2^t,t\in N\)，不足的位置上系数置为\(0\)，令函数\(f(x)=a_0+a_1x+a_2x^2+...+a_nx^n\)
那么\(f(x)=(a_0+a_2x^2+a_4x^4+...+a_nx^n)+(a_1+a_3x^2+a_5x^4+...+a_{n-1}x^{n-2})x\)
令\(A(x)=a_0+a_2x+a_4x^2+...+a_nx^{\frac{n}{2}}\)，\(B(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n-2 }{2}}\)，
则\(f(x)=A(x^2)+B(x^2)x\)。
代入\(x=w_n^k\)，有\(f(w_n^k)=A(w_n^{2k})+B(w_n^{2k})w_n^k=A(w_{\frac{n}{2}}^{k})+B(w_{\frac{n}{2}}^{k})w_n^k\)；
代入\(x=w_n^{k+\frac{n}{2}}\)，有\(f(w_n^{k+\frac{n}{2}})=A(w_n^{2k+n})+B(w_n^{2k+n})w_n^{k+\frac{n}{2}}=A(w_{\frac{n}{2}}^{k})-B(w_{\frac{n}{2}}^{k})w_n^k\)。
我们把两式放在一起：
\(\begin{cases} f(w_n^k)=A(w_{\frac{n}{2}}^{k})+B(w_{\frac{n}{2}}^{k})w_n^k \\ f(w_n^{k+\frac{n}{2}})=A(w_{\frac{n}{2}}^{k})-B(w_{\frac{n}{2}}^{k})w_n^k \end{cases}\)
发现刚好一个加号一个减号，并且下标由\(n\)变为\(\frac{n}{2}\)，因此直接递归下去即可。
这里直接给出结论：
\(FFT\) 表示：\(FFT(k)=\sum_{n=0}^{N-1}IFFT(n) w_N^{kn}\)；
\(IFFT\) 表示：\(IFFT(k)=\frac{1}{N}\sum_{n=0}^{N-1} FFT(n) w_N^{-kn}\)。
可以发现\(IFFT\)和\(FFT\)相比，只在上指标加了个负号，并且乘上了\(\frac{1}{N}\)。
时间复杂度 \(O(nlog(n))\)。

蝴蝶变换

如果每次处理的时候我们都传一个数组，很容易就 \(RE\) 或 \(MLE\) ，我们想办法尽可能避免递归。
神奇的蝴蝶变换出场了。
我们先观察一下\(n=8\)的情况：

观察可知，变化前后，\(k\)的二进制位刚好做了一次反转操作。
于是，我们可以进行如下方式，预处理出\(r[]\)数组：
for (int i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << L - 1);
好了，这就是所谓的蝴蝶变换。

例题

题目链接：洛谷P3803 多项式乘法

Code

const int N = 4000001;
const double PI = acos(-1.0);
struct Comp {
  double x, y;
  Comp (double _x = 0, double _y = 0) {
    x = _x, y = _y;
  }
} a[N], b[N];
Comp operator + (Comp a, Comp b) {
  return Comp(a.x + b.x, a.y + b.y);
}
Comp operator - (Comp a, Comp b) {
  return Comp(a.x - b.x, a.y - b.y);
} 
Comp operator * (Comp a, Comp b) {
  return Comp(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
Comp operator / (Comp a, double b) {
  return Comp(a.x / b, a.y / b);
}

int n, m;
namespace Poly {
  int r[N], lim = 1, L = 0;
  void getR(int n) {
    lim = 1, L = 0;
    while (lim <= n) lim <<= 1, L++;
    for (rint i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << L - 1);
  }
  void fft(Comp *a, int opt) {
    for (rint i = 0; i < lim; i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (rint mid = 1; mid < lim; mid <<= 1) {
      Comp Wn(cos(PI / mid), opt * sin(PI / mid));
      for (rint R = mid << 1, j = 0; j < lim; j += R) {
        Comp w(1, 0);
        for (rint k = 0; k < mid; k++, w = w * Wn) {
          Comp x = a[j + k], y = w * a[j + k + mid];
          a[j + k] = x + y;
          a[j + k + mid] = x - y;
        }
      }
    }
    if (opt == -1) {
      for (rint i = 0; i < lim; i++) a[i] = a[i] / lim;
    }
  }
  void mul(Comp *a, Comp *b, int n, int m) {
    getR(n + m);
    fft(a, 1), fft(b, 1);
    for (rint i = 0; i < lim; i++) a[i] = a[i] * b[i];
    fft(a, -1);
  }
}

int main() {
  n = read(), m = read();
  rep(i, 0, n) a[i].x = read();
  rep(i, 0, m) b[i].x = read();
  Poly::mul(a, b, n, m);
  rep(i, 0, n + m) print(int(a[i].x + 0.5)), putchar(‘ ‘);
}

快速数论变换（NTT）

概述

快速数论变换，英文名\(Number\ Theoretic\ Transform\)，能解决在特殊模数的模意义下的多项式乘法。

前置技能

• 阶
  若\(a\)与\(p\)互素，则满足\(a^n\mod 1(mod\ p)\)最小的\(n\)，称之为\(a\)模\(p\)的阶，记作\(δ_p(a)\)。
• 原根
  设\(p\)为正整数，\(a\)是整数，则若\(δ_p(a)=\varphi(p)\)，则称\(a\)为模\(p\)的一个原根。
  由此可见原根的个数是不唯一的。若模数\(p\)有原根，那么它一定有\(\varphi(\varphi(p))\)个原根。
• 原根的性质

1. FFT中单位根的4个性质，原根均满足；
2. 若\(p\)为素数，设\(g\)为\(p\)的原根，那么\(g^i \mod \ p (0<i<p)\)的结果恰好构成一个全排列\(\{ 1,2,3,...,p-1 \}\)；

推导

与\(FFT\)类似，有一个结论：\(g^{\frac{p-1}{n}}\)与\(Wn\)在模\(p\)意义下同余。
所以我们利用原根来代替单位根，实现整数层面的多项式乘法。

质数原根表

例题

题目链接：洛谷P3803 多项式乘法
没错，还是这题233

Code

const int N = 4000001;
const int mod = 998244353;
const int G = 3;
const int Gi = 332748118;
int a[N], b[N];
int n, m;
int qpow(int a, int b) {
  int ret = 1;
  while (b > 0) {
    if (b & 1) ret = 1ll * ret * a % mod;
    a = 1ll * a * a % mod;
    b >>= 1;
  }
  return ret;
}
namespace Poly {
  int r[N], lim = 1, L = 0;
  void getR(int n) {
    lim = 1, L = 0;
    while (lim <= n) lim <<= 1, L++;
    for (rint i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << L - 1);
  }
  void NTT(int *a, int opt) {
    for (rint i = 0; i < lim; i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (rint mid = 1; mid < lim; mid <<= 1) {
      int Wn = qpow(opt == 1 ? G : Gi, (mod - 1) / (mid << 1));
      for (rint R = mid << 1, j = 0; j < lim; j += R) {
        int w = 1;
        for (rint k = 0; k < mid; k++, w = 1ll * w * Wn % mod) {
          int x = a[j + k], y = 1ll * w * a[j + k + mid] % mod;
          a[j + k] = (x + y) % mod;
          a[j + k + mid] = (x - y + mod) % mod;
        }
      }
    }
    if (opt == -1) {
      int linv = qpow(lim, mod - 2);
      for (rint i = 0; i < lim; i++) a[i] = 1ll * a[i] * linv % mod;
    }
  }
  void mul(int *a, int *b, int n, int m) {
    getR(n + m);
    NTT(a, 1), NTT(b, 1);
    for (rint i = 0; i < lim; i++) a[i] = 1ll * a[i] * b[i] % mod;
    NTT(a, -1);
  }
}

int main() {
  n = read(), m = read();
  rep(i, 0, n) a[i] = read();
  rep(i, 0, m) b[i] = read();
  Poly::mul(a, b, n, m);
  rep(i, 0, n + m) print(a[i]), putchar(‘ ‘);
}

2

2

[算法学习] 多项式全家桶

原文：https://www.cnblogs.com/wlzhouzhuan/p/12657174.html