普通有限多项式笔记

合集 - 算法竞赛相关(18)

1.数论分块笔记2023-12-072.线性空间笔记2023-12-073.SG定理证明2023-12-074.2-SAT 问题入门2023-12-075.The Main Idea of Basic Dynamic Programming Side A2023-12-076.四边形不等式笔记2023-12-087.杨氏矩阵笔记2023-12-088.网络流复习笔记2023-12-089.Kirchhoff 矩阵树定理的无向图情况2023-12-0810.网络流部分结论性质及证明2023-12-2011.SA&SAM 小记2024-01-1712.置换小记2024-06-2713.欧拉五边形数定理小记2024-04-0314.QOJ #1280.Fibonacci Partition/Fibonacci性质大杂烩2024-04-0515.生成函数小记2024-04-11

16.普通有限多项式笔记2024-04-25

17.二项式反演小记2024-06-0718.下降幂及斯特林数杂谈2024-06-10

收起

普通多项式笔记

\(\textrm{Newton's Method}\) ，牛顿迭代

应用于解决已知 \(g(x)\) 的情况下，求出 \(g(f(x))\equiv 0\mod x^n\)。

首先通过列出方程显然，\(f(x) \mod x^n\) 在此时是唯一的。

那么我们假设已知 \(g(f_0(x))\equiv 0\mod x^{n/2}\)，显然此时 \(f_0(x)\mod x^{n/2}\) 也是唯一的，且 \(\deg f_0(x)=\frac{n}{2},f(x)\equiv f_0(x) \mod x^{n/2}\)。

这时我们把 \(g(f(x))\) 在 \(f_0(x)\) 处展开，可以得到：

\[g(f(x))= \sum_{k=0}^{+\infty}\frac{g^{(k)}(f_0(x))}{k!}(f(x)-f_0(x))^k \]

因此有：

\[\sum_{k=0}^{+\infty}\frac{g^{(k)}(f_0(x))}{k!}(f(x)-f_0(x))^k\equiv 0 \mod x^n \]

因为 \(f(x)-f_0(x)\equiv 0 \mod x^{n/2}\)，所以有： \((f(x)-f_0(x))^2 \equiv 0 \mod x^n\)。

当然 \(n=1\) 时应该要特殊地求出 \(f(x)\)。

所以：

\[\sum_{k=0}^{+\infty}\frac{g^{(k)}(f_0(x))}{k!}(f(x)-f_0(x))^k\equiv g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\equiv 0 \mod x^n \\ f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))} \mod x^n \]

最终的结果看起来没什么令人惊奇的地方，但是实际上应用却是非常多的。

\(\textrm{inverse}\)，逆

已知 \(h(x)\)，求 \(f(x)=h^{-1}(x)\)。

考虑使用 \(\textrm{Newton's Method}\) 解决这个问题。

构造 \(g(f(x))=f^{-1}(x)-h(x)\)，此时我们把 \(h(x)\) 当作一个常数。

因此有：

\[\begin{aligned} f(x)&\equiv f_0(x)-\frac{f^{-1}_0(x)-h(x)}{-f_0^{-2}(x)} &\mod x^n \\ &\equiv f_0(x)+f_0^2(x)\left(f_0^{-1}(x)-h(x)\right) &\mod x^n \end{aligned} \]

时间复杂度 \(O(n\log_2 n)\)。

\(\textrm{division}\)，带余除法

已知 \(A(x),B(x)\)，求 \(Q(x),R(x)\) 使得 \(A(x)=B(x)\cdot Q(x)+R(x)\)。

记 \(\deg A(x)=n,\deg B(x)=m,n\ge m\)，那么钦定 \(\deg Q(x)=n-m,\deg R(x)=m-1\)。

记对多项式 \(f(x)\) 的 \(\textrm{Reverse}\) 操作为 \(f^R(x)=x^{\deg f(x)}f\left(\frac1 x\right)\)。

因此有：

\[A(\frac{1}{x})=B(\frac{1}{x})\cdot Q(\frac{1}{x}) +R(\frac{1}{x}) \\ x^nA(\frac{1}{x})=x^{m}B(\frac{1}{x})\cdot x^{n-m}Q(\frac{1}{x}) +x^{n-m+1}\cdot x^{m-1}R(\frac{1}{x}) \\ A^R(x)=B^R(x)\cdot Q^R(x)+x^{n-m+1} R^R(x) \\ A^R(x)\equiv B^R(x)\cdot Q^R(x) \mod x^{n-m+1} \]

可以求出 \(Q^R(x)\) 进而求出 \(Q(x),R(x)\)。

\(\textrm{sqrt}\)，开方

已知 \(h(x)\)，求 \(f(x)=\sqrt {h(x)}\)。

考虑使用 \(\textrm{Newton's Method}\) 解决这个问题。

构造 \(g(f(x))=f^2(x)-h(x)\)，此时我们把 \(h(x)\) 当作一个常数。

因此有：

\[\begin{aligned} f(x)\equiv f_0(x)-\frac{f_0^2(x)-h(x)}{2f_0(x)} \mod x^n \end{aligned} \]

时间复杂度 \(O(n\log_2n)\)。

\(\ln\)，自然对数

已知 \(f(x)\)，求 \(\ln f(x)\)。

考虑求导：

\[(\ln f(x))'=\frac{f'(x)}{f(x)} \]

积分还原后：

\[\ln f(x)=\int \frac{f'(x)}{f(x)} \]

有一个问题是最终答案的常数项，但通常认为 \(\ln\) 存在的充要条件是常数项为 \(1\)。

\(\exp\)，自然常数的幂

已知 \(h(x)\)，求 \(f(x)=\exp h(x)\)。

考虑使用 \(\textrm{Newton's Method}\) 解决这个问题。

构造 \(g(f(x))=\ln f(x)-h(x)\)，此时我们把 \(h(x)\) 当作一个常数。

因此有：

\[\begin{aligned} f(x)&\equiv f_0(x)-\frac{\ln f_0(x)-h(x)}{f_0^{-1}(x)} &\mod x^n \\ &\equiv f_0(x)(1-\ln f_0(x)+h(x)) &\mod x^n \end{aligned} \]

时间复杂度 \(O(n\log_2 n)\)。

需要注意的是 \(h(x)\) 的常数项一定为 \(0\)，否则如果使用快速数论变换时会出现问题，因为我们无法得到 \(\exp [x^0]h(x)\) 在取模意义下的取值。

\(\textrm{power}\)，幂

已知 \(h(x)\)，求 \(h^m(x)\)。

显然有：

\[h^m(x)=\exp (m \ln h(x)) \]

当然，如果出现 \([x^0]h(x)=0\) 或者 \([x^0]h(x)\ne 1,0\) 的情况要先对 \(h(x)\) 做处理，最后再乘回来。

时间复杂度 \(O(n\log_2 n)\)。

全家桶代码

// Not afraid to dark.

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

clock_t start_time, end_time;
#define GET_START start_time = clock ();
#define GET_END end_time = clock (); fprintf (stderr, "TIME COSSEMED : %0.3lf\n", 1.0 * (end_time - start_time) / CLOCKS_PER_SEC);

#define int long long

namespace io {
    int read_pos, read_dt; char read_char;
    __inline__ __attribute__ ((always_inline)) int read (int &p = read_pos){
        p = 0, read_dt = 1; read_char = getchar ();
        while (! isdigit (read_char)){
            if (read_char == '-')
                read_dt = - 1;
            read_char = getchar ();
        }
        while (isdigit (read_char)){
            p = (p << 1) + (p << 3) + read_char - 48;
            read_char = getchar ();
        }
        return p = p * read_dt;
    }
    int write_sta[65], write_top;
    __inline__ __attribute__ ((always_inline)) void write (int x){
        if (x < 0)
            putchar ('-'), x = - x;
        write_top = 0;
        do
            write_sta[write_top ++] = x % 10, x /= 10;
        while (x);
        while (write_top)
            putchar (write_sta[-- write_top] + 48);
    }
    int llen;
    __inline__ __attribute__ ((always_inline)) int get_string (char c[], int &len = llen){
        len = 0;
        read_char = getchar ();
        while (read_char == ' ' || read_char == '\n' || read_char == '\r')
            read_char = getchar ();
        while (read_char != ' ' && read_char != '\n' && read_char != '\r'){
            c[++ len] = read_char;
            read_char = getchar ();
        }
        return len;
    }
}

#define log2_(x) ((x) ? (31 ^ __builtin_clz (x)) : 0)

namespace Region {
    const int mod = 998244353, G = 3, Gi = 332748118;
    __inline__ __attribute__ ((always_inline)) int power (int a, int p){
        int res = 1;
        while (p > 0){
            if (p & 1)
                res = res * a % mod;
            a = a * a % mod;
            p >>= 1;
        }
        return res;
    }
    __inline__ __attribute__ ((always_inline)) int Inverse (int x){
        return power (x, mod - 2);
    }
}

namespace Polynomial {
    // ordinary Number-Theoretic Tranform

    using namespace Region;

    const int N = 5e5;

    __inline__ __attribute__ ((always_inline)) void Reverse (vector < int > &a, int len){
        static int rev[N + 5];
        rev[0] = 0;
        for (int i = 1;i < len;++ i){
            rev[i] = rev[i >> 1] >> 1;
            if (i & 1)
                rev[i] |= len >> 1;
            if (i < rev[i])
                swap (a[i], a[rev[i]]);
        }
    }

    struct pnm {
        int len;
        vector < int > a;
        pnm (int _len_ = 0, int _c_ = 0){
            len = _len_;
            a.resize (len, 0);
            if (len)
                a[0] = _c_;
        }
        __inline__ __attribute__ ((always_inline)) int & operator [] (const int x) {
            return a[x];
        }
        __inline__ __attribute__ ((always_inline)) int operator [] (const int x) const {
            return a[x];
        }
        __inline__ __attribute__ ((always_inline)) int deg (){
            for (int i = len - 1;i > 0;-- i)
                if (a[i] != 0)
                    return i;
            return 0;
        }
        __inline__ __attribute__ ((always_inline)) void flush (int x = 0){
            if (x < len)
                x = len;
            for (int i = 2 << log2_ (x - 1);i > len;-- i)
                a.push_back (0);
            len = (int) a.size ();
        }
        __inline__ __attribute__ ((always_inline)) void release (){
            len = 0, a.clear ();
        }
        __inline__ __attribute__ ((always_inline)) pnm Cut (int l) const {
            pnm res (l);
            for (int i = min (l, len) - 1;i >= 0;-- i)
                res[i] = a[i];
            return res;
        }
        __inline__ __attribute__ ((always_inline)) void Delete (int l){
            while (len > l){
                -- len;
                a.pop_back ();
            }
        }
        __inline__ __attribute__ ((always_inline)) void NTT (int swi){
            Reverse (a, len);
            int om, mid, x, y;
            for (int l = 2;l <= len;l <<= 1){
                om = power ((swi == - 1) ? Gi : G, (mod - 1) / l), mid = l >> 1;
                for (int i = 0;i < len;i += l)
                    for (int j = i, pw = 1;j < i + mid;++ j, pw = pw * om % mod){
                        x = a[j], y = a[j + mid] * pw % mod;
                        a[j] = (x + y) % mod;
                        a[j + mid] = (x - y + mod) % mod;
                    }
            }
            if (swi == - 1){
                int iv = Inverse (len);
                for (int i = 0;i < len;++ i)
                    a[i] = a[i] * iv % mod;
            }
        }
    };

    __inline__ __attribute__ ((always_inline)) pnm operator + (const pnm a, const pnm b){
        pnm res (max (a.len, b.len));
        for (int i = 0;i < a.len;++ i)
            res[i] = a[i];
        for (int i = 0;i < b.len;++ i)
            res[i] = (res[i] + b[i]) % mod;
        return res;
    }
    __inline__ __attribute__ ((always_inline)) pnm operator - (const pnm a, const pnm b){
        pnm res (max (a.len, b.len));
        for (int i = 0;i < a.len;++ i)
            res[i] = a[i];
        for (int i = 0;i < b.len;++ i)
            res[i] = (res[i] - b[i] + mod) % mod;
        return res;
    }
    __inline__ __attribute__ ((always_inline)) pnm operator ^ (const pnm a, const pnm b){
        // assert (a.len == b.len);
        pnm res = a;
        for (int i = 0;i < b.len;++ i)
            res[i] = res[i] * b[i];
        return res;
    }
    __inline__ __attribute__ ((always_inline)) pnm operator * (pnm a, pnm b){
        static int l = 0;
        l = a.len + b.len - 1;
        a.flush (l), b.flush (l);
        l = a.len;
        a.NTT (1), b.NTT (1);
        a = a ^ b;
        a.NTT (- 1);
        return a;
    }
    __inline__ __attribute__ ((always_inline)) pnm inv (pnm f, int le = 0){
        // assert (f[0] != 0);
        f.flush ();
        if (! le)
            le = f.len;
        pnm iv (1, Inverse (f[0]));
        for (int l = 2;l <= le;l <<= 1)
            iv = iv * (pnm (1, 2) - (f.Cut (l) * iv).Cut (l)), iv.Delete (l);
        return iv;
    }
    __inline__ __attribute__ ((always_inline)) pnm operator / (const pnm a, pnm b){
        // assert (b[0] != 0);
        return a * inv (b);
    }
    __inline__ __attribute__ ((always_inline)) void Rev (pnm &f, int l = 0){
        if (! l)
            l = f.len;
        f.flush (l);
        for (int i = 0;i < l - i - 1;++ i)
            swap (f[i], f[l - i - 1]);
    }
    __inline__ __attribute__ ((always_inline)) pnm operator | (pnm a, pnm b){
        int nn = a.deg (), mm = b.deg ();
        Rev (a, nn + 1), Rev (b, mm + 1);
        pnm c = a.Cut (nn - mm + 1) * inv (b.Cut (nn - mm + 1));
        c.Delete (nn - mm + 1);
        Rev (c, nn - mm + 1);
        return c;
    }
    __inline__ __attribute__ ((always_inline)) pnm operator % (const pnm a, const pnm b){
        return a - ((a | b) * b);
    }

    __inline__ __attribute__ ((always_inline)) pnm inte (pnm f){
        if (f.len == 0)
            return f;
        for (int i = 1;i < f.len;++ i)
            f[i - 1] = f[i] * i % mod;
        -- f.len, f.a.pop_back ();
        return f;
    }
    __inline__ __attribute__ ((always_inline)) pnm deri (pnm f, int _c = 0){
        static int jc[N + 5], inj;
        jc[0] = 1;
        for (int i = 1;i <= f.len;++ i)
            jc[i] = jc[i - 1] * i % mod;
        inj = Inverse (jc[f.len]);
        ++ f.len;
        f.a.push_back (0);
        for (int i = f.len - 1;i > 0;inj = inj * (i --) % mod)
            f[i] = f[i - 1] * jc[i - 1] % mod * inj % mod;
        f[0] = _c;
        return f;
    }
    __inline__ __attribute__ ((always_inline)) pnm LN (pnm f, int m = 0){
        // assert (f[0] == 1);
        f.flush ();
        if (! m)
            m = f.len;
        return deri (inte (f) * inv (f, m)).Cut (m);
    }
    __inline__ __attribute__ ((always_inline)) pnm EXP (pnm f, int m = 0){
        // assert (f[0] == 0);
        f.flush ();
        static int p;
        pnm g (1, 1);
        p = 1;
        if (! m)
            m = f.len;
        while (p < m){
            p <<= 1;
            g = g * (pnm (1, 1) - LN (g, p) + f.Cut (p));
            g.Delete (p);
        }
        g.Delete (m);
        return g;
    }
    __inline__ __attribute__ ((always_inline)) pnm SQRT (pnm f, int m = 0){
        // assert (f[0] == 1)
        f.flush ();
        static int p;
        pnm g (1, 1);
        p = 1;
        if (! m)
            m = f.len;
        while (p < m){
            p <<= 1;
            g = g - ((g * g).Cut (p) - f.Cut (p)) * pnm (1, Inverse (2)) * inv (g, p);
            g.Delete (p);
        }
        g.Delete (m);
        return g;
    }
    __inline__ __attribute__ ((always_inline)) pnm POWER (pnm f, int m, int n = 0){
        f.flush ();
        if (! n)
            n = f.len;
        int lowx = 0;
        while (lowx < f.len && ! f[lowx])
            ++ lowx;
        if (lowx == f.len)
            return pnm (n);
        if (lowx && (m >= mod || m * lowx > n))
            return pnm (n);
        int v = f[lowx];
        int iv = Inverse (v);
        for (int i = 0;i < f.len;++ i){
            if (i + lowx < f.len)
                f[i] = f[i + lowx] * iv % mod;
            else
                f[i] = 0;
        }
        f = EXP (pnm (1, m % mod) * LN (f), n).Cut (n);
        if (lowx){
            pnm g (m % mod * lowx + 1);
            g[m % mod * lowx] = power (v, m % (mod - 1));
            return f * g;
        }
        return f * pnm (1, power (v, m % (mod - 1)));
    }
}

using Polynomial::EXP;
using Polynomial::LN;
using Polynomial::POWER;
using Polynomial::SQRT;

signed main (){
    // GET_START

    

    // GET_END
    return 0;
}

循环卷积

就是求：

\[[x^m]h(x)=\sum_{p+q\equiv r\mod n}[x^p]f(x)[x^q]g(x) \]

实际上十分简单。

令 \(h'(x)=f(x)g(x)\)，会发现显然 \(h'(x)\) 和 \(h(x)\) 在 \(\forall 0\le v<n:\omega_n^v\) 处的取值是一样的。原理是：

\[[x^m]h(x)=[x^m]h'(x)+[x^{m+n}]h'(x)+[x^{m+2n}]h'(x)+\cdots \]

而实际上对于 \(x=\omega_n^v\)，等式后面的这些项的贡献 \(\omega_n^{mv},\omega^{mv+nv}\omega^{mv+2nv},\cdots\) 是完全相等的，等于 \(\omega_n^{mv}\)。

也就是说在 \(x=\omega_n^v\) 时分别代表 \(h(x)\) 和 \(h'(x)\) 的等式两边的贡献是一样的。

因此 \(h(\omega_n^v)=h'(w_n^v)\)。

因此我们只需对 \(f(x),g(x)\) 做长度为 \(n\) 的 \(\textrm{DFT}\) 而后各项分别对应乘起来再 \(\textrm{IDFT}\) 得到 \(h(x)\)。

问题在于当 \(n\) 不再是 \(2\) 的幂次时如何解决。

见下文。

长度任意时的快速变换与逆变换

类似 \(\textrm{FFT}\) 地，对于 \(n\) 的一个因数 \(c\)：

\[\begin{aligned} &F(\omega_n^k)=\sum_{l=0}^{c-1} \omega_{n}^{kl}F_l(\omega_{n/c}^k) \end{aligned} \]

因此我们有：

\[T(n)=O(nc)+T(\frac{n}{c}) \]

最优的方法是将 \(n\) 质因数分解后 \(c\) 从小到大依次取质因数，最坏情况可以到达 \(O(n^2)\)。

逆变换亦然。

分治 FFT

给出 \(g_{1,\cdots,n-1}\)，已知 \(f_0=1,f_i=\sum_{j=1}^i f_{i-j}g_j\)，求 \(f_n\)。

显然分治计算贡献的形式是可以的。考虑 \(\texttt{cdq}\) 分治，划分左右两部，首先求出 \(f_{\{L\}}\)，然后 \(\texttt{FFT}\) 求出 \(f_{\{L\}}\to f_{\{R\}}\) 的贡献，然后考虑 \(f_{\{R\}}\) 内部的贡献。

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

然而可能有更简单的做法。

考虑将序列转为多项式 \(f\to f(x),g\to g(x)\)。容易发现如果要满足多项式的乘法，那么 \(f(x)\) 应该是无限项的。特别地，\([x^0]g(x)=0\)。

显然以上形式肯定会有 \(f(x)=f(x)g(x)+\Delta\)，其中 \(\Delta\) 是一个有限的多项式。

其实可以得出 \(\Delta = f_0=1\)。

因此：

\[f(x)=\frac{1}{1-g(x)} \]

得到了一个 \(O(n\log_2 n)\) 的做法。