板子 多项式全家桶

在所有集训中都是前置知识，快哭了。

记录一些又丑又慢的板子，顺便抄写一些推导过程。

辅助

inline void calc_rev(int n) {
    for (int i = 1; i < n; ++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
}
inline int calc_len(int n) {
    int len = 1;
    while (len < n) len <<= 1;
    return len;
}

for NTT：

inline void clear(int *a, int n) {
    for (int i = 0; i < n; ++i) a[i] = 0;
}
inline void copy(int *a, int *b, int n) {
    for (int i = 0; i < n; ++i) a[i] = b[i];
}
inline void px(int *a, int *b, int n) {
    for (int i = 0; i < n; ++i) a[i] = (LL)a[i] * b[i] % P;
}
inline void init_inv(int n) {
    inv[1] = 1;
    for (int i = 2; i <= n; ++i)
        inv[i] = (LL)inv[P % i] * (P - P / i) % P;
}

for FFT：

const double Pi = acos(-1.0);
struct Complex {
    double x, y;
    Complex(double xx = 0, double yy = 0) {
        x = xx, y = yy;
    }
    Complex operator +(const Complex &t) const {
        return Complex(x + t.x, y + t.y);
    }
    Complex operator -(const Complex &t) const {
        return Complex(x - t.x, y - t.y);
    }
    Complex operator *(const Complex &t) const {
        return Complex(x * t.x - y * t.y, y * t.x + x * t.y);
    }
    Complex conj() {
        return Complex(x, -y);
    }
};
inline void clear(Complex *a, int n) {
    for (int i = 0; i < n; ++i) a[i] = Complex(0, 0);
}
inline void copy(Complex *a, Complex *b, int n) {
    for (int i = 0; i < n; ++i) a[i] = b[i];
}
inline void px(Complex *a, Complex *b, int n) {
    for (int i = 0; i < n; ++i) a[i] = a[i] * b[i];
}

FFT & NTT

FFT

inline void FFT(Complex *f, int n, int opt) {
    int lim = calc_len(n);
    calc_rev(lim);
    for (int i = 0; i < lim; ++i)
        if (i < rev[i]) swap(f[i], f[rev[i]]);
    for (int mid = 1; mid < lim; mid <<= 1) {
        Complex wn(cos(Pi / mid), opt * sin(Pi / mid));
        for (int i = 0; i < lim; i += mid << 1) {
            Complex w(1, 0);
            for (int j = 0; j < mid; ++j, w = w * wn) {
                Complex tx = f[i+j], ty = w * f[i+mid+j];
                f[i+j] = tx + ty, f[i+mid+j] = tx - ty;
            }
        }
    }
    if (opt == -1) {
    	for (int i = 0; i < lim; ++i) f[i].x /= lim, f[i].y /= lim;
    }
}

三次变两次优化：

把 \(f_2\) 放到 \(f_1\) 的虚部上，DFT 一遍，平方，再把虚部取出来除以 2 即可。

NTT

const int MOD = 998244353;
const int G = 3;
const int invG = 332748118;

inline void NTT(int *f, int n, int op) {
    int lim = calc_len(n);
    calc_rev(lim);
    for (int i = 0; i < lim; ++i)
        if (i < rev[i]) swap(f[i], f[rev[i]]);
    for (int mid = 1; mid < lim; mid <<= 1) {
        int wn = power(op == 1 ? G : invG, (P - 1) / (mid << 1));
        for (int i = 0; i < lim; i += mid << 1) {
            int w = 1;
            for (int j = 0; j < mid; ++j) {
                int tx = f[i + j], ty = (LL)w * f[i + mid + j] % P;
                f[i + j] = (tx + ty) % P;
                f[i + mid + j] = (tx - ty + P) % P;
                w = (LL)w * wn % P;
            }
        }
    }
    if (op == -1) {
        int invN = power(lim, P - 2);
        for (int i = 0; i < lim; ++i) f[i] = (LL)f[i] * invN % P;
    }
}

多项式乘法

inline void poly_times(int *f, int *g, int n, int lim) {
    static int t[N];
    int len = calc_len(n << 1);
    clear(t, len), copy(t, g, n);
    NTT(f, len, 1), NTT(t, len, 1);
    px(f, t, len), NTT(f, len, -1);
    clear(f + lim, len - lim);
}

多项式求逆

以下是不知道什么时候写的：引用了 autoint 的这篇博客。 \[ \begin{aligned} A\times B&\equiv 1 \pmod{x^n}\\ A\times B'&\equiv 1 \pmod{x^{\frac{n}{2}}}\\ A\times (B-B')&\equiv 0\pmod{x^{\frac{n}{2}}}\\ B-B'&\equiv 0\pmod{x^{\frac{n}{2}}}\\ (B-B')^2&\equiv 0\pmod{x^n}\\ A\times(B^2-2BB'+B'^2)&\equiv 0\pmod{x^n}\\ B-2B'+AB'^2&\equiv 0\pmod{x^n}\\ B&\equiv 2B'-AB'^2\pmod{x^n} \end{aligned} \] 据此可以倍增（可以把 B 数组滚动）。

新笔记： \[ A(x)B(x)-1\equiv 0\pmod{x^n} \] 设 \(G(B(x))=A(x)B(x)-1\)，求导得 \(G'(B(x))=A(x)\)。

已知 \(B_*(x)\equiv A(x)^{-1}\pmod{x^{\frac{n}{2}}}\)，牛顿迭代得 \[ \begin{aligned} B(x)&\equiv B_*(x)-\frac{G(B_*(x))}{G'(B_*(x))}\pmod{x^n}\\ &\equiv B_*(x)-\frac{A(x)B_*(x)-1}{A(x)}\\ &\equiv B_*(x)-B(x)(A(x)B_*(x)-1) \end{aligned} \] \(A(x)B_*(x)-1\) 在 \(\bmod x^n\) 意义下只有次数大于等于 \(\frac{n}{2}\) 的项（因为在 \(\bmod x^{\frac{n}{2}}\) 意义下，\(A(x)B_*(x)-1\) 次数小于 \(\frac{n}{2}\) 的项都为 \(0\)），所以 \(B(x)\) 次数大于等于 \(\frac{n}{2}\) 的项与其相乘后就会被模掉，只有次数小于 \(\frac{n}{2}\) 的项有意义，便可以用 \(B_*(x)\) 代替 \(B(x)\)。 \[ B(x)\equiv 2B_*(x)-A(x)B_*(x)^2\pmod{x^n} \] 然后倍增，复杂度 \(\mathcal O(n\log n)\)，证明就是 \(\mathcal T(n)=\mathcal T(\frac{n}{2})+\mathcal O(n\log n)\)。

inline void poly_inv(int *f, int n) {
    static int b1[N], b2[N], b3[N];
    int lim = calc_len(n);
    clear(b1, lim << 1), clear(b2, lim << 1), clear(b3, lim << 1);
    b1[0] = power(f[0], P - 2);
    for (int len = 2; len <= lim; len <<= 1) {
        for (int i = 0; i < (len >> 1); ++i)
            b2[i] = 2ll * b1[i] % P;
        copy(b3, f, len);
        NTT(b1, len << 1, 1), NTT(b3, len << 1, 1);
        px(b1, b1, len << 1), px(b1, b3, len << 1);
        NTT(b1, len << 1, -1), clear(b1 + len, len);
        for (int i = 0; i < len; ++i)
            b1[i] = (b2[i] - b1[i] + P) % P;
    }
    copy(f, b1, n);
}

多项式开根

要求 \(f[0] = 1\)。 \[ B^2(x)-A(x)= 0 \]

设 \(G(B(x))=B^2(x)-A(x)\)，则 \[ G'(B(x))=2B(x) \] 牛顿迭代 \[ \begin{aligned} B(x)&=B_*(x)-\frac{G(B_*(x))}{G'(B_*(x))}\\ &=B_*(x)-\frac{B_*^2(x)-A(x)}{2B_*(x)}\\ &=\frac{B_*^2(x)+A(x)}{2B_*(x)} \end{aligned} \] 同样地倍增。

inline void poly_sqrt(int *f, int n) {
    static int b1[N], b2[N];
    int lim = calc_len(n);
    clear(b1, lim << 1), clear(b2, lim << 1);
    b1[0] = 1;
    for (int len = 2; len <= lim; len <<= 1) {
        for (int i = 0; i < (len >> 1); ++i) b2[i] = 2ll * b1[i] % P;
        poly_inv(b2, len);
        NTT(b1, len, 1), px(b1, b1, len), NTT(b1, len, -1);
        for (int i = 0; i < len; ++i) b1[i] = (b1[i] + f[i]) % P;
        poly_times(b1, b2, len, len);
    }
    copy(f, b1, n);
}

多项式带余除法

\[ F(x)=Q(x)\times G(x)+R(x) \]

将多项式翻转，即设 \(F^T(x)=x^nF(x^{-1})\)，则 \[ F^T(x)=Q^T(x)\times G^T(x)+x^{n-m+1}R^T(x) \] 想办法消去余数 \(R(x)\)，就对 \(x^{n-m+1}\) 取模 \[ \begin{aligned} F^T(x)&\equiv Q^T(x)\times G^T(x)&\pmod{x^{n-m+1}}\\ Q^T(x)&\equiv F^T(x)\times G^T(x)^{-1}&\pmod{x^{n-m+1}} \end{aligned} \]

求出 \(Q^T(x)\) 后再翻转即可。

那么余数 \(R(x)=F(x)-Q(x)\times G(x)\)。

inline void poly_mod(int *f, int *g, int n, int m) {
    static int a[N], b[N];
    clear(a, n), clear(b, n);
    int L = n - m + 1;
    reverse(f, f + n), copy(a, f, L), reverse(f, f + n);
    reverse(g, g + m), copy(b, g, L), reverse(g, g + m);
    poly_inv(b, L), poly_times(b, a, L, L), reverse(b, b + L);
    poly_times(g, b, n, n);
    for (int i = 0; i < m - 1; ++i) g[i] = (f[i] - g[i] + P) % P;
    clear(g + m - 1, L);
    copy(f, b, L), clear(f + L, n - L);
}
// f -> Q(x), g -> R(x)

多项式求导 & 积分

求导

单项式求导：\((ax^k)'=akx^{k-1}\)，推广一下。

inline void poly_dev(int *f, int n) {
    for (int i = 1; i < n; ++i)
        f[i - 1] = (LL)f[i] * i % P;
    f[n - 1] = 0;
}

积分

逆运算，记得预处理逆元。

inline void poly_int(int *f, int n) {
    for (int i = n; i; --i)
        f[i] = (LL)f[i - 1] * inv[i] % P;
    f[0] = 0;
}

多项式 ln

\[ \ln(A(x))=B(x) \]

已知 \(\ln' x=\frac{1}{x}\)，所以两边分别对 \(x\) 求导，并运用链式法则 \[ \begin{aligned} \frac{d}{dA(x)}\ln(A(x))\frac{dA(x)}{x}&=B'(x)\\ A(x)^{-1}A'(x)&=B'(x)\\ B(x)&=\int(A'(x)A(x)^{-1}) \end{aligned} \] 要求常数项为 \(1\)。

inline void poly_ln(int *f, int n) {
    static int t[N];
    copy(t, f, n);
    poly_dev(t, n), poly_inv(f, n);
    poly_times(f, t, n, n);
    poly_int(f, n - 1);
}

多项式 exp

牛顿迭代

\[ e^{A(x)}=B(x) \]

设 \(G(B(x))=\ln(B(x))-A(x)\)，就可以牛顿迭代了。

求导得 \(G'(B(x))=B(x)^{-1}\)，则 \[ \begin{aligned} B(x)&=B_*(x)-\frac{G(B_*(x))}{G'(B_*(x))}\\ &=B_*(x)-\frac{\ln(B_*(x))-A(x)}{B_*(x)^{-1}}\\ &=B_*(x)(1-\ln(B_*(x))+A(x)) \end{aligned} \] 倍增即可，不过常数很大，大概 \(\mathcal O(48\cdot n\log n)\)。

要求常数项为 \(0\)。

inline void poly_exp(int *f, int n) {
    static int b1[N], b2[N];
    int lim = calc_len(n);
    clear(b1, lim), clear(b2, lim);
    b1[0] = 1;
    for (int len = 2; len <= lim; len <<= 1) {
        copy(b2, b1, len >> 1);
        poly_ln(b2, len);
        for (int i = 0; i < len; ++i) b2[i] = (f[i] - b2[i] + P) % P;
        b2[0] = (b2[0] + 1) % P;
        poly_times(b1, b2, len, len);
    }
    copy(f, b1, n);
}

分治 FFT

设 \(B(x)=e^{A(x)}\)，两边求导得 \[ B'(x)=\frac{d}{dA(x)}e^{A(x)}\frac{dA(x)}{x}=B(x)A'(x) \] 提取第 \(n\) 项系数 \[ (n+1)B[n+1]=\sum_{i=0}^n(i+1)A[i+1]B[n-i] \] 整理一下就是 \[ B[n]=\frac{1}{n}\sum_{i=1}^n iA[i]B[n-i] \] 套上分治 FFT，复杂度 \(\mathcal O(n\log^2 n)\)，但是比牛顿迭代跑的快！

int f[N], g[N], a[N], b[N];
void CDQ(int l, int r) {
    if (r - l <= 1) {
        if (l > 0) f[l] = (LL)f[l] * inv[l] % P;
        return;
    }
    int len = r - l, mid = (l + r) >> 1;
    CDQ(l, mid);
    copy(a, f + l, len >> 1), clear(a + (len >> 1), len >> 1);
    copy(b, g, len);
    NTT(a, len, 1), NTT(b, len, 1);
    px(a, b, len), NTT(a, len, -1);
    for (int i = len >> 1; i < len; ++i)
        f[l + i] = (f[l + i] + a[i]) % P;
    CDQ(mid, r);
}
// main
calc_inv(n);
for (int i = 0; i < n; ++i) g[i] = (LL)g[i] * i % P;
int len = calc_len(n);
f[0] = 1;
CDQ(0, len);

多项式快速幂

在保证 \(f[0]=1\) 的情况下。

先 \(\ln\) 再 \(\exp\)，很好理解。

inline void poly_pow(int *f, int n, int k) {
    poly_ln(f, n);
    for (int i = 0; i < n; ++i) f[i] = (LL)f[i] * k % P;
    poly_exp(f, n);
}

在不保证 \(f[0]=1\) 的情况下。

\(f[0]\not=1\) 时不能直接 \(\ln,\exp\)，想办法让 \(f[0]=1\)。

那就直接让所有项除以 \(f[0]\)。

如果 \(f[0]=0\)？那就找到最低次的非零项 \(x^p\)，所有项提出个 \(x^p\) 的公因子。

设 \(c=f[0]\)，即 \[ A^k(x)=(B(x)\times cx^p)^k=B^k(x)\times c^kx^{pk} \] 根据费马小定理，\(c^k\) 的 \(k\) 是在 \(\bmod \varphi(p)\) 意义下的。

inline void poly_pow2(int *f, int n, char *str) {
    int p = 0, k = 0, k2 = 0;
    while (!f[p]) ++p;
    for (int i = 0; str[i]; ++i) {
        k = (10ll * k + str[i] - '0') % P;
        k2 = (10ll * k2 + str[i] - '0') % (P - 1);
        if ((LL)p * k > (LL)n) {
            for (int i = 0; i < n; ++i) f[i] = 0;
            return;
        }
    }
    int m = n - p * k;
    int c = f[p], invc = power(c, P - 2);
    for (int i = 0; i < m; ++i)
        f[i] = (LL)f[i + p] * invc % P;
    clear(f + m, p * k);
    poly_ln(f, m);
    for (int i = 0; i < m; ++i)
        f[i] = (LL)f[i] * k % P;
    poly_exp(f, m);
    c = power(c, k2);
    for (int i = n - 1; i >= p * k; --i)
        f[i] = (LL)f[i - p * k] * c % P;
    for (int i = 0; i < p * k; ++i) f[i] = 0;
}

拆系数 FFT（MTT）

4 次 FFT 的做法。

先拆系数，\(f(x)=base\times A(x)+B(x)\)，\(g(x)=base\times C(x)+D(x)\)，一般 \(base=\sqrt{p}\)，选择 \(2^{15}\) 便于提取系数。

两个多项式拆成四个后直接 DFT 很浪费，发现虚部都是空的，两两打包一块 DFT，即 \[ \begin{aligned} P(x)&=A(x)+iB(x)\\ Q(x)&=A(x)-iB(x) \end{aligned} \] AzusaCat 的博客证明了 \(P[k]\) 和 \(Q[n-k]\) 是共轭的，所以只需要 DFT \(P(x)\)。

DFT 后解二元一次方程组即可得到 \(A(x),B(x)\) 的点值。

iDFT 一样，将 \(A(x)C(x)\)，\(A(x)D(x)+B(x)C(x)\)，\(B(x)D(x)\) 相乘后打包 iDFT 即可。

inline void MTT(int *f, int *g, int n, int m) {
    static Complex a[N], b[N], c[N], d[N];
    int lim = calc_len(n + m + 1);
    for (int i = 0; i < lim; ++i) {
        a[i].x = f[i] >> 15, a[i].y = f[i] & 32767;
        c[i].x = g[i] >> 15, c[i].y = g[i] & 32767;
    }
    FFT(a, lim, 1), FFT(c, lim, 1);
    for (int i = 1; i < lim; ++i) b[i] = a[lim - i].conj();
    b[0] = a[0].conj(); /* 特判 */
    for (int i = 1; i < lim; ++i) d[i] = c[lim - i].conj();
    d[0] = c[0].conj();
    for (int i = 0; i < lim; ++i) {
        Complex aa = (a[i] + b[i]) * Complex(0.5, 0),
            bb = (a[i] - b[i]) * Complex(0, -0.5),
            cc = (c[i] + d[i]) * Complex(0.5, 0),
            dd = (c[i] - d[i]) * Complex(0, -0.5); /* 乘(0,-0.5)就是除以2i */
        a[i] = aa * cc + Complex(0, 1) * (aa * dd + bb * cc);
        b[i] = bb * dd;
    }
    FFT(a, lim, -1), FFT(b, lim, -1);
    for (int i = 0; i <= n + m; ++i) {
        int t1 = (LL)(a[i].x + 0.5) % P, 
            t2 = (LL)(a[i].y + 0.5) % P, 
            t3 = (LL)(b[i].x + 0.5) % P; /* 注意精度T_T */
        f[i] = (((1ll << 30) * t1 + (1ll << 15) * t2 + t3) % P + P) % P;
    }
}

分治 FFT

给出 \(g[1]\dots g[n-1]\)，求 \(f[0]\dots f[n-1]\)。 \[ f[i]=\sum_{j=1}^if[i-j]g[i] \] \(f\) 每项会被比它小的项贡献，考虑 CDQ 分治，中序遍历处理左边对右边的贡献。

预先把长度补成 \(2\) 的次幂，会好写些。（你这明明是分治 NTT 啊。）

int f[N], g[N], a[N], b[N];
void CDQ(int l, int r) {
    if (r - l <= 1) return;
    int len = r - l, mid = (l + r) >> 1;
    CDQ(l, mid);
    copy(a, f + l, len >> 1), clear(a + (len >> 1), len >> 1);
    copy(b, g, len);
    NTT(a, len, 1), NTT(b, len, 1);
    px(a, b, len), NTT(a, len, -1);
    for (int i = len >> 1; i < len; ++i)
        f[l + i] = (f[l + i] + a[i]) % P;
    CDQ(mid, r);
}
// main
int len = calc_len(n);
f[0] = 1;
CDQ(1, len);

RenaMoe's Blog

板子多项式全家桶

辅助

FFT & NTT

FFT

NTT

多项式乘法

多项式求逆

多项式开根

多项式带余除法

多项式求导 & 积分

求导

积分

多项式 ln

多项式 exp

牛顿迭代

分治 FFT

多项式快速幂

拆系数 FFT（MTT）

分治 FFT