多项式 学习笔记

我是真的没想到，有很多地方不是独立完成的，但这篇博客仍写了 1 天。

还有它貌似是 22 年的第一篇博客。

多项式的定义与数学中的定义相同，这里指的多项式均为一元多项式。

多项式基本知识

多项式的表达

在计算机中， 多项式有两种常见的表达方式，点值表达和系数表达。

• 系数表达：指用 $n + 1$ 个数来表示一个 $n$ 次多项式。其中 $a_k,k\in[0, n]$ 分别表示该多项式的 $k$ 次项系数。
• 点值表达：代数基本定理中说明：$n$ 个不同点可以表示一个唯一确定的 $n - 1$ 次多项式。这里只讨论用一些有特性的点来表示一个多项式，所以横坐标是可以根据下标推出来的，$a_0, \dots a_n$ 分别表示第 $k$ 个点的纵坐标。

多项式的界

大多数情况下我们只需要在意低次项，一个多项式的界就代表最高有用位。

多项式的四则运算

加减法运算：

• 系数表达下：直接系数相加减即可
• 点值表达下：直接纵坐标相加减即可。

乘法运算：

• 系数表达下：本质上是加法卷积，令 $f[k]$ 表示多项式 $f$ 的 $k$ 次项系数。式子是：

  $$C[k] = \sum_{i + j = k} A[i] \times b[j]$$

• 点值表达下：直接 $c_k = a_k \times b_k$ 即可。

求逆运算：

较为复杂，待会分析。

多项式半家桶

从 $FFT$ 到 $NTT$ ，点值表达与系数表达的转换。

引入

从多项式的乘法就可以看出，某些情况下点值表达比系数表达要方便，但有些情况下又必须用系数表达。所以如何将这两种表达方式串联起来，做到快速转化便是主要问题。

考虑到点值表达和系数表达都可以看作 1 行 n 列的矩阵，转化相当于对他们进行线性变换。如果可以找到一些特殊的点，那么就可以用分治来代替暴力以降低复杂度。

数学界某知名大佬 Fourier 找到了一类符合条件的数：单位根（$\omega$）。

单位根

单位根是指复数域单位圆上的点。$n$ 次单位根 $\omega_n^1$ 本质上是 $\omega_n^1 = \sqrt[n]{1}$ 。

这个东西的某些性质和三角函数很像。（毕竟都是圆上的东西）

具体来说：

1. $\omega_n^m = (\omega_n^1)^{m\%n}$
2. $\omega_n^j \times \omega_n^k = w^{j + k}$
3. $\omega_{2n}^{2m} = w_n^m$
4. $\omega_n^{k+n/2} = -\omega_n^k$

有了上面的性质，就可以暴力推式子了。

DFT 的推导

设 $F(x)$ 的项数是 $2$ 的整次幂（不是的暴力补高位为 0 就行） 将 $F(x)$ 按照奇偶性劈成两半：$F(x) = (F[0] + F[2]x^2 + F[n - 2]x^{n - 2}) + (F[1]x + f[3]x^3 + f[n - 1]x^{n - 1})$ 。

设：

$$Fl(x) = \sum F[2k]x^{k}, Fr(x) = \sum F[2k + 1]x^k ,k \in [0, n / 2 - 1]$$

可以发现：

$$F(x) = Fl(x ^ 2) + x\cdot Fr(x^2)$$

按照上下半圆分别带入单位根：

• 当 $k < n / 2$ 时，代入 $\omega^k_n$

  $$F(\omega^k_n) = Fl(\omega^k_{n / 2}) + \omega^k_n Fr(\omega^k_{n / 2})$$

• 当 $k > n / 2$ 时，设 $t = k - n / 2$ ，代入 $\omega^t_n$

  $$F(w^{t + n / 2}_n) = Fl(w^t_{n / 2}) - w^t_n Fr(\omega^t_{n / 2})$$

单看这个式子并不能看出什么，但是如果无限分治下去就可以得出答案。

IDFT 的实现

本质上是单位根反演，没学那么多也写不了那么详细，直接把式子放在这里：

设 $G = DFT(F)$ ，则 $G[k] = \sum_{i = 0}^ {n - 1} (\omega^k_n)^i F[i]$

$$n\cdot F[k] = \sum_{i = 0}^{n - 1} (\omega^{-k}_n)G[i]$$

细节优化

可以发现 $DFT$ 和 $IDFT$ 分治起来只差一个负号，所以可以共用一个函数，并附带一个 $bool$ 型参数。

同时可以发现每个长度为 $n$ 的多项式对应的每个点的位置在开始之前都是唯一确定且已知的。所以可以线性预处理一下。

代码实现

P3803 【模板】多项式乘法（FFT）

代码实现

// 多项式 学习笔记
// Code By CloudySky
#include <bits/stdc++.h>
// #define int long long
const double Pi = acos(-1);
namespace IO {
inline int read() {
  int x = 0, f = 1; char c = getchar();
  while (c < '0' || c > '9') { if (c == '-') f = -1; c = getchar(); }
  while (c >= '0' && c <= '9') { x = x * 10 + (c ^ 48); c = getchar(); }
  return x * f;
}
void print_n(int x) {
  if (x > 9) print_n(x / 10);
  putchar(x % 10 + '0');
}
inline void print(int x, char s = '\n') {
  if (x < 0) putchar('-'), x = -x;
  print_n(x), putchar(s);
}
}  // namespace IO
using namespace IO;
const int Maxn = 1.5e6 + 10;
using namespace std;

struct CP {
  double x, y;
  CP (double tx = 0, double ty = 0) {
    x = tx, y = ty;
  }
  CP operator+ (CP const& B) const {
    return CP (x + B.x, y + B.y);
  }
  CP operator- (CP const& B) const {
    return CP (x - B.x, y - B.y);
  }
  CP operator* (CP const& B) const {
    return CP (x * B.x - y * B.y, x * B.y + y * B.x);
  }
} f[Maxn << 1], g[Maxn << 1];

int tr[Maxn << 1], n, m;

void fft (CP *f, bool flag) {
  for (int i = 0; i < n; ++i)
    if (i < tr[i]) swap (f[i], f[tr[i]]);
  for (int p = 2; p <= n; p <<= 1) {
    int len = p >> 1;
    CP tG(cos(2 * Pi / p), sin (2 * Pi / p));
    if (!flag) tG.y *= -1;
    for (int k = 0; k < n; k += p) {
      CP buf(1, 0);
      for (int l = k; l < k + len; ++l) {
        CP tt = buf * f[len + l];
        f[len + l] = f[l] - tt;
        f[l] = f[l] + tt;
        buf = buf * tG;
      }
    }
  }
}

signed main() {
  n = read () + 1, m = read () + 1;
  for (int i = 0; i < n; ++i)
    scanf ("%lf", &f[i].x);
  for (int i = 0; i < m; ++i)
    scanf ("%lf", &g[i].x);
  for (m += n, n = 1; n < m; n <<= 1);
  for (int i = 0; i < n; ++i)
    tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0);
  fft (f, 1), fft (g, 1);
  for (int i = 0; i < n; ++i)
    f[i] = f[i] * g[i];
  fft (f, 0);
  for (int i = 0; i < m - 1; ++i)
    printf ("%d ", (int) (f[i].x / n + 0.49));
  return 0;
}

剩余系，原根与 $NTT$

前置知识略解：

• 阶：最小的 $k$ 满足 $a^k \equiv 1\bmod p$，记为 $\delta_pa$

• 原根：最大的 $g$ 满足 $\delta_pg = \varphi(p)$

$FFT$ 好是好，但是由于用到了复数，常数巨大且容易丢精度。

所以就要找一些东西来代替单位根。但是实数域已经找不到可以替代的东西了。但好在大多数时候研究问题是在剩余系中，然后就有东西可以代替了。且一般模数都是 $998244353 = 2^{23} \cdot 7 \cdot 17 + 1$ 考虑到原根的性质，可以支持找到一个单位根满足条件，即 $g^{\small\tfrac{p - 1}{n}}$ 。

代码实现

代码实现：

// 多项式 学习笔记
// NTT 多项式转换
void ntt (int *g, int n, bool op) {
  static ull f[Maxn << 1];                      // 用宏定义来避免命名冲突
  tpre (n);                                     // 预处理顺序
  for (int i = 0; i < n; ++i)
    f[i] = g[tr[i]];                            // 调整数组顺序
  for (int p = 2; p <= n; p <<= 1) {
    int len = p >> 1;
    ull tG = Pow (op ? G : invG, (P - 1) / p);
    for (int k = 0; k < n; k += p) {
      int buf = 1;
      for (int l = k; l < k + len; ++l) {
        int tt = buf * f[l | len] % P;
        f[l | len] = (f[l] - tt + P) % P;
        f[l] = (f[l] + tt + P) % P;
        buf = buf * tG % P;
      }
    }
  }
  if (op == 0) {                                 // 根据 DFT 还是 IDFT 来进行下一步操作
    ull invn = Pow (n);
    for (int i = 0; i < n; ++i) g[i] = f[i] * invn % P;
  } else 
    for (int i = 0; i < n; ++i) g[i] = f[i] % P;
}

多项式的求导，积分

一定要记住，这两个的结果都是一个多项式，而不是一个值

多项式求导（通俗理解就是求斜率）

它的数学求法为 $\forall x1, x2 \in M , \text{令} \Delta x = x2 - x1, \text{则} f'(x1) = \dfrac{f(x1 + \Delta x) - f(x1)}{\Delta x}$

用系数就是：

$$F'(x) = \sum_{i = 0} F[i]x^i\cdot \dfrac{i}{x} \\ \Rightarrow F'(x) = \sum_{i = 0} F[i]x^{i - 1}\cdot i$$

可以得出的结论是一个 $n$ 次多项式的导数是一个 $n - 1$ 次多项式。

多项式积分（通俗理解就是求面积）

$$\int(x) = C + \sum_{i = 0}F[i]x^i \cdot \dfrac{x}{i} \\ \Rightarrow\int(x) = C + \sum_{i = 0} \dfrac{F[i]x^{i + 1}}{i}$$

更好地理解

如果上面不理解，可以类比一下物理学 $x - t， v - t, a - t$ 关系式。

一个 $v - t$ 关系式的积分是 $x - t$ 关系式，导数是 $a - t$ 关系式。如果 $v - t$ 为 1 次，那么 $x - t$ 为 2 次；如果 $v - t$ 关系式为 2 次，那么 $a - t$ 为 1 次。

代码实现

代码实现：

// 多项式 学习笔记
// 多项式求导
void der (int *f, int n) {
  for (int i = 1; i < n; ++i) {
    f[i - 1] = 1ll * f[i] * i % P;
  } f[n - 1] = 0;
}
// 多项式积分
void inte (int *f, int n) {
  for (int i = n; i; --i) {
    f[i] = 1ll * f[i - 1] * Pow (i, P - 2) % P;
  } f[0] = 0;
}

多项式的复合 - 牛顿迭代

这个东西不会证。他不是某个函数，更像是帮你推导式子的一种工具，需要配合倍增使用。

具体来说就是当 $G$ 已知且 $G(f(x)) = 0$ 时，$F(x) \equiv F_*(x) - \dfrac{G(F_*(x))}{G'(F_*(x))}\bmod x^n$ 。

其中 $F_*(x)$ 表示上一次倍增的结果，界为 $F(x)$ 的 $\dfrac{1}{2}$ 的多项式，$G$ 一般是构造出来的。

多项式求逆

这个东西有三种求法：朴素求法，倍增大力推式子，倍增牛顿迭代。

朴素求法

考虑逆元的性质，设 $F(x)$ 为已知多项式 $G(x)$ 为要求的逆多项式。

$$C(k) = \sum_{i = 0}^k F[i]G[k - i]$$

只有 $C$ 的第 $0$ 位为 1 。即当目前不是常数时

$$\sum_{i = 0}^n F[i]G[n - i] = 0$$

将 $F[0]G[n]$ 单拎出来，

$$\sum_{i = 1}^n F[i]G[n - i] + F[0]G[n] = 0$$

移向消元：

$$G[n] = -\dfrac{1}{F[0]}\cdot \sum_{i = 1}^n F[i]G[n - i]$$

这样每次求末项，逐步求出的多项式 $G$ 就是 $F$ 的逆多项式。

倍增解法

考虑当前要求界为 $n$ 意义下的逆多项式 $R(x)$ ，已经求出了界为 $\dfrac{n}{2}$ 的逆多项式 $R_*(x)$ 。

则当前两个多项式满足

$$R(x) \equiv R_*(x)\bmod x^{\tfrac{n}{2}}$$

移项

$$R(x) - R_*(x) \equiv 0 \bmod x^ {\tfrac{n}{2}}$$

乘方，扩界

$$(R(x) - R_*(x))^2 \equiv 0 \bmod x^n$$

拆括号，移相，两边同乘 $F(x)$

$$R(x) \equiv 2R_*(x) - R_*^2F(x) \bmod x^n$$

牛顿迭代解法

这个方法推出来的式子和上面是一样的，只是换一种推导思路而已。

$$F(x)R(x) - 1 \equiv 0 \bmod P$$

设 $F(x)R(x) - 1$ 为 $G$ ， 其中 $F(x)$ 已知，可以作为一个系数。下面逐步推导即可。

代码实现

代码实现：

// 多项式 学习笔记
// 多项式求逆
void invp (int *f, int n) {
  int m = n;
  static int w[Maxn << 1], r[Maxn << 1];
  for (n = 1; n < m; n <<= 1);
  w[0] = Pow (f[0]);
  for (int len = 2; len <= n; len <<= 1) {
    for (int i = 0; i < (len >> 1); ++i)
      r[i] = w[i];
    cpy (sav, f, len);
    ntt (sav, len, 1), ntt (r, len, 1);
    px (r, sav, len), ntt (r, len, 0);
    clear (r, 0, len >> 1);
    cpy (sav, w, len);
    ntt (sav, len, 1), ntt (r, len, 1);
    px (r, sav, len), ntt (r, len, 0);
    for (int i = len >> 1; i < len; ++i)
      w[i] = (2ll * w[i] - r[i] + P) % P;
  }
  cpy (f, w, m);
  clear (sav, 0, n), clear (w, 0, n), clear (r, 0, n);
}

本文作者：CloudySky
写作时间： 2022-01-09
最后更新时间： 2022-01-09
遵循协议： BY-NC-SA