多项式

多项式

多项式基础

数域的定义

定义 复数集的子集 \(K\) ,满足 \(0,1 \in K\) 且元素的和差积商(除数不为 \(0\) )都属于 \(K\) ,那么称 \(K\) 是一个数域。

关于群环域的详细定义可以看看抽代,这里只提及必须知道的。

例如有理数集、复数集、模素数 \(m\) 剩余系都是数域,但整数集不是数域。

多项式的定义与基本性质

定义\(a_0, a_1, \cdots ,a_n\) 是数域 \(K\) 上的数,其中 \(n\) 为非负整数,那么称 \(f(x) = \displaystyle \sum_{i = 0}^n a_i x^i\) 是数域 \(K\) 上的一元多项式,其中 \(a_ix^i\)\(f(x)\)\(i\) 次项,\(a_i\) 则是 \(f(x)\)\(i\) 次项系数。

此外,也可以用 \([x^i]f(x)\) 表示多项式 \(f(x)\)\(i\) 次项系数。

注意,这里的 \(x\) 是形式上的记号,而非真正的数。换句话说,单独写出系数序列也能代表一个多项式, \(x\) 的次数更多时候只是用来区分系数。

一元多项式环的定义 数域 \(K\) 上多项式的全体,称为一元多项式环,记作 \(K[x]\) ,而 \(K\) 称为 \(K[x]\) 的系数域。

次数的定义 对于多项式 \(f(x)\) ,其系数非零的最高次项的次数称为多项式的次数,记作 \(\partial(f(x))\)

相等的定义 若两个多项式对应次数的各系数均相等,那么这两个多项式相等。

零多项式的定义 系数全为 \(0\) 的多项式,记为 \(0\) ,其次数不考虑。

多项式加法的定义 对于两个多项式 \(\displaystyle f(x) = \sum_{i = 0}^n a_i x^i, g(x) = \sum_{i = 0}^m b_i x^i\) ,则 \(\displaystyle f(x) + g(x) = \sum_{i=0}^{\max(n,m)} (a_i + b_i)x^i\)

多项式乘法的定义 对于两个多项式 \(\displaystyle f(x) = \sum_{i = 0}^n a_i x^i, g(x) = \sum_{i = 0}^m b_i x^i\) ,则 \(\displaystyle f(x) \cdot g(x) = \sum_{i=0}^n \sum_{j=0}^m a_ib_jx^{i+j}\)

多项式乘法有一个等价形式,\(\displaystyle f(x) \cdot g(x) = \sum_{k=0}^{n+m} x^{k} \sum_{i = 0}^k a_ib_{k-i}\) ,即求两个多项式系数的加法卷积(下标之和相等的位置的值的乘积之和),这是今后许多问题可以转化为多项式计算的关键。

多项式复合的定义 对于两个多项式 \(\displaystyle f(x) = \sum_{i = 0}^n a_i x^i, g(x) = \sum_{i = 0}^m b_i x^i\) ,则 \(\displaystyle f(g(x)) = a_0 + \sum_{i=1}^n a_ig^i(x)\)

性质1 数域 \(K\) 上的两个多项式经过加、减、乘等运算后,所得结果仍然是数域 \(K\) 上的多项式。

性质2 对于两个多项式 \(f(x),g(x)\) ,满足加法结合律、加法交换律、乘法结合律、乘法交换律、乘法对加法分配律、乘法消去律。

性质3 任意 \(n+1\) 个不同的采样点,就可以唯一确定一个 \(n\) 次多项式。

多项式带余式除法

定理1(带余式除法) 在一元多项式环 \(K[x]\) 中,任意两个多项式 \(A(x),B(x)\)\(B(x) \neq 0\) ,一定存在唯一的两个多项式 \(Q(x),R(x)\) 满足 \(A(x) = Q(x)B(x) + R(x)\) ,其中 \(\partial(R(x)) < \partial(B(x))\)\(R(x) = 0\) ,称 \(Q(x)\)\(A(x)\) 除以 \(B(x)\) 的商, \(R(x)\)\(A(x)\) 除以 \(B(x)\) 的余式。

大部分数论整除同余的性质都可以类似地应用到多项式上,后面就不展开了。

形式幂级数的定义

定义\(a_0, a_1, \cdots ,a_n\) 是数域 \(K\) 上的数,那么称 \(f(x) = \displaystyle \sum_{i = 0}^{\infin} a_i x^i\) 是数域 \(K\) 上的形式幂级数,简称幂级数。

形式幂级数环的定义 数域 \(K\) 上形式幂级数的全体,称为形式幂级数环,记作 \(K[[x]]\)

幂级数可以看作是一元多项式的扩展,其具有更多良好的性质,如形式导数和形式不定积分等。

在模意义下,幂级数可等价为有限项的多项式,因此通常会把多项式扩展到幂级数上,借由幂级数的诸多性质得到许多有用的结论,例如模意义下多项式的初等函数。

幂级数的导数和不定积分

注意,极限在环上可能并不存在,但依然可以在形式上的定义导数与积分。

形式导数 设形式幂级数 \(\displaystyle f(x) = \sum_{i = 0}^{\infin}a_ix^i\) ,其形式导数 \(\displaystyle f'(x) = \sum_{i = 1}^{\infin}ia_ix^{i-1}\)

此外,我们也可将 \(f(x)\)\(t\) 阶导数记作 \(f^{(t)}(x)\)

形式不定积分 设形式幂级数 \(\displaystyle f(x) = \sum_{i = 0}^{\infin}a_ix^i\) ,其形式不定积分 \(\displaystyle \int f(x) \text{d} x = \sum_{i = 1}^{\infin}ia_ix^{i-1} + C\)

其他的基本求导法则皆可适用,就不再展开了。

常见幂级数展开

\[\begin{aligned} e^x &= \sum_{i = 0}^{\infin} \frac{1}{i!}x^i \\ \sin x &= \sum_{i = 0}^{\infin} \frac{(-1)^{i}}{(2i+1)!}x^{2i+1} \\ \cos x &= \sum_{i = 0}^{\infin} \frac{ (-1)^{i}}{(2i)!}x^{2i} \\ \frac{1}{1+x} &= \sum_{i = 0}^{\infin} (-1)^ix^i \\ (1+x)^{\alpha} &= \sum_{i = 0}^{\infin} \frac{\alpha^{\underline i}}{i!}x^i \\ \ln(1+x) &= \sum_{i = 1}^{\infin} \frac{(-1)^{i-1}}{i}x^i \\ \arctan x &= \sum_{i = 0}^{\infin} \frac{(-1)^{i}}{2i+1}x^{2i+1} \\ \end{aligned} \]

多项式插值

多项式插值的定义

定义 给定 \(n\) 个点 \((x_1,y_1), \cdots, (x_n, y_n)\) ,其中 \(x_i\) 互不相同,求这些点确定的 \(n-1\) 次多项式函数 \(f(x)\) 的过程,称为多项式插值。

多项式插值的方法

拉格朗日插值法

考虑构造 \(n\) 个辅助函数 \(\displaystyle g_i(x) = y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}\) ,显然 \(g_i(x)\) 满足 \(g_i(x_i) = y_i\)\(g_i(x_j) = 0, j \neq i\)

因此我们令 \(\displaystyle f(x) = \sum_{i = 1}^n g_i(x) = \sum_{i = 1}^n y_i \prod_{j \neq i} \frac{x-x_j}{x_i-x_j}\) 即可唯一确定所求多项式,此为拉格朗日插值公式。

其中,若 \(x_i = i\) ,可以预处理阶乘以及 \(x - x_j\) 的前后缀积,将公式化简为 \(O(n)\)

若我们只需要求出在 \(x = k\) 处的值,那么代入即可。

若要求出具体的系数则设计多项式运算的模拟。

这里只给出求单点值的代码。

时间复杂度 \(O(n^2)\)

空间复杂度 \(O(n)\)

int lagrange_poly(const vector<pair<int, int>> &point, int x) {
    int n = point.size() - 1;
    int ans = 0;
    for (int i = 1;i <= n;i++) {
        int res1 = point[i].second;
        int res2 = 1;
        for (int j = 1;j <= n;j++) {
            if (i == j) continue;
            res1 = 1LL * res1 * (x - point[j].first + P) % P;
            res2 = 1LL * res2 * (point[i].first - point[j].first + P) % P;
        }
        (ans += 1LL * res1 * qpow(res2, P - 2) % P) %= P;
    }
    return ans;
}

重心拉格朗日插值法

若插值点会新增 \(q\) 次,每次重新计算 \(f(k)\) 都是 \(O(n^2)\) ,我们需要对公式做些调整。

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

我们设 \(\displaystyle A = \prod_{i=1}^n (x - x_i) , B(i)=\displaystyle \prod_{j \neq i} (x_i-x_j)\)

每次新增插值点时 \(O(1)\) 更新 \(A\)\(O(n)\) 更新 \(B(i)\) 后,即可在 \(O(n)\) 内得到新的 \(\displaystyle f(k) = A \sum_{i = 1}^n \frac{y_i}{(k-x_i)B(i)}\)

代码可借鉴的拉格朗日插值法做修改,这里就不给出了。

时间复杂度 \(O(n^2 + qn)\)

空间复杂度 \(O(n)\)

加法卷积

加法卷积的定义

定义 对于两个序列 \(f,g\) ,它们的加法卷积序列 \(h\) 满足 \(\displaystyle h_k = \sum_{i + j = k} f_ig_j\)

把序列当作多项式系数理解, \(h\) 其实就是 \(f,g\) 表示的多项式的乘积的系数,因此可以用加法卷积的算法加速多项式乘积,下面也会用多项式的角度介绍加法卷积的算法。

加法卷积的变换

快速傅里叶变换(FFT)

多项式在系数域直接加法卷积是 \(O(n^2)\) 的,但如果我们在若干个不同位置取两个多项式的点值(采样点),容易发现这些点值相乘后得到的新点值就落在所求的多项式上,最后只要把点值变换回系数,就得到目标多项式。

换句话说,系数域的加法卷积可以变换为点值域的对应相乘,而对应相乘这个过程是 \(O(n)\) 的,现在我们只需要一个能够快速在系数域和点值域之间变换算法即可。

这也是大多数变换加速卷积的核心思路,即找到一个快速的可逆变换,使得两个序列变换后的对应位置做运算的结果,恰好为两个序列卷积的变换,最后逆变换回去就是两个序列的卷积。

接下来就是加法卷积的需要的变换,离散傅里叶变换。

离散傅里叶变换(DFT)

首先,点值域的点不能随便取,我们要取 \(n\) 次单位根 \(\omega_n\)\(n\) 个幂 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\)\(n\) 要大于等于多项式的项数。为了方便,我们通常需要令 \(n\)\(2\) 的幂。

\(n\) 次单位根等价于将复平面单位圆弧划分为 \(n\) 等分,其中第 \(k\) 份即 \(\omega_n^k = \cos \dfrac{2k\pi}{n} + \text{i}\sin \dfrac{2k\pi}{n}\)

单位根具有许多有用的性质:

  1. 互异性:若 \(i \neq j\) ,则 \(\omega_n^i \neq \omega_n^j\)
  2. 折半律:\(\omega_{2n}^{2i} = \omega_{n}^{i}\)
  3. 周期律:\(\omega_n^{i + n} = \omega_n^i\)
  4. 半周期律: \(\omega_n^{i + \frac{n}{2}} = -\omega_n^i\)

互异性保证了 \(n\) 个点一定互不相同,接下来考虑如何求值。

\(\displaystyle f(x) = \sum_{i = 0}^{n-1} a_i x^i\) ,那么显然有 \(f(x) = a_0 + a_2x^2 + \cdots + a_{n-1}x^{n-1} + x(a_1 + a_3x^2 + \cdots + a_nx^n)\)

\(f_1(x) = a_0 + a_2x + \cdots a_{n-1}x^{n-1} ,f_2(x) = a_1 + a_3x + \cdots + a_nx^n\) ,那么有 \(f(x) = f_1(x^2) + xf_2(x^2)\)

\(i \in [0, \dfrac{n}{2} - 1]\) 时,我们代入单位根 \(\omega_n^i\) 可得

\[\begin{aligned} f(\omega_n^i) &= f_1((\omega_n^i)^2) +\omega_n^if_2((\omega_n^i)^2) \\ &= f_1(\omega_n^{2i}) +\omega_n^if_2(\omega_n^{2i}) \\ &= f_1(\omega_{\frac{n}{2}}^{i}) +\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ \end{aligned} \]

我们代入单位根 \(\omega_n^{i + \frac{n}{2}}\) 可得

\[\begin{aligned} f(\omega_n^{i + \frac{n}{2}}) &= f_1((\omega_n^{i + \frac{n}{2}})^2) +\omega_n^{i + \frac{n}{2}}f_2((\omega_n^{i + \frac{n}{2}})^2) \\ &= f_1(\omega_n^{2i}) -\omega_n^if_2(\omega_n^{2i}) \\ &= f_1(\omega_{\frac{n}{2}}^{i}) -\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ \end{aligned} \]

注意到 \(f_1(\omega_{\frac{n}{2}}^{i})\)\(f_2(\omega_{\frac{n}{2}}^{i})\) 正是子问题的答案。

因此一个大小为 \(n\) 的问题,变成两个大小为 \(\dfrac{n}{2}\) 的子问题外加 \(O(n)\) 复杂度计算,递归下去总体复杂度是 \(O(n\log n)\) 的。(如果随便取点,问题规模不会折半,也就不能快速了)

于是,多项式系数域到点值域的快速正变换就找到了。

位逆序置换

递归的常数较大,我们希望改为迭代,考虑 \(2^3\) 项多项式的变换过程:

  1. 初始层:\(\{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\}\)
  2. 第一层:\(\{a_0, a_2, a_4, a_6\},\{a_1, a_3, a_5, a_7\}\)
  3. 第二层:\(\{a_0, a_4\},\{a_2, a_6\},\{a_1, a_5\}, \{a_3, a_7\}\)
  4. 第三层:\(\{a_0\}, \{a_4\},\{a_2\}, \{a_6\},\{a_1\}, \{a_5\}, \{a_3\}, \{a_7\}\)

我们需要从下往上迭代,那么就需要知道最后一层的顺序。

容易知道,\(a_i\) 最后会出现在 \(a_{rev_i}\) ,其中 \(rev_i\) 表示 \(i\) 的二进制逆序,例如 \(110\) 的逆序就是 \(011\)

根据 \(rev\) 的定义,我们可以递推它在 \(n\) 项多项式的情况:

\[\begin{aligned} rev_i = \dfrac{rev_{\frac{i}{2}}}{2} + [2 \not\mid i] \cdot 2^{\log n - 1} \end{aligned} \]

因此,我们预处理 \(rev\) 后,将对应系数置换即可迭代。

蝶形运算优化

在上面,当我们求出 \(f_1(x)\)\(f_2(x)\) 的各 \(\dfrac{n}{2}\) 个点值后,设 \(i \in [0, \dfrac{n}{2}-1]\) ,那么

\[\begin{aligned} f(\omega_n^i) &= f_1(\omega_{\frac{n}{2}}^{i}) +\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ f(\omega_n^{i + \frac{n}{2}}) &= f_1(\omega_{\frac{n}{2}}^{i}) -\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ \end{aligned} \]

注意到 \(f(\omega_n^i)\) 和 $ f(\omega_n^{i + \frac{n}{2}})$ 分别在 \(i\)\(i + \frac{n}{2}\) 位置上,而 \(f_1(\omega_{\frac{n}{2}}^{i})\)\(f_2(\omega_{\frac{n}{2}}^{i})\) 也恰好在 \(i\)\(i + \frac{n}{2}\) 位置上,因此我们不需要额外空间,只需要原地覆盖即可。

离散傅里叶逆变换(IDFT)

现在,我们尝试推导一下逆变换。

我们定义点值多项式 \(\displaystyle F(x) = \sum_{i = 0}^{n-1} f(\omega_n^{i})x^i\) ,即 \(f(x)\) 点值当作系数的多项式。

我们代入 \(x = \omega_n^k\) ,那么 \(\displaystyle F(\omega_n^k) = \sum_{i = 0}^{n-1} \omega_n^{ik}\sum_{j = 0}^{n-1} a_j\omega_n^{ij} = \sum_{j = 0}^{n-1} a_j\sum_{i = 0}^{n-1} \omega_n^{i(k+j)}\)

构造辅助多项式 \(\displaystyle G(x) = \sum_{i = 0}^{n-1} x^i\) ,那么 \(\displaystyle F(\omega_n^k) = \sum_{j=0}^{n-1}a_jG(\omega_n^{j+k})\)

考虑 \(G(\omega_n^i)\) ,当 \(i = 0\)\(G(\omega_n^i) = n\) ,否则单位根两两配对 \(G(\omega_n^i) = 0\)

因此 \(F(\omega_n^k) = na_{n-k}\) ,特别地当 \(k = 0\)\(F(\omega_n^0) = na_0\) ,所以我们只需要对点值多项式进行一次DFT,随后将 \([1,n-1]\) 项反转,最后对所有项除以 \(n\) 即可还原多项式。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

const db PI = acos(-1.0);

vector<int> rev;
vector<complex<db>> Wn = { {0, 0}, {1, 0} }; // 0, w20, w40, w41, w80, w81, w82, w83, ...
void FFT(vector<complex<db>> &A, bool is_inv) {
    int n = A.size();

    if (rev.size() != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0;i < n;i++) rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
    }
    for (int i = 0;i < n;i++) if (i < rev[i]) swap(A[i], A[rev[i]]);

    if (Wn.size() < n) {
        int k = Wn.size();
        Wn.resize(n);
        while (k < n) {
            complex<db> w(cos(PI / k), sin(PI / k));
            for (int i = k >> 1;i < k;i++) {
                Wn[i << 1] = Wn[i];
                Wn[i << 1 | 1] = Wn[i] * w;
            }
            k <<= 1;
        }
    }

    for (int k = 1;k < n; k <<= 1) {
        for (int i = 0;i < n;i += k << 1) {
            for (int j = 0;j < k;j++) {
                complex<db> u = A[i + j];
                complex<db> v = A[i + k + j] * Wn[k + j];
                A[i + j] = u + v;
                A[i + k + j] = u - v;
            }
        }
    }

    if (is_inv) {
        reverse(A.begin() + 1, A.end());
        for (int i = 0;i < n;i++) A[i] /= n;
    }
}

vector<complex<db>> poly_mul(vector<complex<db>> A, vector<complex<db>> B) {
    if (A.empty() || B.empty()) return vector<complex<db>>();
    int n = 1, sz = A.size() + B.size() - 1;
    while (n < sz) n <<= 1;
    A.resize(n);
    B.resize(n);
    FFT(A, 0);
    FFT(B, 0);
    for (int i = 0;i < n;i++) A[i] *= B[i];
    FFT(A, 1);
    A.resize(sz);
    return A;
}

快速数论变换(NTT)

考虑在素域内的多项式变换,我们发现原根刚好能代替单位根。

考虑一个素数 \(P\) ,表达为 \(P = r \cdot 2^k + 1\) ,其原根 \(G\) 的阶为 \(\varphi(P) = P-1 = r \cdot 2^k\) ,当多项式含有 \(n = 2^{k'}\) 项时,我们令 \(G_n = G^{\frac{P-1}{n}}\) 那么 \(G_n\) 等价为 \(n\) 次单位根。

注意到,一个素数能支持的多项式长度为 \(2^k\) ,因此 \(k\) 越大越好,不过常见的 \(10^9 + 7\) 就比较鸡肋,因为 \(k = 1\)

NTT其他部分和FTT完全一致。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

const int P = 998244353, G = 3;

int qpow(int a, ll k) {
    int ans = 1;
    while (k) {
        if (k & 1) ans = 1LL * ans * a % P;
        k >>= 1;
        a = 1LL * a * a % P;
    }
    return ans;
}

std::vector<int> rev;
std::vector<int> Wn = { 0,1 }; // 0, w20, w40, w41, w80, w81, w82, w83, ...
/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
    explicit Poly(int n = 0, int val = 0) : std::vector<int>(n, val) {}
    explicit Poly(const std::vector<int> &src) : std::vector<int>(src) {}
    explicit Poly(const std::initializer_list<int> &src) : std::vector<int>(src) {}

    Poly modx(int k) const {
        assert(k >= 0);
        if (k > size()) {
            Poly X = *this;
            X.resize(k);
            return X;
        }
        else return Poly(std::vector<int>(begin(), begin() + k));
    }
    Poly mulx(int k) const {
        assert(k >= 0 || -k < size());
        Poly X = *this;
        if (k >= 0) X.insert(X.begin(), k, 0);
        else X.erase(X.begin(), X.begin() + (-k));
        return X;

    }
    Poly derive(int k = 0) const {
        if (k == 0) k = std::max((int)size() - 1, 0);
        Poly X(k);
        for (int i = 1;i < std::min((int)size(), k + 1);i++) X[i - 1] = 1LL * i * (*this)[i] % P;
        return X;
    }
    Poly integral(int k = 0) const {
        if (k == 0) k = size() + 1;
        Poly X(k);
        for (int i = 0;i < std::min((int)size(), k - 1);i++)  X[i + 1] = 1LL * qpow(i + 1, P - 2) * (*this)[i] % P;
        return X;
    }
    
    Poly &operator+=(const Poly &X) {
        resize(std::max(size(), X.size()));
        for (int i = 0;i < size();i++) ((*this)[i] += X[i]) %= P;
        return *this;
    }
    Poly &operator-=(const Poly &X) {
        resize(std::max(size(), X.size()));
        for (int i = 0;i < size();i++) ((*this)[i] -= X[i] - P) %= P;
        return *this;
    }
    Poly &operator*=(int x) {
        for (int i = 0;i < size();i++) (*this)[i] = 1LL * (*this)[i] * x % P;
        return *this;
    }
    Poly &operator/=(int x) {
        int val = qpow(x, P - 2);
        for (int i = 0;i < size();i++) (*this)[i] = 1LL * (*this)[i] * val % P;
        return *this;
    }
    friend Poly operator+(Poly A, const Poly &B) { return A += B; }
    friend Poly operator-(Poly A, const Poly &B) { return A -= B; }
    friend Poly operator*(Poly A, int x) { return A *= x; }
    friend Poly operator*(int x, Poly A) { return A *= x; }
    friend Poly operator/(Poly A, int x) { return A /= x; }
    friend Poly operator-(Poly A) { return (P - 1) * A; }
    friend std::ostream &operator<<(std::ostream &os, const Poly &X) {
        for (int i = 0;i < X.size();i++) os << X[i] << ' ';
        return os;
    }

    static void NTT(Poly &X, bool is_inv) {
        int n = X.size();

        if (rev.size() != n) {
            int k = __builtin_ctz(n) - 1;
            rev.resize(n);
            for (int i = 0;i < n;i++) rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
        for (int i = 0;i < n;i++) if (i < rev[i]) std::swap(X[i], X[rev[i]]);

        if (Wn.size() < n) {
            int k = __builtin_ctz(Wn.size());
            Wn.resize(n);
            while (1 << k < n) {
                int w = qpow(G, P - 1 >> k + 1);
                for (int i = 1 << k - 1;i < 1 << k;i++) {
                    Wn[i << 1] = Wn[i];
                    Wn[i << 1 | 1] = 1LL * Wn[i] * w % P;
                }
                k++;
            }
        }

        for (int k = 1;k < n; k <<= 1) {
            for (int i = 0;i < n;i += k << 1) {
                for (int j = 0;j < k;j++) {
                    int u = X[i + j];
                    int v = 1LL * X[i + k + j] * Wn[k + j] % P;
                    X[i + j] = (u + v) % P;
                    X[i + k + j] = (u - v + P) % P;
                }
            }
        }

        if (is_inv) {
            std::reverse(X.begin() + 1, X.end());
            int inv = qpow(n, P - 2);
            for (int i = 0;i < n;i++) X[i] = 1LL * X[i] * inv % P;
        }
    }
    Poly &operator*=(Poly X) {
        if (empty() || X.empty()) return *this = Poly();
        int n = 1, sz = size() + X.size() - 1;
        while (n < sz) n <<= 1;
        resize(n);
        X.resize(n);
        NTT(*this, 0);
        NTT(X, 0);
        for (int i = 0;i < n;i++) (*this)[i] = 1LL * (*this)[i] * X[i] % P;
        NTT(*this, 1);
        resize(sz);
        return *this;
    }
    friend Poly operator*(Poly A, const Poly &B) { return A *= B; }
};

常用NTT模数:

\(r\cdot 2^k + 1\) \(r\) \(k\) \(g\)
5767169 11 19 3
7340033 7 20 3
23068673 11 21 3
104857601 25 22 3
167772161 5 25 3
469762049 7 26 3
998244353 119 23 3
1004535809 479 21 3
2013265921 15 27 31
2281701377 17 27 3
3221225473 3 30 5
75161927681 35 31 3
77309411329 9 33 7
206158430209 3 36 22
2061584302081 15 37 7
2748779069441 5 39 3
6597069766657 3 41 5
39582418599937 9 42 5
79164837199873 9 43 5
263882790666241 15 44 7
1231453023109121 35 45 3
1337006139375617 19 46 3
3799912185593857 27 47 5
4222124650659841 15 48 19
7881299347898369 7 50 6
31525197391593473 7 52 3
180143985094819841 5 55 6
1945555039024054273 27 56 5
4179340454199820289 29 57 3

任意模数NTT(MTT)

暂时不学。

多项式初等函数

初等函数 公式 方法 备注
乘法 \(f(x) \cdot g(x)\) NTT/FTT
求逆 \(f^{-1}(x) \equiv f^{-1}_0(x)(2 - f^{-1}_0(x)f(x)) \pmod{x^n}\) 牛顿迭代、乘法 常数项逆元存在
整除 \(\left[\dfrac{f(x)}{g(x)} \right]_R \equiv f_R(x)g_R^{-1}(x) \pmod{x^{n-m+1}}\) 求逆 除式非零
取模 \(f(x) \bmod g(x) = f(x) - g(x)\left[\dfrac{f(x)}{g(x)} \right]\) 整除 除式非零
开方 \(\sqrt{f(x)} \equiv \dfrac{1}{2} \left(\left( \sqrt{f(x)} \right)_0 + f(x)\left( \sqrt{f(x)} \right)_0^{-1} \right) \pmod{x^n}\) 牛顿迭代、求逆 首非零项是偶次项,且二次剩余存在
对数函数 \(\displaystyle \ln f(x) \equiv \int f'(x)f^{-1}(x) \text{d}x \pmod{x^n}\) 求逆 常数项为 \(1\)
指数函数 \(\text{e}^{f(x)} \equiv \left(\text{e}^{f(x)}\right)_0 \left(1-\ln \left(\text{e}^{f(x)}\right)_0 + f(x) \right) \pmod{x^n}\) 牛顿迭代、对数函数 常数项为 \(0\)
幂函数 \(f^k(x) \equiv e^{k \ln f(x)} \pmod{x^n}\) 指数函数
三角函数 欧拉公式 指数函数 常数项为 \(0\)
反三角函数 求导积分 开方 常数项为 \(0\)

多项式牛顿迭代

给定多项式 \(g(x)\) ,求 \(f(x)\) ,满足 \(g(f(x)) \equiv 0 \pmod{x^n}\)

考虑倍增法。

\(n = 1\) 时, \([x^0]g(f(x)) = 0\) 需要单独解出。

假设在模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 时的 \(f(x)\) 的解是 \(f_0(x)\) ,那么我们对 \(g(f(x))\)\(f_0(x)\) 处泰勒展开有

\[\displaystyle \sum_{i=0}^{\infin} \dfrac{g^{(i)}(f_0(x))}{i!}(f(x) - f_0(x))^i \equiv 0 \pmod {x^n} \]

显然,当 \(i \geq 2\) 时, \((f(x) - f_0(x))^i \equiv 0 \pmod{x^n}\) ,因此有

\[\displaystyle \sum_{i=0}^{\infin} \dfrac{g^{(i)}(f_0(x))}{i!}(f(x) - f_0(x))^i \equiv g(f_0(x)) + g'(f_0(x))(f(x) - f_0(x))) \equiv 0 \pmod {x^n} \]

最后化简得

\[f(x) \equiv f_0(x) - \frac{g(f_0(x))}{g'(f_0(x))} \pmod{x^n} \]

这就是模意义下的牛顿迭代,每次都能倍增多项式有效系数,一些关键多项式初等函数需要由此推导。

\(x^n\) 是因为精确解实际上大概率是幂级数,但大部分时候我们的操作只需要前几项,就能保证覆盖所有有意义的部分,因此幂级数是不必要的,求出模意义下的就够了。

多项式求逆

给定多项式 \(f(x)\) ,求 \(f^{-1}(x)\) ,满足 \(f(x)f^{-1}(x) \equiv 1 \pmod{x^n}\)

\(g(f^{-1}(x)) = \dfrac{1}{f^{-1}(x)} - f(x) \equiv 0 \pmod{x^n}\)

\(n = 1\) 时,\([x^0]f^{-1}(x) = ([x^0]f(x))^{-1}\) ,因此需要保证常数项逆元存在。

假设模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 的解为 \(f_0(x)\)

根据牛顿迭代

\[f^{-1}(x) \equiv f^{-1}_0(x) - \dfrac{g(f^{-1}_0(x))}{g'(f^{-1}_0(x))} \equiv f^{-1}_0(x) - \dfrac{\dfrac{1}{f^{-1}_0(x)} - f(x)}{-\dfrac{1}{f^{-2}_0(x)}} \equiv f^{-1}_0(x)(2 - f^{-1}_0(x)f(x)) \pmod{x^n} \]

因此,我们可以用这个公式迭代出多项式的逆,复杂度由主定理 \(T(n) = T(n/2) + O(n \log n) = O(n \log n)\)

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
	Poly inv(int n = 0) const {
        assert(size() && (*this)[0] > 0); // assert [x^0]f(x) inverse exists
        if (n == 0) n = size();
        Poly X{ qpow((*this)[0], P - 2) };
        int k = 1;
        while (k < n) {
            k <<= 1;
            X = (X * (Poly{ 2 } - X * modx(k))).modx(k);
        }
        return X.modx(n);
    }
};

多项式除法&取模

给定多项式 \(f(x),g(x)\) ,求 \(q(x),r(x)\) ,满足 \(f(x) = q(x)g(x) + r(x)\)

其中 \(\partial(q(x)) = \partial(f(x)) - \partial(g(x)), \partial(r(x)) < \partial(g(x))\)

\(n = \partial(f(x)), m = \partial(g(x))\) ,不妨设 \(\partial(r(x)) = m-1\)

因为存在余式,我们不能直接使用模 \(x^m\) 的多项式求逆。

\(f_R(x) = x^nf\left( \dfrac{1}{x} \right)\) ,其实就是将系数反转。

我们对原式变形

\[d\begin{aligned} f(x) &= q(x)g(x) + r(x)\\ f\left( \dfrac{1}{x} \right) &= q\left( \dfrac{1}{x} \right)g\left( \dfrac{1}{x} \right) + r\left( \dfrac{1}{x} \right) \\ x^nf\left( \dfrac{1}{x} \right) &= x^nq\left( \dfrac{1}{x} \right)g\left( \dfrac{1}{x} \right) + x^nr\left( \dfrac{1}{x} \right) \\ f_R(x) &= g_R(x)q_R(x) + x^{n - m + 1} r_R(x) \end{aligned} \]

\(\partial(q_R(x)) = \partial(q(x)) = n-m < n-m+1\) ,因此在模 \(x^{n-m+1}\)\(q_R(x)\) 是不会被影响的,而 \(x^{n-m+1}r_R(x)\) 会被模掉。

所以有 \(f_R(x) \cdot g^{-1}_R(x) \equiv q_R(x) \pmod{x^{n-m+1}}\)

求出 \(q_R(x)\) 后,反转系数就是 \(q(x)\) ,最后 \(r(x) = f(x) - q(x)g(x)\)

实现上注意处理除式的后导 \(0\) ,会导致结果出错,虽然一般题目不需要这个处理。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
    Poly &operator/=(Poly X) {
        while (X.size() && X.back() == 0) X.pop_back();
        assert(X.size()); // assert X(x) is not 0-polynomial
        if (size() < X.size()) return *this = Poly();
        std::reverse(begin(), end());
        std::reverse(X.begin(), X.end());
        *this = (modx(size() - X.size() + 1) * X.inv(size() - X.size() + 1)).modx(size() - X.size() + 1);
        std::reverse(begin(), end());
        return *this;
    }
    Poly &operator%=(Poly X) {
        while (X.size() && X.back() == 0) X.pop_back();
        return *this = (*this - *this / X * X).modx(X.size() - 1);
    }
    friend Poly operator/(Poly A, const Poly &B) { return A /= B; }
    friend Poly operator%(Poly A, const Poly &B) { return A %= B; }
};

多项式开方

给定多项式 \(f(x)\) ,求 \(\sqrt{f(x)} \bmod x^n\)

\(g(\sqrt{f(x)}) = \left( \sqrt{f(x)} \right)^2 - f(x) \equiv 0 \pmod {x^n}\)

\(n = 1\) 时, \([x^0]\sqrt{f(x)} = \sqrt{[x^0]f(x)}\) ,因此需要保证常数项二次剩余存在。

假设模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 的解为 \(f_0(x)\)

根据牛顿迭代

\[\sqrt{f(x)} \equiv \left( \sqrt{f(x)} \right)_0 - \frac{\left( \sqrt{f(x)} \right)_0^2 - f(x)}{2\left( \sqrt{f(x)} \right)_0} \equiv \dfrac{1}{2} \left(\left( \sqrt{f(x)} \right)_0 + f(x)\left( \sqrt{f(x)} \right)_0^{-1} \right) \pmod{x^n} \]

代码没实现求二次剩余,目前只能对常数项为 \(1\) 的开方。

特别地,出现前导 \(0\) 时,前导 \(0\) 个数 \(cnt\) 是偶数(即第一个非零项是偶次)则多项式可以整体除以 \(x^{cnt}\) 再开方,最后乘 \(x^{cnt/2}\) ,否则无解。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

struct Poly : public std::vector<int> {
    Poly sqrt(int n = 0) const {
        if (n == 0) n = size();
        int cnt = 0;
        while (cnt < size() && (*this)[cnt] == 0) cnt++;
        if (cnt == size()) return Poly(n);
        assert(!(cnt & 1) && (*this)[cnt] == 1); // assert cnt is even and [x^cnt]f(x) exists 2-residue 
        Poly X{ 1 };
        int k = 1;
        while (k < n) {
            k <<= 1;
            X = (P + 1 >> 1) * (X + mulx(-cnt).modx(k) * X.inv(k)).modx(k);
        }
        return X.mulx(cnt >> 1).modx(n);
    }
};

多项式对数函数

给定多项式 \(f(x)\) ,求 \(\ln f(x) \bmod x^n\)

求导再积分后, \(\displaystyle \ln f(x) \equiv \int \frac{f'(x)}{f(x)} \text{d}x \pmod {x^n}\)

注意根据 \(\ln\) 的定义, \(f(x)\) 的常数项必须为 \(1\) ,否则对数无意义。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
	Poly log(int n = 0) const {
        assert(size() && (*this)[0] == 1); // assert [x^0]f(x) = 1
        if (n == 0) n = size();
        return (derive(n) * inv(n)).integral(n);
    }
};

多项式指数函数

给定多项式 \(f(x)\) ,求 \(e^{f(x)} \bmod x^n\)

\(g(e^{f(x)}) = \ln e^{f(x)} - f(x) \equiv 0 \pmod{x^n}\)

\(n = 1\) 时,\([x^0]e^{f(x)} = e^{[x^0]f(x)}\) ,因此需要保证常数项为 \(0\) ,否则无意义。

假设模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 的解为 \(f_0(x)\)

根据牛顿迭代

\[e^{f(x)} \equiv \left( e^{f(x)} \right)_0 - \frac{\ln \left( e^{f(x)} \right)_0 - f(x)}{\frac{1}{\left( e^{f(x)} \right)_0}} \equiv \left( e^{f(x)} \right)_0 \left( 1 - \ln \left( e^{f(x)} \right)_0 + f(x) \right) \pmod{x^n} \]

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
	Poly exp(int n = 0) const {
        assert(empty() || (*this)[0] == 0); // assert [x^0]f(x) = 0
        if (n == 0) n = size();
        Poly X{ 1 };
        int k = 1;
        while (k < n) {
            k <<= 1;
            X = (X * (Poly{ 1 } - X.log(k) + modx(k))).modx(k);
        }
        return X.modx(n);
    }
};

多项式幂函数

给定多项式 \(f(x)\) ,求 \(f^k(x) \bmod x^n\)

显然有 \(f^k(x) \equiv e^{k \ln f(x)} \pmod{x^n}\)

指数并非真正的指数,而是多项式函数的自变量,因此指数上的 \(k\) 也是对 \(P\) 取模。

\([x^0]f(x) \neq 1\) 时,我们找到第一个非零项 \([x^{cnt}]f(x)\) ,多项式可以整体除以 \([x^{cnt}]f(x) \cdot x^{cnt}\) 再用上面的公式,最后乘 \(([x^{cnt}]f(x))^k \cdot x^{k \cdot cnt}\) ,其中 \([x^{cnt}]f(x)\)\(k\) 次方要模 \(\varphi(P)\) ,因为他是真正意义的指数。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
	Poly pow(int k = 0, int n = 0) const {
        if (n == 0) n = size();
        int cnt = 0;
        while (cnt < size() && (*this)[cnt] == 0) cnt++;
        if (cnt == size()) return Poly(n);
        if (1LL * k * cnt >= n) return Poly(n);
        int k1 = k % P, k2 = k % (P - 1);
        return ((k1 * (mulx(-cnt).modx(n) / (*this)[cnt]).log(n)).exp(n) * qpow((*this)[cnt], k2)).mulx(cnt * k1).modx(n);
    }
    Poly pow(const std::string &s, int n = 0) const {
        if (n == 0) n = size();
        int cnt = 0;
        while (cnt < size() && (*this)[cnt] == 0) cnt++;
        if (cnt == size()) return Poly(n);
        int k1 = 0, k2 = 0;
        for (auto ch : s) {
            if ((1LL * 10 * k1 + ch - '0') * cnt >= n) return Poly(n);
            k1 = (1LL * 10 * k1 % P + ch - '0') % P;
            k2 = (1LL * 10 * k2 % (P - 1) + ch - '0') % (P - 1);
        }
        return ((k1 * (mulx(-cnt).modx(n) / (*this)[cnt]).log(n)).exp(n) * qpow((*this)[cnt], k2)).mulx(cnt * k1).modx(n);
    }
};

多项式三角函数

给定多项式 \(f(x)\) ,求 \(\sin f(x) \bmod x^n\)\(\cos f(x) \bmod x^n\)

考虑欧拉公式 \(e^{\text{i}\theta} = \cos \theta + \text{i} \sin \theta\)

因此有 \(\sin \theta = \dfrac{e^{\text{i} \theta} - e^{-\text{i} \theta}}{2 \text{i}} , \cos \theta = \dfrac{e^{\text{i} \theta} + e^{-\text{i} \theta}}{2}\)

\(\theta = f(x)\) 即可,其中 \(\text{i}\) 在模 \(P\) 下等价于 \(g^{\frac{P-1}{4}}\) ,注意 \(f(x)\) 常数项必须为 \(0\) ,否则无意义。

此外,用这两个函数可以推导出其他三角函数(但有些无法在 \(0\) 处展开,且在其他位置展开都存在超越数,就不存在于这个模多项式体系了),这里就不一一列举了。

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
	Poly sin(int n = 0) const {
        if (n == 0) n = size();
        return ((I * modx(n)).exp(n) - (((P - I) % P) * modx(n)).exp(n)).modx(n) * qpow(2LL * I % P, P - 2);
    }
    Poly cos(int n = 0) const {
        if (n == 0) n = size();
        return ((I * modx(n)).exp(n) + (((P - I) % P) * modx(n)).exp(n)).modx(n) * qpow(2, P - 2);
    }
};

多项式反三角函数

给定多项式 \(f(x)\) ,求 \(\arcsin f(x) \bmod x^n\)\(\arctan f(x) \bmod x^n\)

考虑求导再积分回去, \(\displaystyle \arcsin f(x) = \int \frac{f'(x)}{\sqrt{1 - f^2(x)}} \text{d}x, \arctan f(x) = \int \frac{f'(x)}{1 + f^2(x)} \text{d}x\)

为什么没有 \(\arccos\) ?因为他的多项式常数是超越数,在这个体系无意义。上面能求的积分出来的常数是 \(0\)

时间复杂度 \(O(n \log n)\)

空间复杂度 \(O(n)\)

/// 有限域多项式板子(部分)
struct Poly : public std::vector<int> {
	Poly asin(int n = 0) const {
        if (n == 0) n = size();
        return (derive(n) * (Poly{ 1 } - (modx(n) * modx(n))).sqrt(n).inv(n)).integral(n);
    }
    Poly atan(int n = 0) const {
        if (n == 0) n = size();
        return (derive(n) * (Poly{ 1 } + (modx(n) * modx(n))).inv(n)).integral(n);
    }
};

位运算卷积

位运算卷积的定义

快速沃尔什变换(FWT)

异或卷积

或卷积

与卷积

子集卷积

子集卷积的定义

快速莫比乌斯变换(FMT)

并卷积

交卷积

子集卷积

多项式分治

多项式多点求值

多项式快速插值

分治加法卷积

多项式线性齐次递推

多项式平移

posted @ 2024-03-09 17:59  空白菌  阅读(202)  评论(0编辑  收藏  举报