基础多项式算法

Posted by Panda2134's Blog on July 1, 2018

约定:设 $\deg(A(x))$ 为多项式 $A(x)$ 的次数。

FFT 和 NTT

可以参见 Miskcoo’s Blog

多项式求逆

定义: 对于 $n$ 次多项式 $A(x)$ ,求出多项式 $B(x)$ ,满足 $A(x) \times B(x) \equiv 1 \pmod{x^p}$.

用数学归纳法的思路。

在 $ \text{mod } x $ 的意义下, $A(x)$ 的逆多项式为 ${a_0}^{-1}$ . 也就是说,要能多项式求逆,首先常数项必须存在逆。在 $\bmod 998244353$ 意义下这等价于 $a_0 \ne 0$.

假定我们已经知道了 $A(x)$ 在 $\text{mod }x^t$ 意义下的逆 $B_0(x)$ ,想要求它在 $\text{mod }x^{2t}$ 意义下的逆 $B(x)$.

由逆多项式定义, $A(x) \times B(x) \equiv 1 \pmod{x^{2t}}$ ,也就是 $1 \dots x^{2t-1}$ 次项系数都是 $0$ ,只有常数项是 $1$ .

那么显然就有

而由定义又有

两式相减,得

显然 $A(x) \ne 0$ ,否则逆多项式不存在。于是一定有

这等价于 $x^t \mid [B(x) - B0(x)]$ .

平方后,有

同乘 $A(x)$ 后:

稍微变形一下,就是:

倍增即可。(想一想,为什么是上取整)

Code:

void poly_inv(int deg, int A[], int B[]) { // mod x^{deg},调用前要保证B[]为空 
    static int tmp[MAXN+10];
    if(deg == 1) {
        B[0] = mod_inv(A[0]);
    } else {
        poly_inv((deg + 1) >> 1, A, B);
        int len;
        for(len = 1; len < 2 * deg; len <<= 1); // 注意:相乘后次数>=n,所以开2倍
        copy(A, A+deg, tmp);
        fill(tmp+deg, tmp+len, 0);
        NTT(len, tmp); NTT(len, B);
        for(int i = 0; i < len; i++)
            B[i] = mul(B[i], sub(2, mul(tmp[i], B[i])));
        NTT(len, B, true);
        fill(B+deg, B+len, 0);
    }
}

多项式除法 / 多项式取模

定义:对于 $n$ 次多项式 $A(x)$ 和 $m$ 次多项式 $B(x)$ ,求出 $n-m$ 次多项式 $Q(x)$ 和小于 $m$ 次的多项式 $R(x)$,使得它们满足下列关系:

思路是考虑用模去掉余数的影响。

取 $A(x)​$ 的系数翻转 $A^R(x) = x^nA(\frac{1}{x})​$ ,于是:

显然 $x^n \times R(1/x)$ 的系数非 $0$ 项次数在 $[n-m+1, n] \cap \mathbb{Z}$ 内。于是放在 $\bmod{x^{n-m+1}}$ 意义下即可消除其影响。

刚好 $Q^R(x)$ 的次数小于 $n-m+1$ ,故可以完整地求出 $Q^R(x)$ ,进而求出 $Q(x)$,剩下的也好办了。

Code:

void poly_sub(int deg, int A[], int B[], int C[]) {
    for(int i = 0; i <= deg; i++) C[i] = sub(A[i], B[i]);
}

void poly_div(int n, int m, int A0[], int B0[], int Q[], int R[]) {
    static int A[MAXN+10], B[MAXN+10], tmp[MAXN+10];
    CLEAR(A); CLEAR(B); CLEAR(tmp);
    copy(A0, A0+n+1, A); copy(B0, B0+m+1, B);
    reverse(A, A+n+1); reverse(B, B+m+1);
    poly_inv(n-m+1, B, tmp);
    conv(n-m, A, tmp, Q);
    fill(Q+n-m+1, Q+n+1, 0);
    reverse(Q, Q+n-m+1);
    copy(A0, A0+n+1, A); copy(B0, B0+m+1, B);
    CLEAR(tmp);
    conv(n, B, Q, tmp);
    poly_sub(n, A, tmp, R);
    fill(R+m, R+n+1, 0);// n+1! past-end iterator! 
}

多项式牛顿迭代

数域为 $\mathbb{R}$ 的牛顿迭代

用途:求解函数 $f(x)$ 的零点。

我们考虑泰勒展开:

取前 $2$ 项来拟合 $f(x)$,求其零点,即为:

不断迭代即可接近零点。

多项式环上的牛顿迭代

用途:已知 $\deg(G(x)) = t$. 求解满足下列方程的 $F(x)$:

考虑类似上文的迭代,每次精度翻倍,即增加了一倍的系数(对于高次项)

当 $n = 1$ 时必须有 $G(F(x)) \equiv 0 \pmod{x}$,这是数学归纳法的基础。

假设我们已知 $F_0(x)$ ,且 $\deg F_0(x) = \left\lceil\frac{n}{2}\right\rceil$. 考虑在 $F_0(x)$ 处对 $G(x)$ 进行泰勒展开:

注意 $\left(F(x) - F_0(x)\right)^2$ 这一项:由于 $F(x) \equiv F_0(x) \pmod{x^{\lceil n/2 \rceil}}$ ,在 $\bmod{x^n}$ 意义下,$F(x)-F_0(x)$ 系数非 $0$ 项次数属于 $\left[\lceil n/2 \rceil, n\right] \cap \mathbb{Z}$.

于是 $\left(F(x) - F_0(x)\right)^2$ 的非 $0$ 项中,最低的一项次数 $\ge n$. 故它即之后的项对答案没有贡献。

故最终的式子是:

多项式 $\ln$

定义:求出多项式 $F(x)$ ,满足 $F(x) \equiv \ln G(x) \pmod{x^n}$.

这里我们不需要采用多项式算法,转而使用微积分解决。

其中 $\frac{G’(x)}{G(x)} = G’(x)\times G^{-1}(x)$,用多项式求逆解决。下面的问题就是如何实现多项式求导和积分。

求导:

积分:

这些都可以 $O(n)$ 实现。

Code:

void poly_int(int deg, int A[], int B[]) {
    B[0] = 0;
    for(int i = 1; i <= deg + 1; i++)
        B[i] = mul(A[i - 1], mod_inv(i));
}

void poly_d(int deg, int A[], int B[]) {
    for(int i = 0; i <= deg - 1; i++)
        B[i] = mul(i + 1, A[i + 1]);
    B[deg] = 0;
}

void poly_ln(int deg, int A0[], int B[]) { // mod x^{deg}
    static int A[MAXN + 10], A1[MAXN + 10], tmp[MAXN + 10];
    fill(A, A + MAXN, 0); fill(A1, A1 + MAXN, 0); fill(tmp, tmp + MAXN, 0);
    copy(A0, A0+deg, A);
    poly_d(deg - 1, A, A1);
    poly_inv(deg, A, tmp);
    copy(tmp, tmp + deg, A);
    conv(deg-1, A, A1, tmp);
    poly_int(deg-1, tmp, B);
}

多项式 $\exp$

定义:求出多项式 $F(x)$,满足 $F(x)\equiv e^{G(x)} \pmod{x^n} \pmod{x^n}$.

变形为 $\ln F(x) - G(x) = 0$ ,然后直接套用多项式牛顿迭代,则有:

使用多项式 $\ln$ 计算 $\ln F_0(x)$ 即可。

Code:

void poly_exp(int deg, int A[], int B[]) { // mod x^{deg}
    static int tmp[MAXN + 10];
    fill(tmp, tmp + MAXN, 0);
    if(deg == 1) {
        B[0] = 1;
    } else {
        int half = (deg + 1) >> 1;
        poly_exp(half, A, B);
        fill(B+half, B+deg, 0);

        poly_ln(deg, B, tmp);
        poly_sub(deg - 1, A, tmp, tmp);
        conv(deg - 1, tmp, B, tmp);
        poly_add(deg - 1, B, tmp, B);
    }
}

多项式开方 / 多项式任意次幂

定义:求出多项式 $F(x)$,满足 $F(x) \equiv {G(x)}^k \pmod{x^n}$.

特殊地,定义多项式开方为 $k=\dfrac{1}{2}$,即 $F(x) \equiv \sqrt{G(x)}$ 的情形。

显然有 $F(x) \equiv \exp(k\ln(G(x)))$,于是应用多项式 $\ln$ 和多项式 $\exp$ 即可解决。

多项式模板大全

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

namespace polynomial {
    const int MAXN = 100000 << 2, MOD = 998244353, G = 3;
    int p_omega[MAXN + 10], p_inv[MAXN + 10];

    inline int pls(int a, int b) {
        int ret = (a + b);
        if(ret >= MOD) ret -= MOD;
        return ret;
    }

    inline int dec(int a, int b) {
        int ret = (a - b);
        if(ret < 0) ret += MOD;
        return ret;
    }

    inline int mul(int a, int b) { return (1ll * a) * b % MOD; }

    inline int fpow(int a, int x) {
        int ret = 1;
        while(x) {
            if(x & 1) ret = mul(ret, a);
            x >>= 1;
            a = mul(a, a);
        }
        return ret;
    }

    inline int mod_inv(int x) { return fpow(x, MOD - 2); }

    //-------------------NTT----------------------

    void init_omega() {
        for(int step = 1; step <= MAXN; step <<= 1) {
            p_omega[step] = fpow(G, (MOD - 1) / step);
            p_inv[step] = fpow(G, (MOD - 1) / step * (step - 1));
        }
    }

    struct _omega_init_t {
        _omega_init_t() { init_omega(); }
    } _omega_init;

    inline int get_rt(int step, bool idft) {
        return !idft ? p_omega[step] : p_inv[step];
    }

    inline void NTT(int len, int A[], bool idft = false) {
        static int rev[MAXN + 10];
        for(int i = 0; i < len; i++)
            rev[i] = (((rev[i >> 1] >> 1) | (len >> (i & 1))) & (len - 1));
        for(int i = 0; i < len; i++)
            if(rev[i] > i) swap(A[i], A[rev[i]]);

        for(int step = 1; step < len; step <<= 1) {
            int rt = get_rt(step << 1, idft);
            for(int i = 0; i < len; i += (step << 1)) {
                for(int j = 0, omega = 1; j < step;
                    j++, omega = mul(omega, rt)) {
                    int t = mul(A[i + j + step], omega);
                    A[i + j + step] = dec(A[i + j], t);
                    A[i + j] = pls(A[i + j], t);
                }
            }
        }

        if(idft) {
            int t = mod_inv(len);
            for(int i = 0; i < len; i++) A[i] = mul(A[i], t);
        }
    }

    inline void conv(int deg, int A0[], int B0[], int C[]) {
        static int A[MAXN + 10], B[MAXN + 10];
        int len;
        copy(A0, A0+deg+1, A), copy(B0, B0+deg+1, B);
        for(len = 1; len < deg * 2 + 1; len <<= 1);
        fill(A+deg+1, A+len, 0), fill(B+deg+1, B+len, 0); // !!
        NTT(len, A), NTT(len, B);
        for(int i = 0; i < len; i++) C[i] = mul(A[i], B[i]);
        NTT(len, C, true);
    }
    //-------------------Polynomial Add & Substract----------

    void poly_add(int deg, int A[], int B[], int C[]) {
        for(int i = 0; i <= deg; i++) C[i] = pls(A[i], B[i]);
    }

    void poly_sub(int deg, int A[], int B[], int C[]) {
        for(int i = 0; i <= deg; i++) C[i] = dec(A[i], B[i]);
    }

    //------------Polynomial Integral & Derivative-----------

    void poly_int(int deg, int A[], int B[]) {
        B[0] = 0;
        for(int i = 1; i <= deg + 1; i++)
            B[i] = mul(A[i - 1], mod_inv(i));
    }

    void poly_d(int deg, int A[], int B[]) {
        for(int i = 0; i <= deg - 1; i++)
            B[i] = mul(i + 1, A[i + 1]);
        B[deg] = 0;
    }

    //-------------------Polynomial inverse------------------

    void poly_inv(int deg, int A[], int B[]) { // mod x^{deg}
        int tmp[MAXN + 10];
        if(deg == 1) {
            assert(A[0] != 0);
            B[0] = mod_inv(A[0]);
        } else {
            poly_inv((deg + 1) >> 1, A, B);
            int len;
            for(len = 1; len < 2 * deg; len <<= 1); // !!
            copy(A, A+deg, tmp); fill(tmp+deg, tmp+len, 0);
            NTT(len, tmp); NTT(len, B);
            for(int i = 0; i < len; i++)
                B[i] = mul(B[i], dec(2, mul(tmp[i], B[i])));
            NTT(len, B, true);
            fill(B+deg, B+len, 0); // !!
        }
    }

    //--------------------Polynomial division-----------------

    void poly_div(int n, int m, int A0[], int B0[], int Q[], int R[]) { // mod x^{deg}
        int A[MAXN + 10], B[MAXN + 10], tmp[MAXN + 10];
        fill(A, A+MAXN, 0); fill(B, B+MAXN, 0); fill(tmp, tmp+MAXN, 0);
        copy(A0, A0+n+1, A), copy(B0, B0+m+1, B);
        reverse(A, A+n+1), reverse(B, B+m+1);
        poly_inv(n-m+1, B, tmp);
        conv(n-m, A, tmp, Q);
        fill(Q+n-m+1, Q+n+1, 0);
        reverse(Q, Q+n-m+1);
        copy(A0, A0+n+1, A), copy(B0, B0+m+1, B);
        conv(n, B, Q, tmp);
        poly_sub(n, A, tmp, R);
        fill(R+m, R+n+1, 0);
    }

    //--------------------Polynomial ln------------------------

    void poly_ln(int deg, int A0[], int B[]) { // mod x^{deg}
        static int A[MAXN + 10], A1[MAXN + 10], tmp[MAXN + 10];
        fill(A, A + MAXN, 0); fill(A1, A1 + MAXN, 0); fill(tmp, tmp + MAXN, 0);
        copy(A0, A0+deg, A);
        poly_d(deg - 1, A, A1);
        poly_inv(deg, A, tmp);
        copy(tmp, tmp + deg, A);
        conv(deg-1, A, A1, tmp);
        poly_int(deg-1, tmp, B);
    }

    //--------------------Polynomial exp-----------------------

    void poly_exp(int deg, int A[], int B[]) { // mod x^{deg}
        static int tmp[MAXN + 10];
        fill(tmp, tmp + MAXN, 0);
        if(deg == 1) {
            B[0] = 1;
        } else {
            int half = (deg + 1) >> 1;
            poly_exp(half, A, B);
            fill(B+half, B+deg, 0);

            poly_ln(deg, B, tmp);
            poly_sub(deg - 1, A, tmp, tmp);
            conv(deg - 1, tmp, B, tmp);
            poly_add(deg - 1, B, tmp, B);
        }
    }
}