Loading

多项式全家桶

以下若无特殊说明,函数均为多项式函数。系数对 998244353 取模。

多项式乘法

FFT——从入门到入土:6. 快速数论变换

多项式牛顿迭代

设有任意可泰勒展开的函数 \(g\),求满足 \(g(f(x)) \equiv0 \pmod{x^n}\) 的多项式 \(f(x)\),此时可以使用牛顿迭代法求解。

具体的讲,牛顿迭代法是有以下过程:

  1. 首先求出 \(n=1\) 时的 \(f(x)\)

  2. 假设我们已经求出了 \(g[f(x)]\equiv0\pmod {x^n}\) 的解 \(f_0(x)\),则 \(g[f(x)]\equiv0\pmod{x^{2n}}\) 的解为

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

    证:
    我们直接在 \(f(x)=f_0(x)\) 处对 \(g\) 泰勒展开:

    \[g[f(x)] \equiv \sum_{i\ge0}\frac{g^{(i)}[f_0(x)]}{i!}[f(x)-f_0(x)]^i \pmod{x^{2n}} \]

    考虑到 \(f(x) \equiv f_0(x) \pmod{x^n}\),我们有 \(\forall i\ge2,[f(x)-f_0(x)]^i\equiv 0\pmod{x^{2n}}\)

    因此又有:

    \[\begin{aligned} g[f(x)] &\equiv g[f_0(x)] + g'[f_0(x)]\cdot [f(x) - f_0(x)] \equiv0 &\pmod {x^{2n}}\\ f(x) &\equiv f_0(x) - \frac{g[f_0(x)]}{g'[f_0(x)]} &\pmod{x^{2n}} \end{aligned} \]

时间复杂度:\(T(n)=T(\frac{n}{2})+f(n)\),注意递归时 \(n \rightarrow \lceil\frac{n}{2}\rceil\)

关于求导:将 \(h\) 看作一次自变量,\(f\) 看作常数即可,由于待求值为 \(h\) 而不是 \(x\) 因此本质应为求偏导(或形式幂级数)一类的方法,这里不做深究。

多项式求逆

求满足 \(h(x) \cdot f(x) \equiv 1 \pmod{x^n}\)\(h(x)\)

由题意,\(\frac{1}{h(x)} \equiv f(x) \pmod {x^n}\),考虑使用牛顿迭代,构造函数 \(g(h(x)) = \frac{1}{h(x)}-f(x)\),即求解 \(g[h(x)] \equiv 0 \pmod{x^n}\)

  1. \(n=1\) 时,\(h(x) \equiv \text{inv}([x^0]f(x))\)(显而易见的,\([x^0]f(x)=0\) 时不存在逆元)。

  2. \(n \ge 2\) 时,

    \[\begin{aligned} h(x) &\equiv h_0(x)-\frac{g[h_0(x)]}{g'(h_0(x))} &\pmod{x^n}\\ &\equiv h_0(x) -\frac{h_0^{-1}(x)-f(x)}{-h_0^{-2}(x)} &\pmod{x^n}\\ &\equiv 2h_0(x)-h_0^2(x)f(x) &\pmod{x^n} \end{aligned} \]

多项式 ln

求满足 \(h(x) \equiv \ln f(x) \pmod{x^n}\)\(h(x)\)

考虑求导再积分:

\[\begin{aligned} h'(x) &\equiv \ln'f(x) \cdot f'(x) &\pmod{x^n} \\ &\equiv f^{-1}(x) \cdot f'(x) &\pmod{x^n} \\ h(x) &\equiv \int \frac{f'(x)}{f(x)} \mathrm{d}x + C &\pmod{x^n}\\ \end{aligned} \]

\(x=0\),可以求出 \(C=[x^0]\ln f(x)=h(0)=\ln f(0)\)

由于 \(e\) 是超越数,因此 \(\nexists x \in \Q\setminus\{0\},\ln x\in \Q\)

因此 \(f(0)=1,[x^0]h(x)=h(0)=0\)

多项式 exp

求满足 \(h(x) \equiv e^{f(x)} \pmod{x^n}\)\(h(x)\)

由于 \(\ln h(x)=f(x)\),使用牛顿迭代法构造 \(g(h(x))=\ln h(x)-f(x)\)

首先由上面得 \(h(x) \equiv 1 \pmod{x^n}\),考虑若当前已经求出 \(g[h_0(x)]\equiv0\pmod{x^{n/2}}\)

\[\begin{aligned} h(x)&\equiv h_0(x)-\frac{\ln h_0(x)-f(x)}{\ln'h_0(x)} &\pmod{x^n} \\ &\equiv h_0(x)-\frac{\ln h_0(x)-f(x)}{h_0^{-1}(x)}&\pmod{x^n} \\ &\equiv h_0(x)-h_0(x) \cdot \ln h_0(x) - f(x) &\pmod{x^n} \end{aligned} \]

多项式幂函数

求满足 \(h(x) \equiv f^k(x) \pmod{x^n}\)\(h(x)\)

首先不妨假设 \([x^0]f(x) =1\),对于其他取值可以通过左右移位、乘系数等变换来将常数项变成 1。

利用对数函数的性质:

\[h(x) \equiv f^k(x) \equiv \exp\ln f^k(x) \equiv \exp[k\ln f(x)] \pmod{x^n} \]

多项式三角函数

中学阶段学过的:

\[\cos x= \frac{e^{ix}+e^{-ix}}{2},\sin x=\frac{e^{ix}-e^{ix}}{2i} \]

使用多项式 exp 直接代入即可。(注意 \(\Z_{998244353}\)\(i\) 的取值)。

多项式开方

求满足 \(h^k(x) \equiv f(x) \pmod{x^n}\)\(h(x)\)

首先假设你已经通过一些方法求出了 \(n=1\) 时的解(如二次剩余、BSGS)。

之后牛顿迭代,过程读者自证不难。

\[h(x)\equiv\frac{(k-1)h_0^k(x)-f(x)}{kh_0^{k-1}(x)} \pmod{x^n} \]

多项式除法&取余

设有 \(n\) 项式 \(f(x)\)\(m\) 项式 \(g(x)\),求 \(\frac{f(x)}{g(x)}\) 的商 \(q(x)\) 和余式 \(r(x)\)

考虑到若 \(r(x)=0\),则 \(q(x)=f(x)g^{-1}(x)\),又因为 \(r\) 的次数有特殊性,所以尝试将 \(r(x)\) 利用模运算消掉。

定义 \(\text{rev}_k f\) 表示将 \(f\)\(0 \to k\) 系数翻转,忽略 \(k\) 表示全部翻转,则由多项式乘法的定义不难看出 \(A=B \cdot C \Leftrightarrow \text{rev}A=\text{rev}B \cdot \text{rev}C,\text{rev}_k(A+B)=\text{rev}_kA+\text{rev}_kB\)

则我们有

\[\text{rev}f=\text{rev}(gq+r)=\text{rev}g\cdot\text{rev}q+\text{rev}_n r \]

发现 \(\text{rev}q\) 最高次为 \(n-m\),而 \(\text{rev} r\) 非零最低次至少为 \(n-m+1\),因此直接 \(\bmod x^{n-m+1}\),即变成

\[\text{rev}f \equiv \text{rev}g\cdot\text{rev}q \pmod{x^{n-m+1}} \]

其中 \(g\) 取模不受影响,使用多项式求逆即可算出 \(g\)

板子速查

#include <bits/stdc++.h>
using namespace std;
#define int long long

const int N = 1e6 + 5, MOD = 998244353, G = 3, Gi = 332748118;

namespace Cipolla { // 二次剩余
int n, p, oyds;
struct comp {
    int r, i;
    comp(int _r = 0, int _i = 0) { r = _r, i = _i; }
    comp operator*(const comp &rhs) const {
        return comp((r * rhs.r % p + i * rhs.i % p * oyds) % p,
                    (r * rhs.i % p + i * rhs.r) % p);
    }
};
int qpow(int a, int b) {
    int res = 1;
    for (; b; b >>= 1, a = a * a % p)
        if (b & 1) res = res * a % p;
    return res;
}
comp cpow(comp a, int b) {
    comp res = comp(1, 0);
    for (; b; b >>= 1, a = a * a)
        if (b & 1) res = res * a;
    return res;
}
int check(int x) {
    if (x % p == 0) return 0;
    return qpow(x, (p - 1) >> 1);
}
int work() {
    if (check(n) == p - 1) return -1;
    int a;
    do {
        a = rand() % p;
        oyds = (a * a % p - n + p) % p;
    } while (check(oyds) != p - 1);
    return cpow(comp(a, 1), p + 1 >> 1).r % p;
}
int Cip(int _n, int _p) {
    n = _n, p = _p;
    n %= p;
    if (n == 0) {
        return 0;
    }
    int x1 = (work() + p) % p, x2 = p - x1;
    return min(x1, x2);
}
} // namespace Cipolla

int inv(int a, int b = MOD - 2, int p = MOD) { // 快速幂
    int res = 1;
    for (; b; b >>= 1, a = a * a % p)
        if (b & 1) res = res * a % p;
    return res % p;
}

int init(int n) { // 预处理
    int len = 1;
    while (len < n) len <<= 1;
    return len << 1;
}
int rev[N];
void initrev(int len) {
    for (int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
}

void ntt(int *a, int n, int tp) { // 快速数论变换
    for (int i = 0; i < n; i++)
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int block = 1; block < n; block <<= 1) {
        int siz = block << 1;
        int g1 = inv((tp == 1) ? G : Gi, (MOD - 1) / siz);
        for (int l = 0, r, gn; l < n; l += siz) {
            gn = 1, r = l + block - 1;
            for (int i = l; i <= r; i++, gn = gn * g1 % MOD) {
                int tmp = a[i];
                a[i] = (tmp + gn * a[i + block] % MOD + MOD) % MOD;
                a[i + block] = (tmp - gn * a[i + block] % MOD + MOD) % MOD;
            }
        }
    }
    if (tp == 1) return;
    int orz = inv(n);
    for (int i = 0; i < n; i++) a[i] = (a[i] * orz + MOD) % MOD;
}

void polymul(int *f, int *g, int n, int m, int *tar) { // 多项式乘法
    int len = init(n + m - 1) >> 1;
    initrev(len);
    static int a[N], b[N];
    for (int i = 0; i < n; i++) a[i] = f[i];
    for (int i = 0; i < m; i++) b[i] = g[i];
    for (int i = n; i < len; i++) a[i] = 0;
    for (int i = m; i < len; i++) b[i] = 0;
    ntt(a, len, 1), ntt(b, len, 1);
    for (int i = 0; i < len; i++) tar[i] = a[i] * b[i] % MOD;
    ntt(tar, len, -1);
    for (int i = n + m - 1; i < len; i++) tar[i] = 0;
}

void polyinv(int n, int *f, int *g) { // 多项式乘法逆
    static int tmp[N];
    if (n == 1) {
        g[0] = inv(f[0]);
        return;
    }
    polyinv((n + 1) >> 1, f, g);
    polymul(f, g, n, n, tmp), polymul(tmp, g, n, n, tmp);
    for (int i = 0; i < n; i++) g[i] = (2ll * g[i] - tmp[i] + MOD) % MOD;
}

void polydiv(int *f, int *g, int n, int m, int *q, int *r) { // 多项式除法 and 取余
    static int a[N], b[N], invb[N];
    for (int i = 0; i < n; i++) a[i] = f[i];
    for (int i = 0; i < m; i++) b[i] = g[i];
    reverse(a, a + n), reverse(b, b + m);
    int siz = n - m + 1;
    polyinv(siz, b, invb);
    polymul(a, invb, n, siz, q);
    reverse(q, q + siz);
    int len = init(n);
    for (int i = siz; i < 2ll * len; i++) q[i] = 0;
    polymul(g, q, m, siz, r);
    for (int i = 0; i < m; i++) r[i] = (f[i] - r[i] + MOD) % MOD;
}

void polyderi(int *f, int n, int *g) { // 多项式求导
    for (int i = 1; i < n; i++) g[i - 1] = i * f[i] % MOD;
    g[n - 1] = 0;
}

void polyinte(int *f, int n, int *g) { // 多项式不定积分
    for (int i = 1; i < n; i++) g[i] = f[i - 1] * inv(i) % MOD;
    g[0] = 0;
}

void polyln(int *f, int n, int *g) { // 多项式 ln
    static int a[N], b[N], c[N];
    polyderi(f, n, a), polyinv(n, f, b);
    polymul(a, b, n, n, c);
    polyinte(c, n, g);
    for (int i = n; i < n * 2; i++) g[i] = 0;
}

void polyexp(int *f, int n, int *g) { // 多项式 exp
    static int tmp[N];
    if (n == 1) {
        g[0] = 1;
        return;
    }
    polyexp(f, (n + 1) >> 1, g);
    polyln(g, n, tmp);
    for (int i = 0; i < n; i++) tmp[i] = (f[i] - tmp[i] + MOD) % MOD;
    tmp[0]++;
    polymul(g, tmp, n, n, g);
    for (int i = n; i < n * 2; i++) g[i] = 0;
}

void polypow(int *f, int n, char *k, int *g) { // 多项式快速幂
    if (k[0] == '0') {
        g[0] = 1;
        for (int i = 1; i < n; i++) g[i] = 0;
        return;
    }
    int fst = n, k1 = 0, k2 = 0, k3 = 0, lenk = strlen(k); // 处理巨大恶臭指数
    for (int i = 0; i < lenk; i++) {
        int num = k[i] - '0';
        k1 = (k1 * 10 % MOD + num) % MOD;
        k2 = (k2 * 10 % (MOD - 1) + num) % (MOD - 1);
        if (k3 < n) k3 = k3 * 10 + k[i] - '0';
    }
    if (f[0] == 0 && k3 > n) {
        for (int i = 0; i < n; i++) g[i] = 0;
        return;
    }
    static int tmp[N], a[N];
    for (int i = 0; i < n; i++)
        if (f[i]) {
            fst = i;
            break;
        }
    if (k1 * fst >= n) {
        for (int i = 0; i < n; i++) g[i] = 0;
        return;
    }
    int d = f[fst], invd = inv(d), pd = inv(d, k2);
    for (int i = fst; i < n; i++) a[i - fst] = f[i] * invd % MOD;
    for (int i = n - fst; i < n; i++) a[i] = 0;
    polyln(a, n - fst, tmp);
    for (int i = 0; i < n - fst; i++) tmp[i] = tmp[i] * k1 % MOD;
    polyexp(tmp, n - fst, g + fst * k1);
    for (int i = 0; i < fst * k1; i++) g[i] = 0;
    for (int i = fst * k1; i < n * 2; i++) g[i] = g[i] * pd % MOD;
}

void polysqrt(int *f, int n, int *g) { // 多项式开根
    if (n == 1) {
        g[0] = Cipolla::Cip(f[0], MOD);
        return;
    }
    polysqrt(f, n + 1 >> 1, g);
    static int tmp[N], inv2 = inv(2);
    polyinv(n, g, tmp);
    polymul(f, tmp, n, n, tmp);
    for (int i = 0; i < n; i++) g[i] = (g[i] + tmp[i]) % MOD * inv2 % MOD;
    for (int i = n; i < n * 2; i++) g[i] = 0;
}
posted @ 2021-07-16 14:25  LewisLi  阅读(138)  评论(1)    收藏  举报