多项式全家桶

整理各种多项式函数，其实只不过是一部分而已，写出来方便复习吧(虽然似乎已经快要退役了呢) ## 快速傅里叶变换 fft

快速数论变换 ntt

离散傅里叶变换 `dft(poly& f)`

类似fft，需要在 \(\mathbb{F}_P\)下找到合适的\(n\)次单位根\(w_n\) 根据数论知识，\(n\)次单位根存在当且仅当 \(n \mid \phi(P)\), 即 \(n | (P - 1)\) 而由于我们需要做蝶形变换，因此 \(n\) 是 \(2\) 的高次幂，因此我们希望 \(P\) 可以写成： \[P = q\cdot 2^k + 1\] 这样的形式，且\(k\)尽可能大，以方便我们可以对更长的数组做ntt. \(998244353\)是一个合适的数，其中一个原根 \(g = 3\) 知道原根后，\(n\)次单位根\(\omega_n\)是易求的： \[\omega_n = g^{\frac{P - 1}{n}}\] 式子右边很容易验证符合 \(w_n\) 的性质. 离散傅里叶变换公式： \[F_k = f(\omega^k) = \sum_{j = 0}^{n - 1} w^{jk}f_j \] 于是仿照fft做蝶形变换即可

离散傅里叶逆变换 `idft(poly& f)`

ntt实现的idft可以使用一点技巧 fft中，实现逆变换是使用\(\omega_{n}^{-1}\),这里也可以这么做，不过也可以利用数论知识换一下形式: \[\begin{align*} f_k = F(\omega^{-k}) &= \sum_{j = 0}^{n - 1} w^{-jk}F_j \\ &= \sum_{j = 0}^{n - 1} w^{(n-j)k}F_j \\ &= \sum_{j = 0}^{n - 1} w^{jk}F_{n - j \bmod n} \end{align*}\] 因此我们可以对于\(j \in [1, n - 1]\), 我们翻转这一段，而\(n - 0 = 0 \pmod n\)，不变，然后做一次dft. 于是写出来一个很像 1-base 但是完全不是这么回事的代码：

void idft(vector<int>& f)
{
    int n = f.size();
    int invn = inv(n);
    reverse(f.begin() + 1, f.end());
    dft(f);
    for(int i = 0;i < n;i++){
        f[i] = 1ll * f[i] * invn % P;
    }
    return;
}

原理

为什么fft可以把卷积变成点积？甚至这一策略可以平移到\(\mathbb{F}_P\)？其实可以只关注我们利用到的性质，在最小条件下进行推理。给定一个交换环\(R\)， \(\omega \in R\), \(\operatorname{ord}(\omega) = n\) 定义\(R^n\)上的dft(没有归一化)： \[\hat{a}(k) = \sum_{j = 0}^{n - 1}a_j w^{jk}, k = 0,1,\dots,n-1\] 定义\(R^n\)上的循环卷积： \[(a * b)_t = \sum_{j = 0}^{n - 1} a_j b_{t - j \bmod n}\]

于是可以证明： \[\widehat{a*b} = \hat{a} * \hat{b}\]

给出一个证明：令\(c = a * b\)有： \[\begin{align*} \hat{c}_k &= \sum_{j = 0}^{n - 1} c_j w^{jk} \\ &= \sum_{j = 0}^{n - 1} \sum_{t = 0}^{n - 1} (a_t b_{j - t \bmod n}) w^{jk}\\ &= \sum_{t = 0} a_t \sum_{j = 0}^{n - 1} b_{j - t \bmod n} w^{jk} \\ &= \sum_{t = 0} a_t \sum_{s = j - t \bmod n, s = 0}^{n - 1} b_{s} w^{(s + t)k} \\ &= \sum_{t = 0} a_t w^{tk} \sum_{s = 0}b_{s} w^{sk} \\ &= \hat{a}_k \hat{b}_k \end{align*}\]

所以这个策略在 \(C\), \(\mathbb{F}_P\)上均有用有一个平时注意不到的点是：我们实际上实现的是循环卷积，而非普通卷积，这个性质在一些特殊的场合可以发挥作用。也就是说，我们实现的多项式乘法不是在\(\mathbb{F}_P^n\)上，而是在\(\mathbb{F}_P[x]/(x^n - 1)\)这个商环上。本质上，dft是商环\(R[x]/(x^n - 1)\)到直积\(R^n\)的同构。

截断 `truck(int m)`

从第\(m\)位开始截断，形式化即\(\mod x^m\)

右移 `shift(int m)`

移位函数，在左边插入若干个\(0\)，或者截断一些位形式化即 \(x^m f\) 去掉负数项

翻转 `rev()`

从最高位开始翻转多项式形式化来说，相当于\(f^R = x^{\operatorname{deg}(f)} f(\frac{1}{x})\)

形式导数 `deriv()`

求形式导数

形式不定积分 `integr()`

求形式积分

多项式乘法逆 `inv(int m)`

使用多项式牛顿迭代： \[\begin{align*} Af &= 1\\ A - \frac{1}{f} &= 0\\ \rightarrow f_{2n} = f_n - \frac{A - f_n^{-1}}{-f_n^{-2}} &= f_n(2 + Af_n) \end{align*}\] 从\(f_0 = A_0^{-1}\)开始迭代即可因此要求常数位非\(0\) 另外一种解释是可以看作\(\mod x^\infty\) 下求逆，那么存在逆元当且仅当\(\gcd(f, x^\infty) = 1\),即\(f\)常数项不为\(0\).

多项式对数函数 `log(int m)`

\[\begin{align*} f &= \ln A(x) \\ f'&= \frac{1}{A(x)} \\ f &= \int \frac{1}{A(x)} \end{align*}\] 求逆再积分回去即可根据定义，常数项必须是1

多项式指数函数 `exp(int m)`

有两种方法，\(O(n\log n)\) 的牛顿迭代和 \(O(n \log^2 n)\) 的半在线卷积实际上半在线卷积要比牛顿迭代快不少

牛顿迭代

根据初始条件写出迭代式： \[\begin{align*} f &= e^{A}\\ A - \ln f &= 0\\ \rightarrow f_{2n} = f_n - \frac{A - \ln f_n}{-f_n^{-1}} &= f_n(1 + A - \ln f_n) \end{align*}\] 从\(f_0 = 1\)开始迭代即可根据主定理，时间复杂度为\(O(n\log n)\)

半在线卷积

\[\begin{align*} F &= e^{G}\\ \ln F &= G\\ \frac{F'}{F} &= G'\\ F' &= FG' \\ \rightarrow nf_{n} &= \sum_{k = 0}^{n - 1} kg_k f_{n - k} \end{align*}\] (对求导算子对系数的影响要熟练) 用半在线卷积实现即可

多项式快速幂 `pow(int n, int k)`

如果常数项为\(1\)，直接\(\ln + \exp\)即可处理常数项非\(1\)，那么不能直接取对数函数首先，先去掉前面为\(0\)的项，相当于提一个因子\(x^i\) 其次，如果处理后的常数项非\(1\)，那么再提一个常数因子\(v\)即可求解多项式快速幂后，乘上 \(v^k\)，补上 \(ki\) 个 \(0\) 即可

在\(k\)相当大的时候需要注意的一些点：实际上，可以证明： \[f(x)^k \equiv f(x)^{k \bmod P} \mod x^P\] 这个可以展开\(\exp\)系数得到这个结论然而，在我们上面的处理中，有一步是乘上 \(v^k\),显然这里的 \(k\) 必须模 \(\phi(P) = P - 1\) 注意这两个地方的模数是不一样的！ 此外，我们实际上要补\(ki\) 而非 \((k \bmod P)i\) 个 \(0\)，所以当 \(k \equiv 0 \pmod P\)时，要特判一下，避免没有补应该补的 \(0\) 导致错误。

多项式开根 `sqrt(int m)`

使用多项式牛顿迭代 \[\begin{align*} f &= \sqrt{A}\\ A - f^2 &= 0\\ \rightarrow f_{2n} = f_n - \frac{A - f_n^2}{2f_n} &= \frac{1}{2 f_n}(A + f_n) \end{align*}\] 从\(f_0 = \sqrt{A_0}\)开始迭代即可这部分时间复杂度 \(O(n \log n)\) 这里要求 \(A_0\) 是\(\mod P\) 的二次剩余可以用Cipolla算法求二次剩余不过也可以用慢一些的算法，毕竟主要开销还是牛顿迭代. 选一个原根\(g\)，用bsgs算出阶 \(ord\)，那么\(g^{\frac{ord}{2}}\)就是一个平方根. 这样的话，时间复杂度是\(O(\sqrt P)\)

多项式除法 `/` & 取余 `%`

除朴素除法外，其实有更快计算多项式带余除法的方法. 设\[A(x) = Q(x)B(x) + R(x)\] 其中 \(\operatorname{deg}(A) = n\), \(\operatorname{deg}(B) = m\) 可得：\(\operatorname{deg}(Q) = n - m\), \(\operatorname{deg}(R) < m\) 代入\(\frac{1}{x}\) \[\begin{align*} A(\frac{1}{x}) &= Q(\frac{1}{x})B(\frac{1}{x}) + R(\frac{1}{x}) \\ x^n A(\frac{1}{x}) &= x^{n - m}Q(\frac{1}{x}) x^m B(\frac{1}{x}) + x^n R(\frac{1}{x})\\ A^R(x) &= Q^R(x) B^R(x) + x^{n - \operatorname{deg}(R)} R^R(x) \end{align*}\] 由于\(n - \operatorname{deg}(R) > n - m\) 所以等式两边的小于等于 \(n - m\) 次项与 \(R^R(x)\)无关，因为它乘上了\(x^{n - \operatorname{deg}(R)}\)这个系数即 \[A^R(x) = Q^R(x) B^R(x) \mod x^{n - m + 1}\\ Q^R(x) = \frac{A^R(x)}{B^R(x)} \mod x^{n - m + 1}\] 于是通过一次多项式求逆和多项式乘法就可以算出\(Q(x)\)，而知道\(Q(x)\)后，\(R(x)\)就很容易计算了。

多项式多点求值 `eval(vector<int>x)`

其实就是使用分治去求解，不过可以通过一些科技去显著降低常数 ### 朴素方法求解\(f(x_0),f(x_1) \dots f(x_{n - 1})\),设 \(m = \lfloor\frac{n}{2}\rfloor\), \(M_{l,r} = \prod_{i = l}^r (x-x_i)\) 将\(f(x)\) 拆成两部分: \(f(x) = M_{0,m}(x)P(x) + R(x)\) 显然左边这部分对于\(x = x_0, x_1,\dots,x_m\)的贡献均为 \(0\) 于是我们可以递归下去，对 \(R\) 求解。这里需要一次多项式取模，复杂度\(O(n \log n)\)，不过常数偏大。而对于\(x = x_{m + 1}, \dots , x_{n - 1}\)，我们模上 \(M_{m + 1, n - 1}\) 同样递归求解即可。构建 \(M\) 我们使用分治ntt,时间复杂度是 \(O(n \log^2 n)\) 于是最后的时间复杂度 \(T(n) = O(n \log n) + 2 T(\frac{n}{2}) \Rightarrow T(n) = O(n \log^2 n)\) 实际上，很容易注意到一个点: \[f(x_i) = f(x) \bmod (x - x_i)\] 我们上述做法不过是在分治地取模而已，类似于分治ntt(虽然这个是自顶向下的)

Laurent 级数 / 转置原理(特勒根原理)

这里给出一个形式洛朗级数的式子 \[f(v) = [x^0]\frac{f(x^{-1})}{1 - vx}\] 这是好验证的，展开分母，\(x^0\)的系数显然是： \[\sum_{i = 0}^{\infty} f_i v^i = f(v)\] 原理推导完毕，我们希望对于每一个 \(v = x_i\) 计算出上式。于是还是考虑分治：设\(M_{l, r} = \prod_{i = l} ^ r (1 - x_i x)\) 我们从\(\frac{f(x^{-1})}{M_{0, n - 1}}\)开始递归，从区间 \([l, r]\) 求解区间 \([l, m]\) 时,乘上\(M_{m+1,r}\)即可。此时我们每次递归下去仅需要一次多项式乘法，相比多项式取模常数大大减小。

不过需要注意的是：这里的乘法不是一般的矩阵乘法！因为我们维护的是一个以\(x^{-1}\)为变量的多项式，所以这里实际上是差卷积，也可以从转置原理的角度理解为转置乘法。因为我们希望最终得到\(x^0\)的系数，于是如果这个区间还有\(k\)个点值需要求解，我们将多项式截断到\(x^{-k}\)即可递归过程的时间复杂度是类似的，为\(O(n \log^2 n)\), 加上构建 \(M\) 的时间复杂度 \(O(n \log^2 n)\) 和一次多项式求逆 \(O(n \log n)\)，总时间复杂度\(O(n \log^2 n)\)

当然，这个做法可以从朴素做法的转置这个角度来理解，不过可能需要对线性变换的理解到达一定程度。

多项式快速插值 `lagrange(vector<int>x, vector<int>y)`

设\(M(x) = \prod_{i = 0}^n (x - x_i)\) 考虑lagrange插值公式： \[\begin{align*} f(x) &= \sum_{i = 0}^{n} y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}\\ \end{align*}\] 对于右边式子的分母，使用洛必达法则，有: \[\prod_{j \neq i} \frac{1}{x_i - x_j} = \lim_{x\rightarrow x_i} \frac{x - x_i}{\prod_j x - x_j} = \lim_{x\rightarrow x_i} \frac{x - x_i}{M(x)} = \frac{1}{M'(x_i)}\] 于是 \[\begin{align*} f(x) &= \sum_{i = 0}^{n} \frac{y_i}{M'(x_i)} \prod_{j \neq i} x - x_j\\ &= \sum_{i = 0}^{n} v_i \prod_{j \neq i} x - x_j \end{align*}\] 分治ntt算出\(M(x)\), 用多点求值求出\(M'(x_i)\), 进一步求出 \(v_i\) 剩下的问题就是怎么计算上面那个和式还是考虑分治，假设分治到区间\([l, r]\)有结果 \[\begin{align*} f_{l, r} &= \sum_{i = l}^r v_i\prod_{j \neq i, l \leq j \leq r}x - x_j \\ &= f_{l, m} \prod_{j = m + 1}^{r}x - x_j + f_{m + 1, r} \prod_{j = l}^{m}x - x_j \end{align*}\] 前面分治ntt计算\(M(x)\)时预处理好连乘积多项式，按照上式分治递归下去即可。时间复杂度为\(O(n \log^2 n)\)

多项式平移 `taylorShift(Poly f, int c)`

将系数表示的多项式\(f(x)\) 变为 \(f(x + c)\) \[\begin{align*} f(x + c) &= \sum_i f_i (x + c)^i\\ &= \sum_{i = 0}^{n} f_i \sum_{j = 0}^i \binom{i}{j} x^j c^{i - j}\\ &= \sum_{i = 0}^{n} i! f_i \sum_{j = 0}^i \frac{x^j}{j!} \frac{c^{i - j}}{(i - j)!}\\ &= \sum_{j = 0}^{n} \frac{x^j}{j!} \sum_{i = j}^{n} i! f_i \frac{c^{i - j}}{(i - j)!}\\ \end{align*}\] 显然这是一个差卷积的形式，直接是实现即可。时间复杂度\(O(n \log n)\)

多项式连续点值平移 `valueShift(poly f, int c)`

连续点值表示的多项式\(f(0),f(1)\dots,f(n)\)，计算出\(f(c),f(c + 1),\dots,f(c + n)\), \(c > n\) 考虑Lagrange插值 \[\begin{align*} f(x)&= \sum_{i = 0}^{n} f(i) \prod_{j \neq i}\frac{x - j}{i - j}\\ &= \sum_{i = 0}^{n} f(i) \frac{(-1)^{n - i} x^{\underline{n}}}{(x - i)i!(n - i)!}\\ &= x^{\underline{n + 1}}\sum_{i = 0}^{n} \frac{(-1)^{n - i}f(i)}{i!(n - i)!} \frac{1}{(x - i)}\\ \frac{f(x)}{x^{\underline{n + 1}}} &= \sum_{i = 0}^{n} \frac{(-1)^{n - i}f(i)}{i!(n - i)!} \frac{1}{(x - i)}\\ \end{align*}\] 我们得到了一个恒等式，而我们希望将左边代入\(f(x + t)\)时通过右式求值(实际上是通过分式线性插值)。此时右式是一个卷积形式，我们不妨设 \[A(x) = \sum_{i = 0}^n\frac{(-1)^{n - i}f(i)}{i!(n - i)!}x^i\\ B(x) = \sum_{i}\frac{1}{c + i - n}x^i\] 有： \[\begin{align*} [x^{n + t}]A(x)B(x) &= \sum_{i = 0}^{n}\frac{(-1)^{n - i}f(i)}{i!(n - i)!} \frac{1}{c + (n + t - i) - n}\\ &= \sum_{i = 0}^{n}\frac{(-1)^{n - i}f(i)}{i!(n - i)!} \frac{1}{c + t - i}\\ &= \frac{f(c + t)}{(c + t)^{\underline{n + 1}}} \\ f(c + t) &= (c + t)^{\underline{n + 1}}[x^{n + t}]A(x)B(x) \end{align*}\] 需要注意的是\(B\)序列更长一些(\(2n + 1\))。我们通过\(O(n)\)预处理下降幂，时间复杂度由卷积决定，为\(O(n \log n)\)

常系数齐次线性递推

Kitamasa/Fiduccia 算法 `kitamasa(vector<int>a, vector<int>c, i64 n)`

Todo

表示为有理函数 + Bostan–Mori 算法 `bostanMori(poly p, poly q, poly n)` `linearRecurrence(poly a, poly c, i64 n)`

Todo

最小多项式 Berlekamp Massey算法 `berlekampMassey`

Todo

多项式复合 `comp(poly g, int m)`

\(h(x) = f(g(x)) \equiv [y^{0}] \frac{f(y^{-1})}{1 - yg(x)} \pmod {x^m}\) Todo

多项式复合逆 `revert(int m)`

Todo

代码

附上最后的整体代码

#include<bits/stdc++.h>
using i64 = long long;
using ll = long long;
using uint = unsigned int;
using ull = unsigned long long;

using namespace std;

//polys.cpp Create time : 2026.01.13 22:21
constexpr int P = 998244353;
constexpr int G = 3;

int power(int a, int p)
{
    int res = 1;
    while(p){
        if(p & 1)res = 1ll * res * a % P;
        a = 1ll * a * a % P;
        p >>= 1;
    }
    return res;
}

int power(int a, int p, int M)
{
    int res = 1;
    while(p){
        if(p & 1)res = 1ll * res * a % M;
        a = 1ll * a * a % M;
        p >>= 1;
    }
    return res;
}

int inv(int a)
{
    return power(a, P - 2);
}

struct Comb {
    int n;
    std::vector<int> _fac;
    std::vector<int> _invfac;
    std::vector<int> _inv;
     
    Comb() : n{0}, _fac{1}, _invfac{1}, _inv{0} {}
    Comb(int n) : Comb() {
        init(n);
    }
     
    void init(int m) {
        if (m <= n) return;
        _fac.resize(m + 1);
        _invfac.resize(m + 1);
        _inv.resize(m + 1);
         
        for (int i = n + 1; i <= m; i++) {
            _fac[i] = 1ll * _fac[i - 1] * i % P;
        }
        _invfac[m] = power(_fac[m], P - 2);
        for (int i = m; i > n; i--) {
            _invfac[i - 1] = 1ll * _invfac[i] * i % P;
            _inv[i] = 1ll * _invfac[i] * _fac[i - 1] % P;
        }
        n = m;
    }
     
    int fac(int m) {
        if (m > n) init(2 * m);
        return _fac[m];
    }
    int invfac(int m) {
        if (m > n) init(2 * m);
        return _invfac[m];
    }
    int inv(int m) {
        if (m > n) init(2 * m);
        return _inv[m];
    }
    int binom(int n, int m) {
        if (n < m || m < 0) return 0;
        return fac(n) * invfac(m) * invfac(n - m);
    }
} comb;


int bsgs(int a, int b, int p)
{
    a %= p;
    b %= p;
    if(p == 1) return 0;
    if(b == 1) return 0;

    int sq = ceil(sqrt(p));

    unordered_map<int,int>mp;
    int x = 1;
    for(int i = 0;i <= sq;i++){
        mp[1ll * x * b % p] = i;
        x = 1ll * x * a % p;
    }

    int asq = power(a, sq, p);
    int y = asq;
    for(int i = 1;i <= sq;i++){
        if(mp.find(y) != mp.end()){
            return i * sq - mp[y];
        }
        y = 1ll * y * asq % p;
    }
    return -1;
}

int Fsqrt(int x)
{
    int p = bsgs(G, x, P);
    return power(G, p / 2);
}

constexpr int Q = 23; //998244353 = K 2^Q + 1
constexpr int MAX_LEN = 22;

array<int, 1 << MAX_LEN> omega()
{
    array<int, 1 << MAX_LEN> res;
    res.fill(1);
    for (int i = 1; i <= (1 << MAX_LEN - 1); i <<= 1) {
        int wi = power(G, (P - 1) / (i * 2));
        for (int j = 1; j < i; j++) {
            res[i + j] = 1ll * res[i + j - 1] * wi % P;
        }
    }
    return res;
}

const array<int, 1 << MAX_LEN> W = omega();

void dft(vector<int>& f)
{
    int n = (int)f.size();

    // DIF: input natural order, output bit-reversed order
    for (int k = n >> 1; k; k >>= 1) {
        for (int i = 0; i < n; i += k << 1) {
            for (int j = 0; j < k; j++) {
                int x = f[i + j];
                int y = f[i + j + k];

                // f[i+j] = x + y
                int s = x + y;
                if (s >= P) s -= P;
                f[i + j] = s;

                // f[i+j+k] = (x - y) * W[k+j]
                int d = x - y;
                if (d < 0) d += P;
                f[i + j + k] = 1ll * d * W[k + j] % P;
            }
        }
    }
    return;
}

void idft(vector<int>& f)
{
    int n = (int)f.size();
    int invn = inv(n);

    // IDIT: take bit-reversed order (from DIF), output natural order
    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 x = f[i + j];
                int y = 1ll * f[i + j + k] * W[k + j] % P;

                int s = x + y;
                if (s >= P) s -= P;

                int d = x - y;
                if (d < 0) d += P;

                f[i + j] = s;
                f[i + j + k] = d;
            }
        }
    }

    // scale by 1/n
    for (int i = 0; i < n; i++) {
        f[i] = 1ll * f[i] * invn % P;
    }

    // reverse except f[0]  (invNTT trick with +roots)
    reverse(f.begin() + 1, f.end());
    return;
}

class poly : public vector<int>
{
public:
    poly() : vector<int>(){}

    template<class F>
    explicit poly(int n, F f) : vector<int>(n)
    {
        for(int i = 0;i < n;i++){
            (*this)[i] = f(i);
        }
    }
    template<class InputIt, class = std::_RequireInputIter<InputIt>>
    explicit poly(InputIt first, InputIt last) : vector<int>(first, last) {}
    
    explicit poly(int n, int val) : vector<int>(n, val){}
    explicit poly(int n) : vector<int>(n){}
    explicit poly(const vector<int> &a) : vector<int>(a) {}
    poly(const std::initializer_list<int> &a) : vector<int>(a){}

    friend poly operator+(const poly &a, const poly &b) {
        poly res(std::max(a.size(), b.size()));
        for (int i = 0; i < a.size(); i++) {
            res[i] = a[i];
        }
        for (int i = 0; i < b.size(); i++) {
            res[i] = res[i] + b[i];
            if(res[i] >= P)res[i] -= P;
        }
        return res;
    }
    friend poly operator-(const poly &a, const poly &b) {
        poly res(std::max(a.size(), b.size()));
        for (int i = 0; i < a.size(); i++) {
            res[i] = a[i];
        }
        for (int i = 0; i < b.size(); i++) {
            res[i] = res[i] - b[i];
            if(res[i] < 0)res[i] += P;
        }
        return res;
    }
    friend poly operator-(const poly &a) {
        std::vector<int> res(a.size());
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] ? P - a[i] : 0;
        }
        return poly(res);
    }
    friend poly operator*(int a, poly b) {
        a = (a % P + P) % P;
        for (int i = 0; i < int(b.size()); i++) {
            b[i] = 1ll * b[i] * a % P;
        }
        return b;
    }
    friend poly operator*(poly a, int b) {
        b = (b % P + P) % P;
        for (int i = 0; i < int(a.size()); i++) {
            a[i] = 1ll * a[i] * b % P;
        }
        return a;
    }
    friend poly operator/(poly a, int b) {
        b = (b % P + P) % P;
        int invb = power(b, P - 2);
        for (int i = 0; i < int(a.size()); i++) {
            a[i] = 1ll * a[i] * invb % P;
        }
        return a;
    }

    friend poly operator/(const poly &a, const poly &b) {
        if(a.size() < b.size())return poly();
        int k = a.size() - b.size() + 1;
        return (a.rev().trunc(k) * b.rev().inv(k).trunc(k)).trunc(k).rev();
    }
    friend poly operator%(const poly &a, const poly &b) {
        return (a - a / b * b).trunc(min(a.size(), b.size() - 1));
    }

    poly &operator+=(poly b) {
        return (*this) = (*this) + b;
    }
    poly &operator-=(poly b) {
        return (*this) = (*this) - b;
    }
    poly &operator*=(poly b) {
        return (*this) = (*this) * b;
    }
    poly &operator*=(int b) {
        return (*this) = (*this) * b;
    }
    poly &operator/=(int b) {
        return (*this) = (*this) / b;
    }

    friend poly operator*(poly a, poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return poly();
        }
        if (a.size() < b.size()) {
            std::swap(a, b);
        }
        int n = 1, tot = a.size() + b.size() - 1;
        while (n < tot) {
            n *= 2;
        }
        if (((P - 1) & (n - 1)) != 0 || b.size() < 128) {
            poly c(a.size() + b.size() - 1);
            for (int i = 0; i < a.size(); i++) {
                for (int j = 0; j < b.size(); j++) {
                    c[i + j] = (c[i + j] + 1ll * a[i] * b[j]) % P;
                }
            }
            return c;
        }
        a.resize(n);
        b.resize(n);
        dft(a);
        dft(b);
        for (int i = 0; i < n; ++i) {
            a[i] = 1ll * a[i] * b[i] % P; 
        }
        idft(a);
        a.resize(tot);
        return a;
    }

    poly shift(int k) const {
        if (k >= 0) {
            auto b = *this;
            b.insert(b.begin(), k, 0);
            return b;
        } else if (this->size() <= -k) {
            return poly();
        } else {
            return poly(this->begin() + (-k), this->end());
        }
    }

    poly trunc(int k) const {
        poly f = *this;
        f.resize(k);
        return f;
    }

    poly deriv() const {
        if (this->empty()) {
            return poly();
        }
        poly res(this->size() - 1);
        for (int i = 0; i < this->size() - 1; ++i) {
            res[i] = 1ll * (i + 1) * (*this)[i + 1] % P;
        }
        return res;
    }

    poly integr() const {
        poly res(this->size() + 1);
        comb.init(this->size() + 1);
        for (int i = 0; i < this->size(); ++i) {
            res[i + 1] = 1ll * (*this)[i] * comb.inv(i + 1) % P;
        }
        return res;
    }

    poly rev() const{
        if(this->empty()){
            return poly();
        }
        poly res = *this;
        while(!res.empty() && res.back() == 0)res.pop_back();
        reverse(res.begin(),res.end());
        return res;
    }

    poly inv(int m) const {
        poly x{power((*this)[0], P - 2)};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (poly{2} - trunc(k) * x)).trunc(k);
        }
        return x.trunc(m);
    }

    poly log(int m) const {
        return (deriv() * inv(m)).integr().trunc(m);
    }

    poly exp(int m) const {
        poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (poly{1} - x.log(k) + trunc(k))).trunc(k);
        }
        return x.trunc(m);
    }

    poly exp2(int m) const {
        poly f(m), h = deriv();
        h.resize(m);
        f[0] = 1;
        // fn = 1/n(\sum h_k f_{n - 1 - k})
        auto cdq = [&](auto&&self, int l ,int r)->void
        {
            if(l == r){
                if(l == 0)return;
                f[l] = 1ll * f[l] * comb.inv(l) % P;
                return;
            }
            int m = l + r >> 1;
            self(self, l, m);
            poly a(m - l + 1),b(r - l);
            for(int i = 0;i <= m - l;i++){
                a[i] = f[i + l];
            }
            for(int i = 0;i < r - l;i++){
                b[i] = h[i];
            }
            auto c = a * b;
            for(int i = 0;i < r - m;i++){
                f[m + i + 1] = f[m + i + 1] + c[i + m - l];
                if(f[m + i + 1] >= P)f[m + i + 1] -= P;
            }
            self(self, m + 1, r);
            return;
        };
        cdq(cdq, 0, m - 1);
        return f;
    }

    poly pow(int k, int m) const {
        int i = 0;
        while (i < this->size() && (*this)[i] == 0) {
            i++;
        }
        if (i == this->size() || 1LL * i * k >= m) {
            return poly(m);
        }
        int v = (*this)[i];
        auto f = shift(-i) * power(v, P - 2);
        return ((f.log(m - i * k) * k).exp2(m - i * k).shift(i * k)) * power(v, k);
    }

    poly pow(int k, int kmodphiP,int big, int m) const {
        int i = 0;
        while (i < this->size() && (*this)[i] == 0) {
            i++;
        }
        if (i == this->size() || 1LL * i * k >= m) {
            return poly(m);
        }
        if(big && i > 0)return poly(m);//x^{big number} f(x)^{k}
        int v = (*this)[i];
        auto f = shift(-i) * power(v, P - 2);
        return ((f.log(m - i * k) * k).exp2(m - i * k).shift(i * k)) * power(v, kmodphiP);
    }

    poly sqrt(int m) const {
        poly x{Fsqrt((*this)[0])};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (trunc(k) * x.inv(k)).trunc(k)) * power(2, P - 2);
        }
        return x.trunc(m);
    }

    poly mulT(poly b) const {
        if (b.size() == 0) {
            return poly();
        }
        int n = b.size();
        std::reverse(b.begin(), b.end());
        return ((*this) * b).shift(-(n - 1));
    }

    vector<int> eval(vector<int> x)const 
    {
        if(this->size() == 0){
            return vector<int>(x.size(), 0);
        }
        if(x.size() == 0){
            return vector<int>(0);
        }
        int n = max(x.size(), this->size()); // !!!
        x.resize(n);
        for(int i = 0;i < n;i++){ //norm P - x[i]
            if(x[i] <= 0)x[i] += P;
            if(x[i] > P)x[i] -= P;
        }
        vector<poly>m(n * 4);//segtree build Q(x)
        auto build = [&](auto&&self, int l, int r, int p)->void
        {
            if(l == r){
                m[p] = poly{1, P - x[l]};
                return;
            }
            int mid = l + r >> 1;
            self(self, l, mid , p * 2);
            self(self, mid + 1, r, p * 2 + 1);
            m[p] = m[p * 2] * m[p * 2 + 1];
            return;
        };
        build(build, 0, n - 1, 1);
        vector<int>ans(n);
        auto work = [&](auto&&self,const poly &f, int l, int r, int p)->void
        {
            if(l == r){
                ans[l] = f[0];
                return;
            }
            int mid = l + r >> 1;
            self(self, f.mulT(m[p * 2 + 1]).trunc(mid - l + 1), l, mid, p * 2);
            self(self, f.mulT(m[p * 2]).trunc(r - mid), mid + 1, r, p * 2 + 1);
            return;
        };
        work(work, mulT(m[1].inv(n)).trunc(n), 0, n - 1, 1);
        return ans;
    }

};

poly lagrange(vector<int>x, vector<int>y)
{
    assert(x.size() == y.size());
    if(x.size() == 0)return poly();
    int n = x.size();
    for(int i = 0;i < n;i++){ //norm P - x[i]
        if(x[i] <= 0)x[i] += P;
        if(x[i] > P)x[i] -= P;
    }
    for(int i = 0;i < n;i++){ //norm P - x[i]
        if(y[i] <= 0)y[i] += P;
        if(y[i] > P)y[i] -= P;
    }
    vector<poly>m(n * 4);//segtree build Q(x)
    auto build = [&](auto&&self, int l, int r, int p)->void
    {
        if(l == r){
            m[p] = poly{P - x[l], 1};
            return;
        }
        int mid = l + r >> 1;
        self(self, l, mid , p * 2);
        self(self, mid + 1, r, p * 2 + 1);
        m[p] = m[p * 2] * m[p * 2 + 1];
        return;
    };
    build(build, 0, n - 1, 1);
    auto dM = m[1].deriv();
    auto v = dM.eval(x);
    for(int i = 0;i < n;i++){
        v[i] = 1ll * power(v[i], P - 2) * y[i] % P;
    }

    auto work = [&](auto&&self, int l, int r, int p)->poly
    {
        if(l == r){
            return poly{v[l]};
        }
        int mid = l + r >> 1;
        return m[p * 2] * self(self, mid + 1, r, p * 2 + 1) + m[p * 2 + 1] * self(self, l, mid, p * 2);
    };
    return work(work, 0, n - 1, 1);
}

vector<int> valueShift(const vector<int>&v, int c)
{
    int n = v.size() - 1;
    poly a(n + 1),b(n * 2 + 1);
    //fac, invfac
    vector<int>fac(n + 1), invfac(n + 1);
    vector<int>cfac(n * 2 + 2), invcfac(n + 1); // base c - n - 1
    int cbase = c - n - 1;
    cfac[0] = fac[0] = 1;
    for(int i = 1;i <= n;i++){
        fac[i] = 1ll * fac[i - 1] * i % P;
    }
    invfac[n] = inv(fac[n]);
    for(int i = n - 1;i >= 0;i--){
        invfac[i] = 1ll * invfac[i + 1] * (i + 1) % P;
    }

    for(int i = 1;i <= n * 2 + 1;i++){
        cfac[i] = 1ll * cfac[i - 1] * (i + cbase) % P;
    }
    invcfac[n] = inv(cfac[n]);
    for(int i = n - 1;i >= 0;i--){
        invcfac[i] = 1ll * invcfac[i + 1] * (i + 1 + cbase) % P;
    }

    for(int i = 0;i <= n;i++){
        a[i] = 1ll * v[i] * invfac[i] % P * invfac[n - i] % P;
        if((n - i) % 2 == 1)a[i] = a[i] ? P - a[i] : 0;
    }
    for(int i = 0;i <= n * 2;i++){
        b[i] = inv(c - n + i);
    }
    auto h = a * b;
    vector<int>vc(n + 1);
    for(int i = 0;i <= n;i++){
        vc[i] = 1ll * h[i + n] * cfac[n + i + 1] % P * invcfac[i] % P;
    }
    return vc;
}

poly taylorShift(poly& f, int c)
{
    int n = f.size();
    poly a(n), b(n * 2 - 1);
    vector<int>fac(n + 1), invfac(n + 1);
    fac[0] = 1;
    for(int i = 1;i <= n;i++){
        fac[i] = 1ll * fac[i - 1] * i % P;
    }
    invfac[n] = inv(fac[n]);
    for(int i = n - 1;i >= 0;i--){
        invfac[i] = 1ll * invfac[i + 1] * (i + 1) % P;
    }
    for(int i = 0;i < n;i++){
        a[i] = 1ll * f[i] * fac[i] % P;
    }
    int cp = 1;
    for(int i = 0;i < n;i++){
        b[i] = 1ll * cp * invfac[i] % P;
        cp = 1ll * cp * c % P;
    }
    auto g = a.mulT(b);
    poly res(n);
    for(int i = 0;i < n;i++){
        res[i] = 1ll * g[i] * invfac[i] % P;
    }
    return res;
}

int bostanMori(poly& p, poly& q, i64 n) // calc [x^n] P(x)/Q(x)
{
    int k = q.size() - 1;
    while (n > 0) {
        poly nq = q;
        for (int i = 1; i <= k; i += 2) {
            nq[i] = -nq[i];
            if(nq[i] < 0)nq[i] += P;
        }
        auto np = p * nq;
        nq = q * nq;
        for (int i = 0; i < k; i++) {
            p[i] = np[i * 2 + n % 2];
        }
        for (int i = 0; i <= k; i++) {
            q[i] = nq[i * 2];
        }
        n /= 2;
    }
    return 1ll * p[0] * inv(q[0]) % P;
}

int linearRecurrence(poly a, poly c, i64 n) //a_0 ... a_{k - 1} , c_1 ... c_k
{
    assert(a.size() + 1 == c.size());
    assert(c[0] == 0 || c[0] == P - 1);
    c[0] = -1;
    int k = a.size();
    poly p = (a * (-c)).trunc(k), q = -c;
    return bostanMori(p, q, n);
}

poly berlekampMassey(const vector<int>& s) {
    if (s.empty()) {
        return poly{P - 1};
    }

    poly C{1};   // c[0] = 1
    poly B{1};   
    int L = 0;   
    int m = 1;   // dist from last update
    int b = 1;   // last discrepancy

    for (int n = 0; n < (int)s.size(); n++) {
        int d = 0;
        for (int i = 0; i <= L; i++) {
            if (n - i < 0) break;
            d = (d + 1ll * C[i] * s[n - i]) % P;
        }
        if (d == 0) {
            m++;
            continue;
        }

        poly T = C; 

        // coef = d / b
        int coef = 1LL * d * inv(b) % P;

        // C = C - coef * x^m * B
        if ((int)C.size() < (int)B.size() + m) C.resize(B.size() + m, 0);
        for (int i = 0; i < (int)B.size(); i++) {
            C[i + m] = C[i + m] - 1ll * coef * B[i] % P;
            if(C[i + m] < 0)C[i + m] += P;
        }

        if (2 * L <= n) {
            L = n + 1 - L;
            B = T;
            b = d;
            m = 1;
        } else {
            m++;
        }
    }
    return -C;
}

namespace compositon{
    using biPoly = vector<poly>; // sum_{j>=0} A[j](x) * y^j

    // (A*B) mod x^nx, y^my
    // packing: x^i y^j -> t^{i + j*block}
    biPoly biMul(const biPoly& A, const biPoly& B, int nx, int my) {
        int Ay = min<int>(A.size(), my);
        int By = min<int>(B.size(), my);
        if (Ay == 0 || By == 0) return {};
        int Cy = min<int>(Ay + By - 1, my);

        int block = 2 * nx - 1;           
        poly PA(Ay * block), PB(By * block);

        for (int j = 0; j < Ay; j++) {
            int lim = min<int>(nx, A[j].size());
            for (int i = 0; i < lim; i++) PA[j * block + i] = A[j][i];
        }
        for (int j = 0; j < By; j++) {
            int lim = min<int>(nx, B[j].size());
            for (int i = 0; i < lim; i++) PB[j * block + i] = B[j][i];
        }

        poly PC = PA * PB;

        biPoly C(Cy, poly(nx));
        for (int j = 0; j < Cy; j++) {
            int base = j * block;
            for (int i = 0; i < nx; i++) {
                int idx = base + i;
                if (idx < (int)PC.size()) C[j][i] = PC[idx];
            }
        }
        return C;
    }

    // Q(-x,y)
    biPoly neg_x(const biPoly& Q, int nx) {
        biPoly R = Q;
        for (auto& fx : R) {
            if ((int)fx.size() < nx) fx.resize(nx);
            for (int i = 1; i < nx; i += 2){
                fx[i] = -fx[i];
                if(fx[i] < 0)fx[i] += P;
            }
        }
        return R;
    }

    // V(x^2,y) = A(x,y)
    // 输入 A mod x^nx，输出 V mod x^{ceil(nx/2)}
    biPoly even_x(const biPoly& A, int nx) {
        int nx2 = (nx + 1) / 2;
        biPoly V(A.size(), poly(nx2));
        for (int j = 0; j < (int)A.size(); j++) {
            int lim = min<int>(nx, A[j].size());
            for (int i = 0; i < nx2; i++) {
                int k = 2 * i;
                if (k < lim) V[j][i] = A[j][k];
            }
        }
        return V;
    }

    // W(x^2,y) mod x^nx
    biPoly subst_x2(const biPoly& W, int nx) {
        biPoly R(W.size(), poly(nx));
        for (int j = 0; j < (int)W.size(); j++) {
            for (int i = 0; i < (int)W[j].size(); i++) {
                int k = 2 * i;
                if (k < nx) R[j][k] = W[j][i];
            }
        }
        return R;
    }

    // 取 y 的 [l, r) 段并左移：返回 size=r-l
    biPoly slice_y(const biPoly& A, int l, int r, int nx) {
        int len = max(0, r - l);
        biPoly R(len, poly(nx));
        for (int i = 0; i < len; i++) {
            int j = l + i;
            if (0 <= j && j < (int)A.size()) {
                R[i] = A[j].trunc(nx);
            }
        }
        return R;
    }

    // 忠实版 Comp：返回 F_{d,m}( P(y) / Q(x,y) ) mod x^nx
    // P 是 y 多项式（长度 m）；Q 是 biPoly（y>=0）
    // 返回 biPoly，size = m-d，且每项是 x 多项式（长度 nx）
    biPoly Comp(const poly& P, const biPoly& Q,int nx, int d, int m) {
        // base: x 只剩常数项
        if (nx == 1) {
            // Q0(y) = Q(0,y)
            poly Q0(m);
            int limy = min<int>(m, Q.size());
            for (int j = 0; j < limy; j++) {
                if (!Q[j].empty()) Q0[j] = Q[j][0];
            }
            // C(y) = P(y) / Q0(y) mod y^m
            poly C = (P * Q0.inv(m)).trunc(m);

            // 输出 y 的 [d,m)
            biPoly res(m - d, poly(1));
            for (int i = d; i < m; i++) res[i - d][0] = C[i];
            return res;
        }

        int degyQ = min(m - 1, (int)Q.size() - 1);
        int e = max(0, d - degyQ);

        biPoly Qneg = neg_x(Q, nx);

        // A = Q * Q(-x)  (mod x^nx, y^m)
        biPoly A = biMul(Q, Qneg, nx, m);

        // V(x,y) = A_even(x,y) with x := x^2 collapsed, so xlen halves
        biPoly V = even_x(A, nx);
        int nx2 = (nx + 1) / 2;

        // W = Comp(nx2, e, m, P, V)   (y窗口从 e 开始)
        biPoly W = Comp(P, V, nx2, e, m);

        // lift: W(x^2,y)
        biPoly W2 = subst_x2(W, nx);

        // B = W(x^2,y) * Q(-x,y) 只需要 y < (m-e)
        biPoly B = biMul(W2, Qneg, nx, m - e);

        // 返回 y 的 [d-e, m-e)
        return slice_y(B, d - e, m - e, nx);
    }

    poly comp(poly f, poly g, int n) {
        f = f.trunc(n);
        g = g.trunc(n);
        if ((int)g.size() < n) g.resize(n);

        int m = min<int>(n, f.size());
        // P(y) = y^{m-1} f(1/y)  <=> P[i] = f[m-1-i]
        poly P(m);
        for (int i = 0; i < m; i++) P[i] = f[m - 1 - i];

        // Q(x,y) = 1 - y g(x)
        biPoly Q(2, poly(n));
        Q[0][0] = 1;
        Q[1] = (-g).trunc(n);
        if ((int)Q[1].size() < n) Q[1].resize(n);

        // 取 y^{m-1}：即 d=m-1, 输出 size=1
        biPoly R = Comp(P, Q, n, m - 1, m);
        return R[0].trunc(n);
    }

    biPoly biMul_fullX(const biPoly& A, const biPoly& B, int nx_in, int my) {
        int Ay = std::min<int>((int)A.size(), my);
        int By = std::min<int>((int)B.size(), my);
        if (Ay == 0 || By == 0) return {};

        int Cy = std::min<int>(Ay + By - 1, my);
        int nx_out = 2 * nx_in - 1;

        // block must be >= nx_out to avoid x-carry crossing y-blocks
        int block = nx_out;

        poly PA(Ay * block), PB(By * block);

        for (int j = 0; j < Ay; j++) {
            int lim = std::min<int>(nx_in, (int)A[j].size());
            for (int i = 0; i < lim; i++) PA[j * block + i] = A[j][i];
        }
        for (int j = 0; j < By; j++) {
            int lim = std::min<int>(nx_in, (int)B[j].size());
            for (int i = 0; i < lim; i++) PB[j * block + i] = B[j][i];
        }

        poly PC = PA * PB;

        biPoly C(Cy, poly(nx_out));
        for (int j = 0; j < Cy; j++) {
            int base = j * block;
            for (int i = 0; i < nx_out; i++) {
                int idx = base + i;
                if (idx < (int)PC.size()) C[j][i] = PC[idx];
            }
        }
        return C;
    }

    // pick x-parity and downsample: keep coeffs x^{2t+parity} -> new[t]
    biPoly pick_parity_x(const biPoly& A, int parity, int nx_new) {
        biPoly R(A.size(), poly(nx_new));
        for (int j = 0; j < (int)A.size(); j++) {
            int lim = (int)A[j].size();
            for (int t = 0; t < nx_new; t++) {
                int idx = 2 * t + parity;
                if (idx < lim) R[j][t] = A[j][idx];
            }
        }
        return R;
    }

    // ans[k] = [x^n] f(x)^k * g(x), k=0..n
    // require f.size()==g.size()==n+1 and f[0]==0
    poly powerProjection(poly f, poly g) {
        int N = (int)f.size();
        int n = N - 1;
        // 当前 xlen = n+1, 当前 ylen = my
        int nx = n + 1;
        int my = std::min(N, 2); // 初始 Q = 1 - y f, 只需要 y^0..y^1

        biPoly P, Q;
        P.assign(1, g.trunc(nx)); // P(y) = g(x)
        Q.assign(2, poly(nx));    // Q(y) = 1 - y f(x)
        Q[0][0] = 1;
        Q[1] = (-f).trunc(nx);

        // 主循环：x-Bostan–Mori，目标提取 [x^n]
        while (n > 0) {
            // ylen 更新：乘法会让 ylen 至多变成 2*my-1，但我们只关心到 N 项
            int my2 = std::min(N, 2 * my - 1);

            // R = Q(-x,y)
            biPoly R = neg_x(Q, nx);

            // S = P * R , T = Q * R  (x full conv, y truncated)
            biPoly S = biMul_fullX(P, R, nx, my2);
            biPoly T = biMul_fullX(Q, R, nx, my2);

            int parity = (n & 1);
            int n2 = n >> 1;
            int nx2 = n2 + 1;

            // P <- pick parity row from S; Q <- even row from T
            P = pick_parity_x(S, parity, nx2);
            Q = pick_parity_x(T, 0, nx2);

            // 截 ylen 到 my2（pick_parity_x 保留了全部 y 行）
            if ((int)P.size() > my2) P.resize(my2);
            if ((int)Q.size() > my2) Q.resize(my2);

            // shrink n, nx, my
            n = n2;
            nx = nx2;
            my = my2;
        }

        // n == 0 => x 只剩常数项：结果是 y-series = P0(y)/Q0(y)
        // 取各 y^j 的 x^0 系数拼成 y 多项式
        poly Py(my), Qy(my);
        for (int j = 0; j < my; j++) {
            if (j < (int)P.size() && !P[j].empty()) Py[j] = P[j][0];
            if (j < (int)Q.size() && !Q[j].empty()) Qy[j] = Q[j][0];
        }

        poly invQy = Qy.inv(my);
        poly ans = (Py * invQy).trunc(N); // 需要 0..N-1 共 N 项（k=0..n）
        return ans;
    }

    // g = f^{<-1>} mod x^m, require f[0]=0 and f[1]!=0
    poly compInv(poly f, int m) {
        if (m <= 0) return poly();
        if (m == 1) return poly{0};

        f = f.trunc(m);
        // assert(f.size() >= 2 && f[0] == 0 && f[1] != 0);

        int c = f[1];
        // powerProjection 需要 size = (m+1): n = m
        poly fext = f;
        fext.resize(m + 1);          // f[m]=0 不影响 mod x^m 的 inverse
        poly one(m + 1);
        one[0] = 1;

        // hk[k] = [x^m] f(x)^k
        poly hk = powerProjection(fext, one); // length m+1

        // A(x) = (x/g)^m mod x^m  (length m)
        poly A(m);
        for (int k = 1; k <= m; k++) {
            int deg = m - k;
            // coeff = m/k * hk[k]
            int coef = (long long)m * inv(k) % P * hk[k] % P;
            A[deg] = coef;
        }

        // normalize by c^m so that B[0]=1
        int cpm = power(c, m);
        poly B = A / cpm;
        // (optional) assert(B[0] == 1);

        // H = x/g = c * exp( log(B)/m )
        poly L = (B.log(m) / m); // /m uses scalar inverse via your operator/(int)
        poly H = L.exp(m) * c;

        // g = x * inv(H)
        poly invH = H.inv(m);
        poly ginv = invH.shift(1).trunc(m);
        return ginv;
    }
}


void solve()
{
    
}

int main()
{
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    int tt = 1;
    //cin >> tt;
    while(tt--){
        solve();
    }
    return 0;
}

多项式全家桶

快速数论变换 ntt

离散傅里叶变换 dft(poly& f)

离散傅里叶逆变换 idft(poly& f)

原理

截断 truck(int m)

右移 shift(int m)

翻转 rev()

形式导数 deriv()

形式不定积分 integr()

多项式乘法逆 inv(int m)

多项式对数函数 log(int m)

多项式指数函数 exp(int m)

牛顿迭代

半在线卷积

多项式快速幂 pow(int n, int k)

多项式开根 sqrt(int m)

多项式除法 / & 取余 %

多项式多点求值 eval(vector<int>x)

Laurent 级数 / 转置原理(特勒根原理)

多项式快速插值 lagrange(vector<int>x, vector<int>y)

多项式平移 taylorShift(Poly f, int c)

多项式连续点值平移 valueShift(poly f, int c)

常系数齐次线性递推

Kitamasa/Fiduccia 算法 kitamasa(vector<int>a, vector<int>c, i64 n)

表示为有理函数 + Bostan–Mori 算法 bostanMori(poly p, poly q, poly n) linearRecurrence(poly a, poly c, i64 n)

最小多项式 Berlekamp Massey算法 berlekampMassey

多项式复合 comp(poly g, int m)

多项式复合逆 revert(int m)

代码

离散傅里叶变换 `dft(poly& f)`

离散傅里叶逆变换 `idft(poly& f)`

截断 `truck(int m)`

右移 `shift(int m)`

翻转 `rev()`

形式导数 `deriv()`

形式不定积分 `integr()`

多项式乘法逆 `inv(int m)`

多项式对数函数 `log(int m)`

多项式指数函数 `exp(int m)`

多项式快速幂 `pow(int n, int k)`

多项式开根 `sqrt(int m)`

多项式除法 `/` & 取余 `%`

多项式多点求值 `eval(vector<int>x)`

多项式快速插值 `lagrange(vector<int>x, vector<int>y)`

多项式平移 `taylorShift(Poly f, int c)`

多项式连续点值平移 `valueShift(poly f, int c)`

Kitamasa/Fiduccia 算法 `kitamasa(vector<int>a, vector<int>c, i64 n)`

表示为有理函数 + Bostan–Mori 算法 `bostanMori(poly p, poly q, poly n)` `linearRecurrence(poly a, poly c, i64 n)`

最小多项式 Berlekamp Massey算法 `berlekampMassey`

多项式复合 `comp(poly g, int m)`

多项式复合逆 `revert(int m)`