多项式合集

Contents

1 多项式合集
- 1.1 拉格朗日插值
- 1.2 快速傅里叶变换（FFT)

多项式合集

[TOC]

拉格朗日插值

问题背景

给出 $n$ 个点 $(x_i,y_i)$ ，令这 $n$ 个点确定的多项式为 $L(x)$ ，求 $L(k)\bmod 998244353$ 的值。

结论

L(x) = \sum_{i=1}^n y_il_i(x)

其中每个 $l_i(x)$ 为拉格朗日基本多项式，表达式为

l_i(x) = \prod_{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}

其特点是 $l_i(x_i)=1$ ， $\forall j\ne i$ 有 $l_i(x_j)=0$

推导

抛开拉插，这道题明显可以列方程组然后使用高斯消元求解，但是复杂度为 $O(n^3)$ 且精度问题明显，所以拉格朗日是这样考虑的：

对于每个点 $P_i(x_i,y_i)$ ，构造一个 $n-1$ 次多项式 $l_i(x)$ 使其在 $x_i$ 上取值为 $1$ ，在其余 $x_j$ 上为 $0$ 。构造的结果就是上面的结论：

l_i(x) = \prod_{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}

这个多项式的正确性还是很显然的。然后我们也知道这个多项式它就是唯一的。

然后考虑构造答案：很显然对于点 $P_i(x_i,y_i)$ ，只有 $l_i(x_i)$ 的取值为 $1$ ，其他的都为 $0$ 。所以答案的正确性也是比较显然的：对于 $x_i$ ，只有 $y_il_i(x_i)$ 产生了贡献，其余的都是 $0$ 。故这个多项式是正确的。

所以回到一开始，我们需要的就是

f(k) = \sum_{i=1}^n y_i\prod_{j=1,j\ne i}^n\frac{k-x_j}{x_i-x_j}

由于模数是质数，所以使用费马小定理求逆元，跑得飞快。

复杂度 $O(n^2)$ ，求逆元就是个很小的常数

#include <cstdio>
#include <cctype>
#define il inline

typedef long long ll;

inline ll read()
{
    char c = getchar();
    ll s = 0;
    bool x = 0;
    while (!isdigit(c))
        x = x | (c == '-'), c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return x ? -s : s;
}

const ll maxn = 2e3 + 5, mod = 998244353;
ll x[maxn], y[maxn];

ll pow(ll base, ll p)
{
    ll ans = 1;
    base = (base + mod) % mod;
    for (; p; p >>= 1)
    {
        if (p & 1)
            ans = ans * base % mod;
        base = base * base % mod;
    }
    return ans;
}

il ll inv(ll n)
{
    return pow(n, mod - 2);
}

int main()
{
    ll n = read(), k = read();
    for (int i = 1 ; i <= n; ++i)
        x[i] = read(), y[i] = read();
    ll ans = 0;
    for (int i = 1; i <= n; ++i)
    {
        ll prod1 = 1, prod2 = 1;
        for (int j = 1; j <= n; ++j)
        {
            if (i == j)
                continue;
            prod1 = prod1 * (k - x[j]) % mod;
            prod2 = prod2 * (x[i] - x[j]) % mod;
        }
        ans = (ans + prod1 * y[i] % mod * inv(prod2) % mod + mod) % mod;
    }
    printf("%lld\n", ans);
    return 0;
}

拉格朗日插值与范德蒙矩阵

可以考虑将这 $n+1$ 个点值表示为如下形式：

\begin{bmatrix} x_0^0 & x_0^1 & x_0^2 &\cdots &x_0^n\\ x_1^0 & x_1^1 & x_1^2 &\cdots & x_1^n\\ \vdots & \vdots & \vdots & &\vdots\\ x_n^0 & x_n^1 & x_n^2 & \cdots & x_n^n \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0\\y_1\\ \vdots\\ y_n \end{bmatrix}

左边这个矩阵就是所谓的范德蒙德矩阵，记作 $\boldsymbol V$ ，系数列向量记作 $\boldsymbol A$ ，右边的记作 $\boldsymbol B$ ，则很明显：

\boldsymbol{VA} = \boldsymbol B

打开来看清楚些实际就是多项式 $f$ 在每个点处的值：

y_j = f(x_j) = \sum_{i = 0}^na_ix_j^i

我们把两边都乘上 $\boldsymbol V^{-1}$ ：

\begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix} = \boldsymbol V^{-1} \begin{bmatrix} y_0\\y_1\\ \vdots\\ y_n \end{bmatrix}

就得到了 $a_i$ 一定可以表示为某种形如

a_k = \sum \begin{bmatrix} \vdots \end{bmatrix}y_k

的形式，即 $a_k$ 只与 $x_i$ 与 $y_k$ 有关。

所以不难发现对于一个要求的 $f(x_\ominus)$ ，都可以被表示为如下形式

f(x_\ominus)=\sum\delta_k(x_\ominus)y_k

$\delta_k(x)$ 构造的过程即需要考虑 $x=x_k$ 时 $\delta_j(x) = 0\land\delta_k(x) = 1$ ，其中 $k\not=j$ 。 $\delta_j(x_k) = 0$ 说明每一个 $\delta_j$ 都要有 $(x-x_k)$ 这个因式，然后又因为 $\delta_k(x_k) = 1$ ，所以最终构造出来就是上面的结果：

f(x) = \sum_{i=1}^n y_i\prod_{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}

我们其实也可以利用拉格朗日插值来求范德蒙矩阵的逆阵，复杂度 $O(n^2)$

快速傅里叶变换（FFT)

复数基础

定义虚数单位 $\mathrm i^2 = \sqrt{-1}$ ，把形如 $a + b\mathrm i\:(a,b\in\mathbb R)$ 的数称为复数，所有复数的集合称为复数集 $\mathbb C$ 。

复数一般使用 $z$ 表示，表示为 $z = a + b\mathrm i$ ，这种形式称为复数的代数形式。 $a$ 被称为复数的实部， $b$ 称为复数的虚部，未加说明的情况下一般认为 $a,b\in\mathbb R$ 。很明显地，当 $a = 0\land b\not=0$ 时，这个复数为纯虚数，当 $b=0$ 时，这个复数为实数。

每个复数 $a + b\mathrm i$ 都能对应平面直角坐标系里面的一个点 $(a,b)$ ，同样的也可以对应一个向量 $(a,b)$ 。故定义复数的模为 $\sqrt{a^2 + b^2}$ 。

定义复数的加法与乘法：

\begin{aligned} &(a + b\mathrm i) + (c + d\mathrm i)\\ =&(a + c) + (b + d)\mathrm i \end{aligned}

\begin{aligned} &(a+b\mathrm i)(c + d\mathrm i)\\ =&ac + ad\mathrm i + cb\mathrm i + bd\mathrm i^2\\ =&(ac - bd) + (ad + bc)\mathrm i \end{aligned}

这都是比较显然的。

容易看出复数满足很多实数的运算律。

定义复数 $z=a+b\mathrm i$ 的共轭复数为 $\overline{z} = a - b\mathrm i$ ，不难发现 $z$ 与 $\overline{z}$ 关于实轴对称。

z\overline z=(a+b\mathrm i)(a-b\mathrm i) = a^2 + b^2=|z|^2

复数既然可以对应平面直角坐标系中的向量，不难发现其可以使用其模长与辐角来表示：

z=a+b\mathrm i\iff z = r(\cos\theta+\mathrm i\sin\theta)

其中 $r$ 为 $z$ 的模长， $\theta$ 为其辐角。即我们可以把一个复数表示成二元组 $(r,\theta)$ 的形式。

现在考虑两个复数 $(r_1,\theta_1)$ 与 $(r_2,\theta_2)$ 相乘得到的结果：

\begin{aligned} (r_1,\theta_1)\times(r_2,\theta_2) &= r_1(\cos\theta_1 + \mathrm i\sin\theta_1)r_2(\cos\theta_2 + \mathrm i\sin\theta_2)\\ &=r_1r_2(\cos\theta_1\cos\theta_2 - \sin\theta_1\sin\theta_2 + \mathrm i\sin\theta_1\cos\theta_2 + \mathrm i\sin\theta_2\cos\theta_1)\\ &=r_1r_2\left(\cos(\theta_1 + \theta_2) + \mathrm i\sin(\theta_1 + \theta_2)\right)\\ &=(r_1r_2,\theta_1 + \theta_2) \end{aligned}

于是我们可以概括复数乘法的法则：模长相乘，辐角相加。（上述推导需要掌握基本的三角恒等变换）

从欧拉公式到单位圆

给出复数指数幂的定义：

\mathrm e^{x +y\mathrm i} = e^x(\cos y + \mathrm i\sin y)

这个公式是由我也不会证明的泰勒展开推导出来的：

\begin{aligned} \sin(x) &= x - \frac{x^3}{3!}+\frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} + \cdots = \sum_{k = 1}^\infin\frac{(-1)^{k - 1}x^{2k - 1}}{(2k-1)!}\\ \cos(x) &= 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!} + \frac{x^8}{8!} + \cdots = \sum_{k = 0}^\infin\frac{(-1)^{k} x^{2k}}{(2k)!}\\ \mathrm e^x &= 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots = \sum_{k = 0}^\infin\frac{x^k}{k!} \end{aligned}

将 $x + y\mathrm i$ 代入进去即可推导。

如果 $x = 0$ ，我们就得到大名鼎鼎的欧拉公式：

\mathrm e^{x\mathrm i} = \cos x + \mathrm i\sin x

更特殊地，如果 $x = \pi$ ，得到的就是下面这个神奇的式子：

\mathrm e^{\pi\mathrm i} = -1

复平面上我们可以定义类似于平面直角坐标系上的单位圆，单位圆上的所有复数构成集合 $\{z: |z| = 1\}$ 。这些复数都可以表示为 $\cos\theta + \mathrm i\sin\theta$ 或 $e^{\theta \mathrm i}$ 的形式。

多项式的表示法

系数表示法：顾名思义

f(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n\iff f(x) = \{a_0,a_1,a_2,\cdots,a_n\} = \begin{bmatrix} x^0 & x^1 & x^2 &\cdots & x^n \end{bmatrix} \begin{bmatrix} a_0\\a_1\\a_2\\\vdots\\a_n \end{bmatrix}

点值表示法：

我们知道由一个多项式在 $n + 1$ 个点上的取值是可以唯一确定一个多项式的，其本质也就是线性方程组的解。所以一个 $n$ 次多项式可以用 $n + 1$ 个点表示：

f(x) = \{(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\}

\begin{bmatrix} x_0^0 & x_0^1 & x_0^2 &\cdots &x_0^n\\ x_1^0 & x_1^1 & x_1^2 &\cdots & x_1^n\\ \vdots & \vdots & \vdots & &\vdots\\ x_n^0 & x_n^1 & x_n^2 & \cdots & x_n^n \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0\\y_1\\ \vdots\\ y_n \end{bmatrix}

通过下面的这个形式我们看得出来其就是一个典型的线性方程组的形式，不难证明其解的唯一性。

并且我们发现点值表示法有一个很明显的优势：可以在 $O(n)$ 的时间内将两个多项式乘起来，只需把对应点的 $y$ 乘起来即可。

通俗的来说，FFT 实现的就是快速求多项式乘法的过程：先把系数表示法转成点值表示法（DFT，离散傅里叶变换），乘完之后再把点值还原为插值（IDFT，离散傅里叶逆变换）。可是朴素的 DFT 需要的时间复杂度为 $O(n^2)$ ，IDFT 还回其系数需要高斯消元是 $O(n^3)$ 的。而 FFT 利用了一些很特殊很特殊的值加速了 DFT 和 IDFT 的过程，使得总时间复杂度降低到了 $O(n\log n)$ 。

单位复数根

解这个方程：

x^n = 1

我们会发现这个方程在实数范围内只有 $1$ 或者 $2$ 个解。然而代数基本定理告诉我们这样的方程有 $n$ 个复数域上的解。由模长相乘辐角相加我们知道因为最终 $x^n = 1$ ，所以这些满足条件的 $x$ 的模长必定也是 $1$ 。然后需要满足他们的辐角的 $n$ 倍能被 $2\pi$ 整除。

不难发现其就是 $n$ 等分单位圆：

我们记 $n$ 次单位根的第 $k$ 个记为 $\omega_n^k$ ，不难发现 $\omega_k^n = \mathrm e^{\frac{2k\pi i}{n}}$ 。由此可见，单位复数根具有一些非常好的性质比如：

\begin{aligned} \omega_n^0 = \omega_n^n &= 1\\ \omega_n^k &= \omega_{2n}^{2k}\\ \omega_{2n}^{k + n} &= -\omega_{2n}^k\\ \left(\omega_{2n}^{k + n}\right)^2 &=\omega_n^k \end{aligned}

利用这些性质，我们可以加速 DFT 的过程。FFT 就是利用分治思想加速求每个 $f(\omega_n^k)$ 的值

DFT

此时 DFT 的分治思想就是分开考虑奇次项和偶次项：

考虑

f(x) = a_0x^0 + a_1x^1 + a_2x^2 + \cdots

将其分为两个多项式

\begin{aligned} f(x) &= a_0x^0 + a_2x^2 + a_4x^4 + a_6x^6 + a_8x^8 + \cdots +a_1x^1 + a_3x^3 + a_5x^5 + a_7x^7 + a_9x^9 + \cdots\\ &= a_0x^0 + a_2x^2 + a_4x^4 + a_6x^6 + a_8x^8+\cdots +x(a_1x^0 + a_3x^2 + a_5x^4 + a_7x^6 + \cdots) \end{aligned}

考虑两个新多项式：

\begin{aligned} f_0(x) &= a_0x^0 + a_2x^1 + a_4x^2 + a_6x^3 + \cdots\\ f_1(x) &= a_1x^0 + a_3x^1 + a_5x^2 + a_7x^3 + \cdots \end{aligned}

不难发现

f(x) = f_0(x^2) + xf_1(x^2)

利用单位复数根的性质：

\begin{aligned} \mathrm{DFT}(f(\omega_n^k)) &= \mathrm{DFT}(f_0(\omega_n^{2k})) + \omega_n^k\mathrm{DFT}(f_1(\omega_n^{2k}))\\ &=\mathrm{DFT}(f_0(\omega_\frac n2^k)) + \omega_n^k\mathrm{DFT}(f_1(\omega_\frac n2^k)) \end{aligned}

\begin{aligned} \mathrm{DFT}(f(\omega_n^{k + \frac n2})) &= \mathrm{DFT}(f_0(\omega_n^{2k + n})) + \omega_{n}^{k + \frac n2}\mathrm{DFT}(f_1(\omega_n^{2k + n}))\\ &=\mathrm{DFT}(f_0(\omega_n^n\omega_n^{2k})) - \omega_n^k\mathrm{DFT}(f_1(\omega_n^n\omega_n^{2k}))\\ &=\mathrm{DFT}(f_0(\omega_\frac n2^k)) - \omega_n^k\mathrm{DFT}(f_1(\omega_\frac n2^k)) \end{aligned}

其中 $k < \displaystyle\frac n2$ 。不难发现只要我们求得出 $\mathrm{DFT}(f_0(\omega_\frac n2^k))$ 与 $\mathrm{DFT}(f_1(\omega_\frac n2^k))$ 的话，就可以同时求出 $\mathrm{DFT}(f(\omega_n^k))$ 和 $\mathrm{DFT}(f(\omega_n^{k + \frac n2}))$ 。接下来再对 $f_0$ 与 $f_1$ 递归 DFT 即可。其时间复杂度函数是形如下面这样的：

T(n) = T(n/2) + O(n)

所以总复杂度为 $\Theta(n\log n)$

实际实现的时候一定要注意传进去的系数一定要是 $2^m$ 个的，不然分治的过程中左右不一样会出问题。第一次传进去的时候就高位补 $0$ ，补成最高项次数为 $2^{m - 1}$ 的多项式。

void dft(int lim, complex *a)
{
    if (lim == 1) return;//常数项直接返回
    complex a1[lim >> 1], a2[lim >> 1];
    for (int i = 0; i < lim; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];//把系数按照奇偶分开
    dft(lim >> 1, a1, type);//求 DFT(f_0())
    dft(lim >> 1, a2, type);//求 DFT(f_1())
    complex Wn = complex(cos(2.0 * pi / lim), sin(2.0 * pi / lim)), w = complex(1, 0);
    for (int i = 0; i < (lim >> 1); ++i, w = w * Wn)
    {
        a[i] = a1[i] + w * a2[i];//求 DFT(f(\omega_n^k))
        a[i + (lim >> 1)] = a1[i] - w * a2[i];//求 DFT(f(\omega_n^{k+\fracn2}))
    }
    return;
}

IDFT

好了现在假装我们已经求出了两个多项式的点值表达并已经将他们乘起来，但是我们最终还是要把他还原回去到系数表示的。这个过程就叫做 IDFT。

其实就是我们需要求解下面关于 $a$ 的线性方程组：

\begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 &\cdots &(\omega_n^0)^n\\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 &\cdots & (\omega_n^1)^n\\ \vdots & \vdots & \vdots & &\vdots\\ (\omega_n^{n})^0 & (\omega_n^{n})^1 & (\omega_n^{n})^2 & \cdots & (\omega_n^n)^n \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0\\y_1\\ \vdots\\ y_n \end{bmatrix}

我们将其乘上左边矩阵的逆：

\begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 &\cdots &(\omega_n^0)^n\\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 &\cdots & (\omega_n^1)^n\\ \vdots & \vdots & \vdots & &\vdots\\ (\omega_n^{n})^0 & (\omega_n^{n})^1 & (\omega_n^{n})^2 & \cdots & (\omega_n^n)^n \end{bmatrix}^{-1} \begin{bmatrix} y_0\\y_1\\ \vdots\\ y_n \end{bmatrix}

模相同的正交列向量构成的矩阵的逆是转置的模分之一倍，所以：

\begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 &\cdots &(\omega_n^0)^n\\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 &\cdots & (\omega_n^1)^n\\ \vdots & \vdots & \vdots & &\vdots\\ (\omega_n^{n})^0 & (\omega_n^{n})^1 & (\omega_n^{n})^2 & \cdots & (\omega_n^n)^n \end{bmatrix}^{-1} =\frac {1}{n+1} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & (\omega_n^{-0})^2 &\cdots &(\omega_n^{-0})^n\\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & (\omega_n^{-1})^2 &\cdots & (\omega_n^{-1})^n\\ \vdots & \vdots & \vdots & &\vdots\\ (\omega_n^{-n})^0 & (\omega_n^{-n})^1 & (\omega_n^{-n})^2 & \cdots & (\omega_n^{-n})^n \end{bmatrix}

所以不难发现，IDFT 其实就是再做了一遍 DFT，只不过是反起来的。只是算出来最后的系数结果都要除以点值的个数，反应在代码里面就是那个 lim 变量。

不难发现 $\omega_n^k$ 的共轭就是虚部取反，所以可以在 DFT 函数里面传一个参数表示是否为 IDFT。

这样子一个递归版的 FFT 就写完了，总体的代码如下：

#include <cstdio>
#include <cctype>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

const int maxn = 2e6 + 5;
const double pi = acos(-1.0);

inline int read()
{
    char c = getchar();
    int s = 0;
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return s;
}

struct complex
{
    double x, y;
    complex(double xx = 0, double yy = 0)
    {
        x = xx, y = yy;
    }
} a[maxn], b[maxn];

complex operator+(const complex &a, const complex &b) {return complex(a.x + b.x, a.y + b.y);}
complex operator-(const complex &a, const complex &b) {return complex(a.x - b.x, a.y - b.y);}
complex operator*(const complex &a, const complex &b) {return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void dft(int lim, complex *a, int type)//type = 1 DFT；type = -1 IDFT
{
    if (lim == 1) return;//返回常数项
    complex a1[lim >> 1], a2[lim >> 1];
    for (int i = 0; i < lim; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    dft(lim >> 1, a1, type);
    dft(lim >> 1, a2, type);
    complex Wn = complex(cos(2.0 * pi / lim), type * sin(2.0 * pi / lim)), w = complex(1, 0);
    for (int i = 0; i < (lim >> 1); ++i, w = w * Wn)
    {
        a[i] = a1[i] + w * a2[i];
        a[i + (lim >> 1)] = a1[i] - w * a2[i];
    }
    return;
}

int main()
{
    int n = read(), m = read();
    FOR(i, 0, n) a[i].x = read();
    FOR(i, 0, m) b[i].x = read();
    int lim = 1;
    while (lim <= n + m) lim <<= 1;//lim一定要大于 n + m
    dft(lim, a, 1);
    dft(lim, b, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i];//点值乘起来
    dft(lim, a, -1);//IDFT还回去
    FOR(i, 0, n + m)
        printf("%d ", (int)(a[i].x / lim + 0.5));//最后要除那个数然后还原回去，四舍五入
    return 0;
}

位逆序置换

然而，上面的代码连模板都跑不过去……

考虑继续优化 DFT 的过程。递归的过程中开了大量的空间并且常数巨大，考虑非递归写法。

只考虑我们对 $0$ 到 $7$ 操作：

递归的过程：

original     0   1   2   3   4   5   6   7
recursion#1     0   2   4   6   1   3   5   7
recursion#2     0   4   2   6   1   5   3   7
recursion#3     0   4   2   6   1   5   3   7
original bin    000 001 010 011 100 101 110 111
now bin         000 100 010 110 001 101 011 111

可见递归到最后的结果无非就是一个二进制反转。

所以我们可以考虑非递归，一开始就先把所有的数放到最后的位置，然后迭代的时候一步步还回去即可。这个过程就是位逆序置换（蝴蝶变换）

考虑处理出 $x$ 二进制位翻转之后的数 $R(x)$ 。易知 $R(0) = 0$ 。我们可以从小到大求 $R(x)$ 。很明显， $\lfloor x/2\rfloor$ 的二进制位是 $x$ 右移一位，那么如果知道了 $R(\lfloor x/2\rfloor)$ 就可以很容易的求出 $R(x)$ ，再分 $x$ 的奇偶性判断就可以了。

R(x) = \left\lfloor\frac{R(\lfloor x/2\rfloor)}{2}\right\rfloor + (x\bmod 2)\times\frac{len}2

举个例子：翻转 $(10101110)_2$ ，首先我们知道它的二分之一倍为 $(01010111)_2$ ，其翻转结果为 $(11101010)_2$ ，除以二变为 $(01110101)_2$ ，由于它是偶数所以前面不用补 $1$ 。不难发现其就是一开始要求的翻转结果。

预处理翻转结果的代码：

while (lim <= n + m) lim <<= 1;
FOR(i, 0, lim - 1)
    rev[i] = ((rev[i >> 1] >> 1) | (((i & 1) ? (lim >> 1) : 0)));

然后在处理翻转的时候只需要下面几行：

FOR(i, 0, lim - 1)
    if (i < rev[i])
        myswap(a[i], a[rev[i]]);

不难验证其正确性。

而且观察我们在求 $\mathrm{DFT}(f(\omega_n^k))$ 时我们需要算两遍 $\omega_n^k\mathrm{DFT}(f_1(\omega_\frac n2^k))$ ，复数的乘法常数很大，考虑使用临时变量记录以降低常数。

这样子的话迭代版的 DFT 过程就很好写了：

void DFT(int lim, complex *a, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(a[i], a[rev[i]]);//先预处理翻转完了的结果
    for (int p = 2; p <= lim; p <<= 1)//模拟合并答案的过程，即为所谓的 n
    {
        int len = p >> 1;//即上面的 n / 2
        complex Wp = complex(cos(2 * pi / p), type * sin(2 * pi / p));//处理出 p 次单位根
        for (int k = 0; k < lim; k += p)//对每一个进行合并
        {
            complex w = complex(1, 0);//处理 \omega_p^0
            for (int l = k; l < k + len; ++l, w = w * Wp)//开始合并
            {
                //此时的 a[l] 就是之前的 a1[i]，a[len + l] 就是之前的 a2[i]
                complex tmp = w * a[len + l];
                a[len + l] = a[l] - tmp;//相当于上面的 a[i + (lim >> 1)] = a1[i] - w * a2[i]
                a[l] = a[l] + tmp;//相当于上面的 a[i] = a1[i] + w * a2[i]
            }
        }
    }
}

例题

总的一个非递归版 FFT 的实现如下（洛谷 P3803）：

#include <cstdio>
#include <cctype>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

const int maxn = 3e6 + 5;
const double pi = acos(-1.0);

inline int read()
{
    char c = getchar();
    int s = 0;
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
        s = 10 * s + c - '0', c = getchar();
    return s;
}

template<typename T> inline void myswap(T &a, T &b)
{
    T t = a;
    a = b;
    b = t;
    return;
}

struct complex
{
    double x, y;
    complex(double xx = 0, double yy = 0)
    {
        x = xx, y = yy;
    }
} a[maxn], b[maxn];

int rev[maxn];

complex operator+(const complex &a, const complex &b) {return complex(a.x + b.x, a.y + b.y);}
complex operator-(const complex &a, const complex &b) {return complex(a.x - b.x, a.y - b.y);}
complex operator*(const complex &a, const complex &b) {return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void DFT(int lim, complex *a, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(a[i], a[rev[i]]);//先预处理翻转完了的结果
    for (int p = 2; p <= lim; p <<= 1)//模拟合并答案的过程，即为所谓的 n
    {
        int len = p >> 1;//即上面的 n / 2
        complex Wp = complex(cos(2 * pi / p), type * sin(2 * pi / p));//处理出 p 次单位根
        for (int k = 0; k < lim; k += p)//对每一个进行合并
        {
            complex w = complex(1, 0);//处理 \omega_p^0
            for (int l = k; l < k + len; ++l, w = w * Wp)//开始合并
            {
                //此时的 a[l] 就是之前的 a1[i]，a[len + l] 就是之前的 a2[i]
                complex tmp = w * a[len + l];
                a[len + l] = a[l] - tmp;//相当于上面的 a[i + (lim >> 1)] = a1[i] - w * a2[i]
                a[l] = a[l] + tmp;//相当于上面的 a[i] = a1[i] + w * a2[i]
            }
        }
    }
}

int main()
{
    int n = read(), m = read();
    FOR(i, 0, n) a[i].x = read();
    FOR(i, 0, m) b[i].x = read();
    int lim = 1;
    while (lim <= n + m) lim <<= 1;//补齐高位
    FOR(i, 0, lim - 1)
        rev[i] = ((rev[i >> 1] >> 1) | (((i & 1) ? (lim >> 1) : 0)));//先处理翻转完的结果
    DFT(lim, a, 1);//DFT
    DFT(lim, b, 1);//DFT
    FOR(i, 0, lim)
        a[i] = a[i] * b[i];//对处理出来的点值进行乘法
    DFT(lim, a, -1);//IDFT
    FOR(i, 0, n + m)
        printf("%d ", (int)(a[i].x / lim + 0.5));
    return 0;
}

使用 FFT 来求高精度整数乘法的实现（洛谷 P1919）：

#include <cstdio>
#include <cstring>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)
#define DEC(i, a, b) for (int i = a; i >= b; --i)

template<typename T> inline void myswap(T &a, T &b) {T t = a; a = b; b = t; return;}

typedef double db;

const int maxn = 3000000 + 5;
const db pi = acos(-1.0);

struct cmplx
{
    db x, y;
    cmplx(db xx = 0, db yy = 0) {x = xx, y = yy;}
} a[maxn], b[maxn];

cmplx operator+(const cmplx &a, const cmplx &b) {return cmplx(a.x + b.x, a.y + b.y);}
cmplx operator-(const cmplx &a, const cmplx &b) {return cmplx(a.x - b.x, a.y - b.y);}
cmplx operator*(const cmplx &a, const cmplx &b) {return cmplx(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

char s1[maxn], s2[maxn];
int rev[maxn], ans[maxn];

void DFT(cmplx *f, int lim, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            myswap(f[i], f[rev[i]]);
    for (int p = 2; p <= lim; p <<= 1)
    {
        int len = p >> 1;
        cmplx Wp(cos(2.0 * pi / p), type * sin(2.0 * pi / p));
        for (int k = 0; k < lim; k += p)
        {
            cmplx w(1, 0);
            for (int l = k; l < k + len; ++l, w = w * Wp)
            {
                cmplx tmp = w * f[l + len];
                f[l + len] = f[l] - tmp;
                f[l] = f[l] + tmp;
            }
        }
    }
    return;
}

int main()
{
    scanf("%s\n%s", s1, s2);
    int n1 = -1, n2 = -1;
    DEC(i, strlen(s1) - 1, 0) a[++n1].x = s1[i] - '0';
    DEC(i, strlen(s2) - 1, 0) b[++n2].x = s2[i] - '0';
    int lim = 1;
    while (lim <= n1 + n2) lim <<= 1;
    FOR(i, 0, lim - 1)
        rev[i] = ((rev[i >> 1] >> 1) | (((i & 1) ? (lim >> 1) : 0)));
    DFT(a, lim, 1);
    DFT(b, lim, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i];
    DFT(a, lim, -1);
    FOR(i, 0, lim)
        ans[i] = (int)(a[i].x / lim + 0.5);
    FOR(i, 0, lim)
        if (ans[i] >= 10) ans[i + 1] += ans[i] / 10, ans[i] %= 10, lim += (i == lim);
    while (!ans[lim] && lim > -1) --lim;
    if (lim == -1) puts("0");
    else DEC(i, lim, 0) printf("%d", ans[i]);
    return 0;
}

当然，千万要记得 IDFT 还回去的时候要除以 lim，实在怕记不住就在 DFT 函数里面加几句话直接处理好

if (type == -1)
    FOR(i, 0, lim - 1)
        f[i].x /= lim;

可以应用 FFT 的一道题：洛谷 P3338 [ZJOI2014]力

题意：给定 $\{q\}$ ，定义

F_i = \sum_{j = 1}^{i - 1}\frac{q_iq_j}{(i - j)^2} - \sum_{j = i + 1}^n\frac{q_iq_j}{(i - j)^2}

求

E_i=\frac{F_i}{q_i}

考虑暴力的话，这道题是 $O(n^2)$ 的，过不去，考虑转化式子：

\begin{aligned} E_i &= \frac{F_i}{q_i}\\ &=\sum_{j = 1}^{i - 1}\frac{q_j}{(i - j)^2} - \sum_{j = i + 1}^n\frac{q_j}{(i - j)^2}\\ \end{aligned}

我们尝试将其化为卷积的形式，令 $f_i = q_i$ ，且 $f_0 = 0$ ； $g_i =\dfrac{1}{i^2}$ ，且 $g_0 = 0$ ，回代：

\begin{aligned} E_i &= \sum_{j = 0}^{i}f_jg_{i - j} - \sum_{j = i}^nf_jg_{j - i} \end{aligned}

左边的部分已经是一个卷积的形式了，考虑继续化简右边。此时我们可以使用一个翻转的技巧，令 $f'_i = f_{n - i}$ ， $t = n - i$ ，则右半边的式子可以化为 $\displaystyle\sum_{j = 0}^{t}f'_{t - j}g_j$ 。现在两边都化为卷积的形式了，可以愉快的使用 FFT 加速了。

即我们设多项式 $A(x) =\displaystyle\sum_{i = 0}^nf_ix^n$ ， $B(x) = \displaystyle\sum_{i = 0}^ng_ix^n$ ， $C(x) = \displaystyle\sum_{i = 0}^nf'_i$ 。再令 $L(x) = A(x)B(x)$ ， $R(x) = B(x)C(x)$ ，不难发现答案 $E_i = l_i - r_{n - i}$ ，其中 $l_i$ 和 $r_i$ 分别为 $L(x)$ 和 $R(x)$ 中 $x^i$ 的系数。

#include <cstdio>
#include <cmath>
#define FOR(i, a, b) for (int i = a; i <= b; ++i)

typedef double db;

const int maxn = 3e5 + 5;
const db pi = acos(-1.0);

struct cmplx
{
    db x, y;
    cmplx(db xx = 0, db yy = 0) {x = xx, y = yy;}
} a[maxn], b[maxn], c[maxn];

cmplx operator+(const cmplx &a, const cmplx &b) {return cmplx(a.x + b.x, a.y + b.y);}
cmplx operator-(const cmplx &a, const cmplx &b) {return cmplx(a.x - b.x, a.y - b.y);}
cmplx operator*(const cmplx &a, const cmplx &b) {return cmplx(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

int rev[maxn];

template<typename T> inline void swap(T &a, T &b) {T t = a; a = b; b = t; return;}

void DFT(cmplx *f, int lim, int type)
{
    FOR(i, 0, lim - 1)
        if (i < rev[i])
            swap(f[i], f[rev[i]]);
    for (int p = 2; p <= lim; p <<= 1)
    {
        int len = p >> 1;
        cmplx Wp(cos(pi / len), type * sin(pi / len));
        for (int k = 0; k < lim; k += p)
        {
            cmplx w(1, 0);
            for (int l = k; l < k + len; ++l, w = w * Wp)
            {
                cmplx tmp = w * f[l + len];
                f[l + len] = f[l] - tmp;
                f[l] = f[l] + tmp;
            }
        }
    }
    if (type == -1)
        FOR(i, 0, lim - 1)
            f[i].x /= lim;
    return;
}

int main()
{
    int n; scanf("%d", &n);
    FOR(i, 1, n)
    {
        scanf("%lf", &a[i].x);
        b[i].x = (1.0 / i / i);
        c[n - i].x = a[i].x;
    }
    int lim = 1;
    while (lim <= (n << 1)) lim <<= 1;
    FOR(i, 0, lim)
        rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0));
    DFT(a, lim, 1), DFT(b, lim, 1), DFT(c, lim, 1);
    FOR(i, 0, lim)
        a[i] = a[i] * b[i], c[i] = b[i] * c[i];
    DFT(a, lim, -1), DFT(c, lim, -1);
    FOR(i, 1, n)
        printf("%.3lf\n", a[i].x - c[n - i].x);
    return 0;
}

一	二	三	四	五	六	日
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30	31

【算法笔记】多项式合集

多项式合集

拉格朗日插值

问题背景

结论

推导

拉格朗日插值与范德蒙矩阵

快速傅里叶变换（FFT)

复数基础

从欧拉公式到单位圆

多项式的表示法

单位复数根

DFT

IDFT

位逆序置换

例题

作者

留言

撰写回覆或留言取消回复

【算法笔记】多项式合集

多项式合集

拉格朗日插值

问题背景

结论

推导

拉格朗日插值与范德蒙矩阵

快速傅里叶变换（FFT)

复数基础

从欧拉公式到单位圆

多项式的表示法

单位复数根

DFT

IDFT

位逆序置换

例题

作者

留言

撰写回覆或留言 取消回复

撰写回覆或留言取消回复