Poly

polynomial

Intro

Degree

对于一个多项式 $f(x)$，称其最高次项的次数为该多项式的 度（degree），也称次数，记作 $\operatorname{deg}{f}$。

Composite

定义 $R[[x]]$ 中元素 $f$ 的乘方为

$$
f^1=f,f^k=f^{k-1}\times f
$$

在此基础上，定义 $R[[x]]$ 中元素 $f,g$ 的复合为

$$
(f\circ g)(x)=f(g(x))=f_0+\sum_{k=1}^{+\infty}f_kg^k(x)
$$

Derivative

尽管一般环甚至未必存在极限，
我们依然可以定义形式幂级数的 形式导数（formal derivative）为
$$
\left(\sum_{k=0}^{+\infty}f_kx^k\right)’=\sum_{k=1}^{+\infty}kf_kx^{k-1}
$$

Inverse

对于形式幂级数 $f$，若 $f_0\not=0$，其 乘法逆元（multiplicative inversion）$f^{-1}$ 为另一形式幂级数，满足

$$
f\times f^{-1}=f^{-1}\times f=1
$$

用形式幂级数乘法定义展开该式，可得 $f^{-1}$ 系数的递推式

$$
f^{-1}_0=\dfrac{1}{f_0},f^{-1}n=\dfrac{-1}{f_0}\sum{k=0}^{n-1}f^{-1}kf{n-k}
$$

Quotient & Remainder

对于多项式 $f(x), g(x)$，存在唯一的 $Q(x), R(x)$ 满足：

$$
\begin{aligned}
f(x) &= Q(x) g(x) + R(x) \
\operatorname{deg}{R} &< \operatorname{deg}{g}
\end{aligned}
$$

当 $\operatorname{deg}{f} \ge \operatorname{deg}{g}$ 时有 $\operatorname{deg}{Q} = \operatorname{deg}{f} - \operatorname{deg}{g}$，否则有 $Q(x) = 0$。我们称 $Q(x)$ 为 $g(x)$ 除 $f(x)$ 的 商（quotient），$R(x)$ 为 $g(x)$ 除 $f(x)$ 的 余数（remainder）。

Modulo

模多项式是多项式环的子环，由多项式环除以同余的等价关系得到。

在上文提到的带余除法中，多项式 $f(x)$ 与它的余式 $R(x)$ 在模多项式 $g(x)$ 的意义下同余。

$$
f(x) \equiv R(x) \pmod{g(x)}
$$

这个同余式也意味着，对于多项式 $g(x)$ 的任意一个根 $x_0$，代入 $f(x)$ 和 $R(x)$ 中，得到的点值相同。即：

$$
f(x_0)=R(x_0)
$$

并且，如果根 $x_0$ 在多项式 $g(x)$ 中的重数是 $k$，即 $(x-x_0)^k$ 整除 $g(x)$，则对任意大于等于 $0$ 小于 $k$ 的整数 $t$，有：

$$
f^{t}(x_0)=R^{t}(x_0)
$$

这里的记号表示 $t$ 阶导数。

模多项式同余可以应用于幂级数。一个无限项的幂级数，可以在模具体的多项式情形下，和一个有限项的多项式同余。例如：

$$
1+x+x^2+x^3+\ldots \equiv 1+x+\ldots+x^{n-1} \pmod{x^n}
$$

显然剩余的所有项都被 $x^n$ 整除，因此模 $x^n$ 的操作等价于「截断」，将无穷项的幂级数截断到前 $n$ 项，直接将更高位的信息丢失。

Fast Fourier Transform (FFT)

// n -> 2^p
template<typename T>
void dft(vector<complex<T>>& a) {
    constexpr T Pi = acos(-1.0);
    int n = a.size();
    vector<int> rev(n);
    for (int i = 0; i < n; i++) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) {
            rev[i] |= n >> 1;
        }
        if (i > rev[i]) {
            swap(a[i], a[rev[i]]);
        }
    }
    for (int i = 1; i < n; i <<= 1) {
        complex<T> wn(cos(Pi / i), sin(Pi / i));
        for (int j = 0; j < n; j += i * 2) {
            complex<T> w(1, 0);
            for (int k = 0; k < i; k++) {
                complex<T> u = a[k + j];
                complex<T> v = a[k + j + i] * w;
                a[k + j] = u + v;
                a[k + j + i] = u - v;
                w *= wn;
            }
        }
    }
}

// n -> 2^p
template<typename T>
void idft(vector<complex<T>>& a) {
    int n = a.size();
    reverse(a.begin() + 1, a.end());
    dft(a);
    for (auto& x : a) {
        x /= n;
    }
}

求多项式乘法

vector<double> mul(vector<double> a, vector<double> b) {
    int n = 1, tot = a.size() + b.size() - 1;
    while (n < tot) {
        n *= 2;
    }
    vector<complex<double>> A(n), B(n);
    for (int i = 0; i < a.size(); i++) {
        A[i] = a[i];
    }
    for (int i = 0; i < b.size(); i++) {
        B[i] = b[i];
    }
    dft(A);
    dft(B);
    for (int i = 0; i < n; i++) {
        A[i] *= B[i];
    }
    idft(A);
    vector<double> ans(tot);
    for (int i = 0; i < tot; i++) {
        ans[i] = A[i].real();
    }
    return ans;
}

Fast number-theoretic transform (FNTT)

常见质数及原根
$$
p = 167772161 = 5 \times 2^{25}+1, g=3
$$

$$
p = 469762049 = 7 \times 2^{26}+1, g=3
$$

$$
p = 754974721 = 3^2 \times 5 \times 2^{24}+1, g=11
$$

$$
p = 998244353 = 7 \times 17 \times 2^{23}+1, g=3
$$

$$
p = 1004535809 = 479 \times 2^{21}+1, g=3
$$

模数是质数的情况：化成 $p * 2^n + 1$ 的形式（$p$ 为质数）找到原根 $g$， $2^n$ 表示计算的最长长度

前置 $MInt$

template<class T>
constexpr T power(T a, long long b) {
    T res = 1;
    for (; b; a *= a, b >>= 1) {
        if (b & 1) {
            res *= a;
        }
    }
    return res;
}

template<int P>
struct MInt {
	unsigned x;
	constexpr MInt() : x(0) {}
    constexpr MInt(long long sig) : x(norm(sig % P)) {}
	constexpr int val() const { 
        return x; 
    }
    // -P <= v < 2P
    constexpr long long norm(long long v) const {
        if (v < 0) v += P;
        if (v >= P) v -= P;
        return v;
    }
    constexpr MInt inv() const {
        return power(*this, P - 2);
    }
	constexpr MInt& operator+=(MInt that) { 
        x = norm(x + that.x);
        return *this; 
    }
	constexpr MInt& operator-=(MInt that) {
        x = norm(x + P - that.x);
        return *this; 
    }
	constexpr MInt& operator*=(MInt that) {
        x = (unsigned long long)x * that.x % P; 
        return *this; 
    }
	constexpr MInt& operator/=(MInt that) { 
        return (*this) *= that.inv(); 
    }
    constexpr friend MInt operator+(MInt a, MInt b) { return a += b; }
    constexpr friend MInt operator-(MInt a, MInt b) { return a -= b; }
    constexpr friend MInt operator*(MInt a, MInt b) { return a *= b; }
    constexpr friend MInt operator/(MInt a, MInt b) { return a /= b; }
    constexpr friend MInt operator-(MInt a) { return MInt() -= a; } // 负号
    constexpr friend bool operator<(MInt a, MInt b) { return a.x < b.x; }
    constexpr friend bool operator<=(MInt a, MInt b) { return a.x <= b.x; }
    constexpr friend bool operator==(MInt a, MInt b) { return a.x == b.x; }
    constexpr friend bool operator!=(MInt a, MInt b) { return a.x != b.x; }
	friend ostream& operator<<(ostream& os, MInt a) { 
        os << a.x; 
        return os; 
    }
    friend istream& operator>>(istream& is, MInt& a) { 
        long long x; 
        is >> x;
        a = MInt(x); 
        return is; 
    }
};

$FNTT$

// P 质数
template<int P>
constexpr MInt<P> findPrimitiveRoot() {
    long long x = P - 1;
    vector<int> factors;
    for (long long i = 2; i * i <= x; i++) {
        if (x % i == 0) {
            factors.push_back(i);
            while (x % i == 0) {
                x /= i;
            }
        }
    }
    if (x > 1) {
        factors.push_back(x);
    }
    for (int i = 1; i < P; i++) {
        if (power(MInt<P>(i), P - 1) == 1) {
            bool flag = true;
            for (auto& f : factors) {
                if (power(MInt<P>(i), (P - 1) / f) == 1) {
                    flag = false;
                    break;
                }
            }
            if (flag) {
                return i;
            }
        }
    }
    assert(false);
}

template<int P>
constexpr MInt<P> primitiveRoot = findPrimitiveRoot<P>();

template<int P>
void dft(vector<MInt<P>>& a) {
    int n = a.size();
    vector<int> rev(n);
    for (int i = 0; i < n; i++) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) {
            rev[i] |= n >> 1;
        }
        if (i > rev[i]) {
            swap(a[i], a[rev[i]]);
        }
    }
    for (int i = 1; i < n; i <<= 1) {
        MInt<P> wn = power(primitiveRoot<P>, (P - 1) / (i * 2));
        for (int j = 0; j < n; j += i * 2) {
            MInt<P> w = 1;
            for (int k = 0; k < i; k++) {
                MInt<P> u = a[k + j];
                MInt<P> v = a[k + j + i] * w;
                a[k + j] = u + v;
                a[k + j + i] = u - v;
                w *= wn;
            }
        }
    }
}

template<int P>
void idft(vector<MInt<P>>& a) {
    int n = a.size();
    reverse(a.begin() + 1, a.end());
    dft(a);
    MInt<P> inv = MInt<P>(n).inv();
    for (int i = 0; i < n; i++) {
        a[i] *= inv;
    }
}

Fast Walsh Transform (FWT)

FWT 和 FFT 的核心思想应该是相同的，都是对数组的变换。我们记对数组 $A$ 进行快速沃尔什变换后得到的结果为 $FWT[A]$。

那么 FWT 核心思想就是：

我们需要一个新序列 $C$，由序列 $A$ 和序列 $B$ 经过某运算规则得到，即 $C = A \cdot B$；

我们先正向得到序列 $FWT[A], FWT[B]$，再根据 $FWT[C]=FWT[A] \cdot FWT[B]$ 在 $O(n)$ 的时间复杂度内求出 $FWT[C]$，其中 $\cdot$ 是序列对应位置相乘；

然后逆向运算得到原序列 $C$。时间复杂度为 $O(n \log{n})$。

在算法竞赛中，FWT 是用于解决对下标进行位运算卷积问题的方法。

公式：$C_{i} = \sum_{i=j \oplus k}A_{j} B_{k}$

（其中 $\oplus$ 是二元位运算中的某一种）

or

$$
C_{i} = \sum\limits_{i = j | k} A_{j} B_{k}
$$

我们按照定义，显然可以构造 $FWT[A]i = A’i = \sum{i=i\cup j}A{j}$，来表示 $j$ 满足二进制中 $1$ 为 $i$ 的子集。

那么有：

$$
\begin{aligned}
FWT[A]i\cdot FWT[B]i&=\left(\sum{i\cup j=i} A_j\right)\left(\sum{i\cup k=i} B_k\right) \
&=\sum_{i\cup j=i}\sum_{i\cup k=i}A_jB_k \
&=\sum_{i\cup(j\cup k)=i}A_jB_k \
&= FWT[C]_i
\end{aligned}
$$
我们令 $A_0$ 表示 $A$ 的前一半，$A_1$ 表示区间的后一半，那么 $A_0$ 就是 A 下标最大值的最高位为 $0$，他的子集就是他本身的子集（因为最高位为 $0$ 了），但是 $A_1$ 的最高位是 $1$，他满足条件的子集不仅仅是他本身，还包最高位为 $0$ 的子集，即

$$
FWT[A] = merge(FWT[A_0], FWT[A_0] + FWT[A_1])
$$

其中 merge 表示像字符串拼接一样把两个数组拼起来，$+$ 就是普通加法，表示对应二进制位相加。

接下来就是反演了，其实反演是很简单的，既然知道了 $A_0$ 的本身的子集是他自己（$A_0 = FWT[A_0]$），$A_1$ 的子集是 $FWT[A_0] + FWT[A_1]$，那就很简单的得出反演的递推式了：

$$
UFWT[A’] = merge(UFWT[A_0’], UFWT[A_1’] - UFWT[A_0’])
$$
$dft$ 系数 $1$，$idft$ 系数 $-1$

template<int P>
void or_dft(vector<MInt<P>>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                a[i + j + k] += a[j + k]; // * 1
            }
        }
    }
}

template<int P>
void or_idft(vector<MInt<P>>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                a[i + j + k] -= a[j + k]; // * -1
            }
        }
    }
}

template<int P>
vector<MInt<P>> or_fwt(vector<MInt<P>> a, vector<MInt<P>> b) {
    int n = 1, tot = max(a.size(), b.size());
    while (n < tot) {
        n *= 2;
    }
    a.resize(n);
    b.resize(n);
    or_dft(a);
    or_dft(b);
    for (int i = 0; i < n; i++) {
        a[i] *= b[i];
    }
    or_idft(a);
    return a;
}

and

$$
C_{i} = \sum\limits_{i = j & k} A_{j} B_{k}
$$

$$
FWT[A] = merge(FWT[A_0] + FWT[A_1], FWT[A_1])
$$

$$
UFWT[A’] = merge(UFWT[A_0’] - UFWT[A_1’], UFWT[A_1’])
$$

$dft$ 系数 $1$，$idft$ 系数 $-1$，但是前后 $merge$ 交换了

template<int P>
void and_dft(vector<MInt<P>>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                a[j + k] += a[i + j + k]; // 与 or 相反
            }
        }
    }
}

template<int P>
void and_idft(vector<MInt<P>>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                a[j + k] -= a[i + j + k];
            }
        }
    }
}

template<int P>
vector<MInt<P>> and_fwt(vector<MInt<P>> a, vector<MInt<P>> b) {
    int n = 1, tot = max(a.size(), b.size());
    while (n < tot) {
        n *= 2;
    }
    a.resize(n);
    b.resize(n);
    and_dft(a);
    and_dft(b);
    for (int i = 0; i < n; i++) {
        a[i] *= b[i];
    }
    and_idft(a);
    return a;
}

xor

$$
C_{i} = \sum\limits_{i = j \oplus k} A_{j} B_{k}
$$

$$
FWT[A] = merge(FWT[A_0] + FWT[A_1], FWT[A_0] - FWT[A_1])
$$

逆变换易得：

$$
UFWT[A’] = merge(\frac{UFWT[A_0’] + UFWT[A_1’]}{2}, \frac{UFWT[A_0’] - UFWT[A_1’]}{2})
$$

template<int P>
void xor_dft(vector<MInt<P>>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                auto u = a[j + k];
                auto v = a[i + j + k];
                a[j + k] = u + v; // 1
                a[i + j + k] = u - v;
            }
        }
    }
}

template<int P>
void xor_idft(vector<MInt<P>>& a) {
    int n = a.size();
    auto inv = MInt<P>(2).inv();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                auto u = a[j + k];
                auto v = a[i + j + k];
                a[j + k] = (u + v) * inv; // 1/2
                a[i + j + k] = (u - v) * inv;
            }
        }
    }
}

template<int P>
vector<MInt<P>> xor_fwt(vector<MInt<P>> a, vector<MInt<P>> b) {
    int n = 1, tot = max(a.size(), b.size());
    while (n < tot) {
        n *= 2;
    }
    a.resize(n);
    b.resize(n);
    xor_dft(a);
    xor_dft(b);
    for (int i = 0; i < n; i++) {
        a[i] *= b[i];
    }
    xor_idft(a);
    return a;
}

$dft$ 系数 $1$，$idft$ 系数 $\frac{1}{2}$

xnor

$$
C_{i} = \sum\limits_{i = j \otimes k} A_{j} B_{k}
$$

$$
FWT[A] = merge(FWT[A_1] - FWT[A_0], FWT[A_1] + FWT[A_0])
$$

$$
UFWT[A’] = merge(\frac{UFWT[A_1’] - UFWT[A_0’]}{2}, \frac{UFWT[A_1’] + UFWT[A_0’]}{2})
$$

template<int P>
void xnor_dft(vector<MInt<P>>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                auto v = a[j + k]; // 与 xor 相反
                auto u = a[i + j + k];
                a[j + k] = u + v; // 1
                a[i + j + k] = u - v;
            }
        }
    }
}

template<int P>
void xnor_idft(vector<MInt<P>>& a) {
    int n = a.size();
    auto inv = MInt<P>(2).inv();
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i * 2) {
            for (int k = 0; k < i; k++) {
                auto v = a[j + k];
                auto u = a[i + j + k];
                a[j + k] = (u + v) * inv; // 1/2
                a[i + j + k] = (u - v) * inv;
            }
        }
    }
}

template<int P>
vector<MInt<P>> xnor_fwt(vector<MInt<P>> a, vector<MInt<P>> b) {
    int n = 1, tot = max(a.size(), b.size());
    while (n < tot) {
        n *= 2;
    }
    a.resize(n);
    b.resize(n);
    xnor_dft(a);
    xnor_dft(b);
    for (int i = 0; i < n; i++) {
        a[i] *= b[i];
    }
    xnor_idft(a);
    return a;
}

$dft$ 系数 $1$，$idft$ 系数 $\frac{1}{2}$

子集卷积

$$
c_{k} = \sum\ a_{i} b_{j} \quad if \quad i & j=0 \quad i \mid j=k
$$

for (int i = 0; i < (1 << n); i++) {
    int x;
    cin >> x;
    a[bitCount(i)][i] = x;
}

for (int i = 0; i < (1 << n); i++) {
    int x;
    cin >> x;
    b[bitCount(i)][i] = x;
}

for (int i = 0; i <= n; i++) {
    xor_dft(a[i]);
    xor_dft(b[i]);
}

for (int i = 0; i <= n; i++) {
    for (int j = 0; j < (1 << n); j++) {
        for (int k = 0; k <= i; k++) {
            c[i][j] += a[k][j] * b[i - k][j];
        }
    }
}

for (int i = 0; i <= n; i++) {
    xor_idft(c[i]);
}

for (int i = 0; i < (1 << n); i++) {
    cout << c[bitCount(i)][i] << " ";
}

$O(n^2 \times 2^{n})$

Elementary-function

Derivative

$$
\left(\sum_{k=0}^{+\infty}f_kx^k\right)’=\sum_{k=1}^{+\infty}kf_kx^{k-1}
$$

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

Integral

$$
\int_{k = 0}^{+\infty} f_kx^k=\sum_{k=1}^{+\infty} \frac{f_kx^{k}}{k + 1} + C
$$

$C$ 这里取了 $0$

constexpr Poly integr() const {
    Poly res(this->size() + 1);
    for (int i = 0; i < this->size(); ++i) {
        res[i + 1] = (*this)[i] / (i + 1);
    }
    return res;
}

Inverse

$$
f(x) \ast f^{-1}(x) \equiv 1 \pmod{x^{n}}
$$

求 $f^{-1}(x)$

假设有 $f_{0}^{-1}(x)$ 满足：
$$
f(x) \ast f_{0}^{-1}(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}
$$
再把原式变形一下：
$$
f(x) \ast f^{-1}(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}
$$
可得：
$$
f(x) \ast (f_{0}^{-1}(x) - f^{-1}(x)) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}
$$
由于逆元存在（$f(x) \neq 0$）：
$$
f_{0}^{-1}(x) - f^{-1}(x) \equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}
$$
两边平方可得：
$$
f^{-2}\left(x\right)-2f^{-1}\left(x\right)f^{-1}{0}\left(x\right)+f^{-2}{0}\left(x\right)\equiv 0 \pmod{x^{n}}
$$
两边同乘 $f\left(x\right)$ 并移项可得：

$$
f^{-1}\left(x\right)\equiv f^{-1}{0}\left(x\right)\left(2-f\left(x\right)f^{-1}{0}\left(x\right)\right) \pmod{x^{n}}
$$

// 获取前 k 个系数
constexpr Poly trunc(int k) const {
    Poly f = *this;
    f.resize(k);
    return f;
}

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

div

发现若能消除 $R\left(x\right)$ 的影响则可直接多项式求逆解决。

考虑构造变换

$$
f^{R}\left(x\right)=x^{\operatorname{deg}{f}}f\left(\frac{1}{x}\right)
$$

观察可知其实质为反转 $f\left(x\right)$ 的系数。

设 $n=\operatorname{deg}{f},m=\operatorname{deg}{g}$。

将 $f\left(x\right)=Q\left(x\right)g\left(x\right)+R\left(x\right)$ 中的 $x$ 替换成 $\frac{1}{x}$ 并将其两边都乘上 $x^{n}$，得到：

$$
\begin{aligned}
x^{n}f\left(\frac{1}{x}\right)&=x^{n-m}Q\left(\frac{1}{x}\right)x^{m}g\left(\frac{1}{x}\right)+x^{n-m+1}x^{m-1}R\left(\frac{1}{x}\right)\
f^{R}\left(x\right)&=Q^{R}\left(x\right)g^{R}\left(x\right)+x^{n-m+1}R^{R}\left(x\right)
\end{aligned}
$$

注意到上式中 $R^{R}\left(x\right)$ 的系数为 $x^{n-m+1}$，则将其放到模 $x^{n-m+1}$ 意义下即可消除 $R^{R}\left(x\right)$ 带来的影响。

又因 $Q^{R}\left(x\right)$ 的次数为 $\left(n-m\right)<\left(n-m+1\right)$，故 $Q^{R}\left(x\right)$ 不会受到影响。

则：

$$
f^{R}\left(x\right)\equiv Q^{R}\left(x\right)g^{R}\left(x\right)\pmod{x^{n-m+1}}
$$

使用多项式求逆即可求出 $Q\left(x\right)$，将其反代即可得到 $R\left(x\right)$。

constexpr friend Poly operator/(Poly a, Poly b) {
    if (b.empty()) {
        assert(false);
    }
    if (a.size() < b.size()) {
        return Poly();
    }
    reverse(a.begin(), a.end());
    reverse(b.begin(), b.end());
    int k = a.size() - b.size() + 1;
    auto c = a.trunc(k) * b.inv(k).trunc(k);
    reverse(c.begin(), c.end());
    return c.shift(k - (int)c.size());
}
// 保留 b.size() - 1 项，高位补 0 
constexpr friend Poly operator%(Poly a, Poly b) {
    return (a - a / b * b).trunc(b.size() - 1);
}

$\ln$

给定多项式 $f(x)$，求模 $x^{n}$ 意义下的 $\ln{f(x)}$

首先，对于多项式 $f(x)$，若 $\ln{f(x)}$ 存在，则由其定义，其必须满足：

$$
[x^{0}]f(x)=1
$$

对 $\ln{f(x)}$ 求导再积分，可得：

$$
\begin{aligned}
\frac{\mathrm{d} \ln{f(x)}}{\mathrm{d} x} & \equiv \frac{f’(x)}{f(x)} & \pmod{x^{n}} \
\ln{f(x)} & \equiv \int \mathrm{d} \ln{x} \equiv \int\frac{f’(x)}{f(x)} \mathrm{d} x & \pmod{x^{n}}
\end{aligned}
$$

1
2
3

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

$\exp$

设给定函数为 $h\left(x\right)$，有方程：

$$
g\left(f\left(x\right)\right)=\ln{f\left(x\right)}-h\left(x\right)\pmod{x^{n}}
$$

应用 Newton’s Method 可得：

$$
\begin{aligned}
f\left(x\right)&\equiv f_{0}\left(x\right)-\frac{\ln{f_{0}\left(x\right)}-h\left(x\right)}{\frac{1}{f_{0}\left(x\right)}}&\pmod{x^{n}}\
&\equiv f_{0}\left(x\right)\left(1-\ln{f_{0}\left(x\right)+h\left(x\right)}\right)&\pmod{x^{n}}
\end{aligned}
$$

constexpr 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);
}

$pow$

$$
g(x) \equiv f^{k}(x) \pmod{ x^{n}}
$$

两边取 $\ln$ ：
$$
ln(g(x)) = k \times ln(f(x)) \pmod{x^{n}}
$$
两边取 $\exp$ ：
$$
g(x) \equiv \exp{(\ln{g(x))}} \equiv \exp{(k \times \ln{f(x))}} \pmod{x^{n}}
$$

k// k <= P - 1
constexpr Poly pow(int k, int m) const {
    int i = 0;
    while (i < this->size() && (*this)[i] == 0) { // 将前面的 0 去除
        i++;
    }
    if (i == this->size() || 1LL * i * k >= m) {
        return Poly(m);
    }
    Value v = (*this)[i];
    auto f = shift(-i) * v.inv();
    return (f.log(m - i * k) * k).exp(m - i * k).shift(i * k) * power(v, k);
}

若输入 $k$ 较大

constexpr Poly pow(const string& s, int m) const {
    bool flag = false;
    long long k = 0, o = 0;
    for (auto& x : s) {
        k = k * 10 + x - '0';
        if (k >= m) {
            flag = true;
        }
        k %= P;
        o = o * 10 + x - '0';
        o %= P - 1; // 欧拉函数
    }

    if (flag && (*this)[0] == 0) {
        return Poly(m);
    }

    int i = 0;
    while (i < this->size() && (*this)[i] == 0) {
        i++;
    }
    if (i == this->size() || 1LL * i * k >= m) {
        return Poly(m);
    }
    Value v = (*this)[i];
    auto f = shift(-i) * v.inv();
    return (f.log(m - i * k) * k).exp(m - i * k).shift(i * k) * power(v, o); // 修改为 o
}

$sqrt$

给定多项式 $g\left(x\right)$，求 $f\left(x\right)$，满足：

$$
f^{2}\left(x\right)\equiv g\left(x\right) \pmod{x^{n}}
$$
首先讨论 $\left[x^0\right]g(x)$ 不为 $0$ 的情况。

易知：

$$
\left[x^0\right]f(x) = \sqrt{\left[x^0\right]g(x)}
$$

若 $\left[x^0\right]g(x)$ 没有平方根，则多项式 $g(x)$ 没有平方根。

$\left[x^0\right]g(x)$ 可能有多个平方根，选取不同的根会求出不同的 $f(x)$。

假设现在已经求出了 $g\left(x\right)$ 在模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 意义下的平方根 $f_{0}\left(x\right)$，则有：

$$
\begin{aligned}
f_{0}^{2}\left(x\right)&\equiv g\left(x\right) &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\
f_{0}^{2}\left(x\right)-g\left(x\right)&\equiv 0 &\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\
\left(f_{0}^{2}\left(x\right)-g\left(x\right)\right)^{2}&\equiv 0 &\pmod{x^{n}}\
\left(f_{0}^{2}\left(x\right)+g\left(x\right)\right)^{2}&\equiv 4f_{0}^{2}\left(x\right)g\left(x\right) &\pmod{x^{n}}\
\left(\frac{f_{0}^{2}\left(x\right)+g\left(x\right)}{2f_{0}\left(x\right)}\right)^{2}&\equiv g\left(x\right) &\pmod{x^{n}}\
\frac{f_{0}^{2}\left(x\right)+g\left(x\right)}{2f_{0}\left(x\right)}&\equiv f\left(x\right) &\pmod{x^{n}}\
2^{-1}f_{0}\left(x\right)+2^{-1}f_{0}^{-1}\left(x\right)g\left(x\right)&\equiv f\left(x\right) &\pmod{x^{n}}
\end{aligned}
$$

constexpr Poly sqrt(int m) const {
    Poly x{1};
    int k = 1;
    while (k < m) {
        k *= 2;
        x = (x + (trunc(k) * x.inv(k)).trunc(k)) * Value(2).inv();
    }
    return x.trunc(m);
}

Multipoint eval interpolation

给出一个多项式 $f\left(x\right)$ 和 $n$ 个点 $x_{1},x_{2},\dots,x_{n}$，求

$$
f\left(x_{1}\right),f\left(x_{2}\right),\dots,f\left(x_{n}\right)
$$
考虑使用分治来将问题规模减半。

将给定的点分为两部分：

$$
\begin{aligned}
X_{0}&=\left{x_{1},x_{2},\dots,x_{\left\lfloor\frac{n}{2}\right\rfloor}\right}\
X_{1}&=\left{x_{\left\lfloor\frac{n}{2}\right\rfloor+1},x_{\left\lfloor\frac{n}{2}\right\rfloor+2},\dots,x_{n}\right}
\end{aligned}
$$

构造多项式

$$
g_{0}\left(x\right)=\prod_{x_{i}\in X_{0}}\left(x-x_{i}\right)
$$

则有 $\forall x\in X_{0}:g_{0}\left(x\right)=0$。

考虑将 $f\left(x\right)$ 表示为 $g_{0}\left(x\right)Q\left(x\right)+f_{0}\left(x\right)$ 的形式，即：

$$
f_{0}\left(x\right)\equiv f\left(x\right)\pmod{g_{0}\left(x\right)}
$$

则有 $\forall x\in X_{0}:f\left(x\right)=g_{0}\left(x\right)Q\left(x\right)+f_{0}\left(x\right)=f_{0}\left(x\right)$，$X_{1}$ 同理。

至此，问题的规模被减半，可以使用分治 + 多项式取模解决。

时间复杂度

$$
T\left(n\right)=2T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log^{2}{n}\right)
$$

constexpr Poly mulT(Poly b) const {
    if (b.size() == 0) {
        return Poly();
    }
    int n = b.size();
    reverse(b.begin(), b.end());
    return ((*this) * b).shift(-(n - 1));
}
constexpr vector<Value> eval(vector<Value> x) const {
    if (this->size() == 0) {
        return vector<Value>(x.size(), 0);
    }
    const int n = max(x.size(), this->size());
    vector<Poly> q(4 * n);
    vector<Value> ans(x.size());
    x.resize(n);
    function<void(int, int, int)> build = 
        [&](int p, int l, int r) {
        if (r - l == 1) {
            q[p] = Poly{1, -x[l]};
        } else {
            int m = (l + r) / 2;
            build(2 * p, l, m);
            build(2 * p + 1, m, r);
            q[p] = q[2 * p] * q[2 * p + 1];
        }
    };
    build(1, 0, n);
    function<void(int, int, int, const Poly &)> work = 
        [&](int p, int l, int r, const Poly &num) {
        if (r - l == 1) {
            if (l < int(ans.size())) {
                ans[l] = num[0];
            }
        } else {
            int m = (l + r) / 2;
            auto lp = num.mulT(q[2 * p + 1]);
            lp.resize(m - l);
            work(2 * p, l, m, lp);
            auto rp = num.mulT(q[2 * p]);
            rp.resize(r - m);
            work(2 * p + 1, m, r, rp);
        }
    };
    work(1, 0, n, mulT(q[1].inv(n)));
    return ans;
}

Newton’s Method

给定多项式 $g\left(x\right)$，已知有 $f\left(x\right)$ 满足：

$$
g\left(f\left(x\right)\right)\equiv 0\pmod{x^{n}}
$$

求出模 $x^{n}$ 意义下的 $f\left(x\right)$。

考虑倍增。

首先当 $n=1$ 时，$\left[x^{0}\right]g\left(f\left(x\right)\right)=0$ 的解需要单独求出。

假设现在已经得到了模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 意义下的解 $f_{0}\left(x\right)$，要求模 $x^{n}$ 意义下的解 $f\left(x\right)$。

将 $g\left(f\left(x\right)\right)$ 在 $f_{0}\left(x\right)$ 处进行泰勒展开，有：

$$
\sum_{i=0}^{+\infty}\frac{g^{\left(i\right)}\left(f_{0}\left(x\right)\right)}{i!}\left(f\left(x\right)-f_{0}\left(x\right)\right)^{i}\equiv 0\pmod{x^{n}}
$$

因为 $f\left(x\right)-f_{0}\left(x\right)$ 的最低非零项次数最低为 $\left\lceil\frac{n}{2}\right\rceil$，故有：

$$
\forall 2\leqslant i:\left(f\left(x\right)-f_{0}\left(x\right)\right)^{i}\equiv 0\pmod{x^{n}}
$$

则：

$$
\sum_{i=0}^{+\infty}\frac{g^{\left(i\right)}\left(f_{0}\left(x\right)\right)}{i!}\left(f\left(x\right)-f_{0}\left(x\right)\right)^{i}\equiv g\left(f_{0}\left(x\right)\right)+g’\left(f_{0}\left(x\right)\right)\left(f\left(x\right)-f_{0}\left(x\right)\right) \equiv 0\pmod{x^{n}}
$$

$$
f\left(x\right)\equiv f_{0}\left(x\right)-\frac{g\left(f_{0}\left(x\right)\right)}{g’\left(f_{0}\left(x\right)\right)}\pmod{x^{n}}
$$

多项式求逆

设给定函数为 $h\left(x\right)$，有方程：

$$
g\left(f\left(x\right)\right)=\frac{1}{f\left(x\right)}-h\left(x\right)\equiv 0\pmod{x^{n}}
$$

应用 Newton’s Method 可得：

$$
\begin{aligned}
f\left(x\right)&\equiv f_{0}\left(x\right)-\frac{\frac{1}{f_{0}\left(x\right)}-h\left(x\right)}{-\frac{1}{f_{0}^{2}\left(x\right)}}&\pmod{x^{n}}\
&\equiv 2f_{0}\left(x\right)-f_{0}^{2}\left(x\right)h\left(x\right)&\pmod{x^{n}}
\end{aligned}
$$

时间复杂度

$$
T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)
$$

多项式开方

设给定函数为 $h\left(x\right)$，有方程：

$$
g\left(f\left(x\right)\right)=f^{2}\left(x\right)-h\left(x\right)\equiv 0\pmod{x^{n}}
$$

应用 Newton’s Method 可得：

$$
\begin{aligned}
f\left(x\right)&\equiv f_{0}\left(x\right)-\frac{f_{0}^{2}\left(x\right)-h\left(x\right)}{2f_{0}\left(x\right)}&\pmod{x^{n}}\
&\equiv\frac{f_{0}^{2}\left(x\right)+h\left(x\right)}{2f_{0}\left(x\right)}&\pmod{x^{n}}
\end{aligned}
$$

时间复杂度

$$
T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)
$$

多项式 exp

设给定函数为 $h\left(x\right)$，有方程：

$$
g\left(f\left(x\right)\right)=\ln{f\left(x\right)}-h\left(x\right)\pmod{x^{n}}
$$

应用 Newton’s Method 可得：

时间复杂度

$$
T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)
$$

多项式平移

多项式平移是简单情况的多项式复合变换，给出 $f(x)=\sum\limits_{i=0}^nf_ix^i$ 的系数和一个常数 $c$，求 $f(x+c)$ 的系数，即 $f(x)\mapsto f(x+c)$。

二项式定理法

考虑二项式定理 $\displaystyle (a+b)^n=\sum _ {i=0}^n\binom{n}{i}a^ib^{n-i}$ 那么

$$
\begin{aligned}
f(x+c)&=\sum _ {i=0}^nf_i(x+c)^i\
&=\sum _ {i=0}^nf_i\left(\sum _ {j=0}^i\binom{i}{j}x^jc^{i-j}\right)\
&=\sum _ {i=0}^nf_ii!\left(\sum _ {j=0}^i\frac{x^j}{j!}\frac{c^{i-j}}{(i-j)!}\right)\
&=\sum _ {i=0}^n\frac{x^i}{i!}\left(\sum _ {j=i}^{n}f_jj!\frac{c^{j-i}}{(j-i)!}\right)
\end{aligned}
$$

$$
\begin{align}
\sum\limits_{j=i}^{n}f_jj!\frac{c^{j-i}}{(j-i)!} &= \sum\limits_{k = 0}^{n - i} f_{k + i}(k + i)! \frac{c^{k}}{k!} \
&= \sum\limits_{k = 0}^{n - i} f^{rev}_{n - i - k}(k + i)! \frac{c^{k}}{k!}
\end{align}
$$

构造：
$$
A(x) = \sum f_{i} i! x^{i}\
B(x) = \sum \frac{c^{i}}{i!} x^{n - i}
$$
$A(x)B(x)$ 最后除 $i!$ 即为平移后的多项式。

$O(n\log(n))$

// f(x) -> f(x + c)
template<int P>
Poly<P> shift(const Poly<P>& a, int c) {
    int n = a.size() - 1; // 最高次数
    vector<MInt<P>> fac(n + 1, 1), invfac(n + 1, 1), p(n + 1, 1);
    for (int i = 1; i <= n; i++) {
        fac[i] = fac[i - 1] * i;
        p[i] = p[i - 1] * c;
    }
    invfac[n] = fac[n].inv();
    for (int i = n; i > 1; i--) {
        invfac[i - 1] = invfac[i] * i;
    }
    Poly<P> A(n + 1), B(n + 1);
    for (int i = 0; i <= n; i++) {
        A[i] = a[i] * fac[i];
        B[i] = p[n - i] * invfac[n - i];
    }
    A *= B;
    for (int i = n; i <= n * 2; i++) {
        A[i] *= invfac[i - n];
    }
    return A.shift(-n);
}

连续点值平移

给定一个不超过 $n$ 次的多项式的 $n+1$ 个点值 $f(0),f(1) \dots f(n)$，和一个正整数 $m$，求 $f(m),f(m+1) \dots f(m+n)$。

考虑拉格朗日插值：
$$
\begin{align}
F(x) &= \sum\limits_{0 \le i \le n} f(i) \prod\limits_{0 \le j \le n \wedge j \neq i} \frac{x - x_{j}}{x_{i} - x_{j}} \
&= \sum\limits_{0 \le i \le n} f(i) \prod\limits_{0 \le j \le n \wedge j \neq i} \frac{x - j}{i - j} \
&= \sum\limits_{0 \le i \le n} \frac{f(i)}{i!} \frac{(-1)^{n - i}}{(n - i)!} \frac{x!}{(x - n - 1)!(x - i)} \
&= \frac{x!}{(x - n - 1)!} \sum\limits_{0 \le i \le n} \frac{f(i) (-1)^{n - i}}{i! (n - i)!} \frac{1}{x - i}
\end{align}
$$
令
$$
A(x) = \sum\limits_{0 \le i \le n} \frac{f(i) (-1)^{n - i}}{i! (n - i)!} x^{i} \
B(x) = \sum\limits_{i \ge 0} \frac{1}{m - n + i} x^{i}
$$

$$
\begin{align}
x^{n + m} &= \sum\limits_{0 \le i \le n} \frac{f(i) (-1)^{n - i}}{i! (n - i)!} \frac{1}{m + x - i} \
&= \frac{(m + x - n - 1)!}{(m + x)!} f(m + x)
\end{align}
$$

int n, m;
cin >> n >> m;
Comb comb(n + 1);
Poly<P> a(n + 1);
auto finv = MInt<P>(-1).inv();
for (int i = 0; i <= n; i++) {
    cin >> a[i];
    a[i] *= comb.invfac(i) * comb.invfac(n - i) * ((n - i) % 2 == 0 ? 1 : finv);
}
Poly<P> b(n * 2 + 1); // n * 2 + 1 项
for (int i = 0; i <= n * 2; i++) {
    b[i] = MInt<P>(m - n + i).inv();
}
a *= b;
vector<MInt<P>> upfac(n * 2 + 2, 1);
for (int i = 1; i <= n * 2 + 1; i++) {
    upfac[i] = upfac[i - 1] * (m - n - 1 + i);
}
for (int i = n; i <= 2 * n; i++) {
    cout << a[i] * upfac[i + 1] / upfac[i - n] << " ";
}

Chirp Z Transform

给出多项式 $A(x)=\sum\limits_{i=0}^na_ix^i\in\mathbb{C}\lbrack x\rbrack$ ，求出 $A(1),A(c),A(c^2),\dots$ 的一种算法。

对于非负整数 $k$ 和 $i$ 考虑

$$
ki=\binom{i+k}{2}-\binom{i}{2}-\binom{k}{2}
$$

其中 $\binom{a}{b}=\frac{a(a-1)\cdots (a-b+1)}{b!}$ 为二项式系数，那么

$$
A(c^k)=c^{-\binom{k}{2}}\sum _ {i=0}^na_ic^{\binom{i+k}{2}-\binom{i}{2}}
$$

令 $A_0(x)=\sum\limits_{i}a_{n-i}c^{-\binom{n-i}{2}}x^i$ 且对于 $\forall j\gt n$ 和 $\forall j\lt 0$ 令 $a_j=0$、$B_0(x)=\sum\limits_{i\geq 0}c^{\binom{i}{2}}x^i$ 那么对于 $t\geq 0$ 有

$$
\begin{aligned}
\lbrack x^{n+t}\rbrack (A_0(x)B_0(x))&=\sum _ {i=0}^{n+t}\left(\lbrack x^{n+t-i}\rbrack A_0(x)\right)\left(\lbrack x^{i}\rbrack B_0(x)\right)\
&=\sum_{i=0}^{n+t}a_{i-t}c^{\binom{i}{2}-\binom{i-t}{2}}\
&=\sum_{i=-t}^na_ic^{\binom{i+t}{2}-\binom{i}{2}}\
&=c^{\binom{t}{2}}\cdot A(c^t)
\end{aligned}
$$

通过计算 $c^{-\binom{t}{2}}\lbrack x^{n+t}\rbrack (A_0(x)B_0(x))$ 可得到 $A(1),A(c),\dots$，该算法需一次卷积。且 $\forall i\geq 0$ 有 $c^{\binom{i+1}{2}}=c^{\binom{i}{2}}\cdot c^i$，可递推计算。

简单来说构造 $A_0(x), B_0(x)$ ：
$$
A_0(x)=\sum\limits_{i}a_{n-i}c^{-\binom{n-i}{2}}x^i \
B_0(x)=\sum\limits_{i\geq 0}c^{\binom{i}{2}}x^i \
A_0(x)B_0(x) = c^{\binom{t}{2}}\cdot A(c^t)
$$
$t$ 为 $B_0(x)$ 的最高项。

给定一个 $n + 1$ 项多项式 $P(x)$ 以及 $c, m$，请计算 $P(c^0),P(c^1),\dots,P(c^{m-1})$。

long long n, c, m;
cin >> n >> c >> m;
Poly<P> a(n + 1);
for (int i = 0; i <= n; i++) {
    cin >> a[i];
}

vector<MInt<P>> p(n + m + 1, 1); // c^i;
vector<MInt<P>> cp(n + m + 1, 1); // c^(binom{i}{2})
for (int i = 1; i <= n + m; i++) {
    p[i] = p[i - 1] * c;
    cp[i + 1] = cp[i] * p[i];
}

Poly<P> A = a;
reverse(A.begin(), A.end());
for (int i = 0; i <= n; i++) {
    A[i] /= cp[n - i];
}

Poly<P> B(n + m + 1);
for (int i = 0; i <= n + m; i++) {
    B[i] = cp[i];
}
A *= B;

for (int i = n; i < n + m; i++) {
    cout << A[i] / cp[i - n] << " ";
}

Poly Template

template<class T>
constexpr T power(T a, long long b) {
    T res = 1;
    for (; b; a *= a, b >>= 1) {
        if (b & 1) {
            res *= a;
        }
    }
    return res;
}

template<int P>
struct MInt {
	unsigned x;
	constexpr MInt() : x(0) {}
    constexpr MInt(long long sig) : x(norm(sig % P)) {}
	constexpr int val() const { 
        return x; 
    }
    // -P <= v < 2P
    constexpr long long norm(long long v) const {
        if (v < 0) v += P;
        if (v >= P) v -= P;
        return v;
    }
    constexpr MInt inv() const {
        return power(*this, P - 2);
    }
	constexpr MInt& operator+=(MInt that) { 
        x = norm(x + that.x);
        return *this; 
    }
	constexpr MInt& operator-=(MInt that) {
        x = norm(x + P - that.x);
        return *this; 
    }
	constexpr MInt& operator*=(MInt that) {
        x = (unsigned long long)x * that.x % P; 
        return *this; 
    }
	constexpr MInt& operator/=(MInt that) { 
        return (*this) *= that.inv(); 
    }
    constexpr friend MInt operator+(MInt a, MInt b) { return a += b; }
    constexpr friend MInt operator-(MInt a, MInt b) { return a -= b; }
    constexpr friend MInt operator*(MInt a, MInt b) { return a *= b; }
    constexpr friend MInt operator/(MInt a, MInt b) { return a /= b; }
    constexpr friend MInt operator-(MInt a) { return MInt() -= a; } // 负号
    constexpr friend bool operator<(MInt a, MInt b) { return a.x < b.x; }
    constexpr friend bool operator<=(MInt a, MInt b) { return a.x <= b.x; }
    constexpr friend bool operator==(MInt a, MInt b) { return a.x == b.x; }
    constexpr friend bool operator!=(MInt a, MInt b) { return a.x != b.x; }
	friend ostream& operator<<(ostream& os, MInt a) { 
        os << a.x; 
        return os; 
    }
    friend istream& operator>>(istream& is, MInt& a) { 
        long long x; 
        is >> x;
        a = MInt(x); 
        return is; 
    }
};

// P 质数
template<int P>
constexpr MInt<P> findPrimitiveRoot() {
    long long x = P - 1;
    vector<int> factors;
    for (long long i = 2; i * i <= x; i++) {
        if (x % i == 0) {
            factors.push_back(i);
            while (x % i == 0) {
                x /= i;
            }
        }
    }
    if (x > 1) {
        factors.push_back(x);
    }
    for (int i = 1; i < P; i++) {
        if (power(MInt<P>(i), P - 1) == 1) {
            bool flag = true;
            for (auto& f : factors) {
                if (power(MInt<P>(i), (P - 1) / f) == 1) {
                    flag = false;
                    break;
                }
            }
            if (flag) {
                return i;
            }
        }
    }
    assert(false);
}

template<int P>
constexpr MInt<P> primitiveRoot = findPrimitiveRoot<P>();

template<int P>
void dft(vector<MInt<P>>& a) {
    int n = a.size();
    vector<int> rev(n);
    for (int i = 0; i < n; i++) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) {
            rev[i] |= n >> 1;
        }
        if (i > rev[i]) {
            swap(a[i], a[rev[i]]);
        }
    }
    for (int i = 1; i < n; i <<= 1) {
        MInt<P> wn = power(primitiveRoot<P>, (P - 1) / (i * 2));
        for (int j = 0; j < n; j += i * 2) {
            MInt<P> w = 1;
            for (int k = 0; k < i; k++) {
                MInt<P> u = a[k + j];
                MInt<P> v = a[k + j + i] * w;
                a[k + j] = u + v;
                a[k + j + i] = u - v;
                w *= wn;
            }
        }
    }
}

template<int P>
void idft(vector<MInt<P>>& a) {
    int n = a.size();
    reverse(a.begin() + 1, a.end());
    dft(a);
    MInt<P> inv = MInt<P>(n).inv();
    for (int i = 0; i < n; i++) {
        a[i] *= inv;
    }
}

template<int P>
struct Poly : public vector<MInt<P>> {
    using Value = MInt<P>;
    
    Poly() : vector<Value>() {}
    explicit constexpr Poly(int n) : vector<Value>(n) {}
    
    explicit constexpr Poly(const vector<Value> &a) : vector<Value>(a) {}
    constexpr Poly(const initializer_list<Value> &a) : vector<Value>(a) {}
    
    template<class InputIt, class = _RequireInputIter<InputIt>>
    explicit constexpr Poly(InputIt first, InputIt last) : vector<Value>(first, last) {}
    
    template<class F>
    explicit constexpr Poly(int n, F f) : vector<Value>(n) {
        for (int i = 0; i < n; i++) {
            (*this)[i] = f(i);
        }
    }
    // 系数移动 (k > 0 右移)
    constexpr 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());
        }
    }
    // 获取前 k 个系数
    constexpr Poly trunc(int k) const {
        Poly f = *this;
        f.resize(k);
        return f;
    }
    constexpr friend Poly operator+(const Poly &a, const Poly &b) {
        Poly res(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] += b[i];
        }
        return res;
    }
    constexpr friend Poly operator-(const Poly &a, const Poly &b) {
        Poly res(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] -= b[i];
        }
        return res;
    }
    constexpr friend Poly operator-(const Poly &a) {
        std::vector<Value> res(a.size());
        for (int i = 0; i < res.size(); i++) {
            res[i] = -a[i];
        }
        return Poly(res);
    }
    constexpr friend Poly operator*(Poly a, Poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return Poly();
        }
        int n = 1, tot = a.size() + b.size() - 1;
        while (n < tot) {
            n *= 2;
        }
        a.resize(n);
        b.resize(n);
        dft(a);
        dft(b);
        for (int i = 0; i < n; ++i) {
            a[i] *= b[i];
        }
        idft(a);
        a.resize(tot);
        return a;
    }
    constexpr friend Poly operator*(Value a, Poly b) {
        for (int i = 0; i < b.size(); i++) {
            b[i] *= a;
        }
        return b;
    }
    constexpr friend Poly operator*(Poly a, Value b) {
        for (int i = 0; i < a.size(); i++) {
            a[i] *= b;
        }
        return a;
    }
    constexpr friend Poly operator/(Poly a, Poly b) {
        if (b.empty()) {
            assert(false);
        }
        if (a.size() < b.size()) {
            return Poly();
        }
        reverse(a.begin(), a.end());
        reverse(b.begin(), b.end());
        int k = a.size() - b.size() + 1;
        auto c = a.trunc(k) * b.inv(k).trunc(k);
        reverse(c.begin(), c.end());
        return c.shift(k - (int)c.size());
    }
    constexpr friend Poly operator/(Poly a, Value b) {
        for (int i = 0; i < a.size(); i++) {
            a[i] /= b;
        }
        return a;
    }
    // 保留 b.size() - 1 项，高位补 0 
    constexpr friend Poly operator%(Poly a, Poly b) {
        return (a - a / b * b).trunc(b.size() - 1);
    }
    constexpr Poly &operator+=(Poly b) {
        return (*this) = (*this) + b;
    }
    constexpr Poly &operator-=(Poly b) {
        return (*this) = (*this) - b;
    }
    constexpr Poly &operator*=(Poly b) {
        return (*this) = (*this) * b;
    }
    constexpr Poly &operator*=(Value b) {
        return (*this) = (*this) * b;
    }
    constexpr Poly &operator/=(Poly b) {
        return (*this) = (*this) / b;
    }
    constexpr Poly &operator/=(Value b) {
        return (*this) = (*this) / b;
    }
    constexpr Poly &operator%=(Poly b) {
        return (*this) = (*this) % b;
    }
    constexpr Poly deriv() const {
        if (this->empty()) {
            return Poly();
        }
        Poly res(this->size() - 1);
        for (int i = 0; i < this->size() - 1; ++i) {
            res[i] = (i + 1) * (*this)[i + 1];
        }
        return res;
    }
    constexpr Poly integr() const {
        Poly res(this->size() + 1);
        for (int i = 0; i < this->size(); ++i) {
            res[i + 1] = (*this)[i] / (i + 1);
        }
        return res;
    }
    constexpr Poly inv(int m) const {
        Poly x{(*this)[0].inv()};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{2} - trunc(k) * x)).trunc(k);
        }
        return x.trunc(m);
    }
    constexpr Poly log(int m) const {
        return (deriv() * inv(m)).integr().trunc(m);
    }
    constexpr 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);
    }
    // k <= P - 1 的情况
    constexpr 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);
        }
        Value v = (*this)[i];
        auto f = shift(-i) * v.inv();
        return (f.log(m - i * k) * k).exp(m - i * k).shift(i * k) * power(v, k);
    }
    constexpr Poly sqrt(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (trunc(k) * x.inv(k)).trunc(k)) * Value(2).inv();
        }
        return x.trunc(m);
    }
    constexpr Poly mulT(Poly b) const {
        if (b.size() == 0) {
            return Poly();
        }
        int n = b.size();
        reverse(b.begin(), b.end());
        return ((*this) * b).shift(-(n - 1));
    }
    constexpr vector<Value> eval(vector<Value> x) const {
        if (this->size() == 0) {
            return vector<Value>(x.size(), 0);
        }
        const int n = max(x.size(), this->size());
        vector<Poly> q(4 * n);
        vector<Value> ans(x.size());
        x.resize(n);
        function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = Poly{1, -x[l]};
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        function<void(int, int, int, const Poly &)> work = [&](int p, int l, int r, const Poly &num) {
            if (r - l == 1) {
                if (l < int(ans.size())) {
                    ans[l] = num[0];
                }
            } else {
                int m = (l + r) / 2;
                auto lp = num.mulT(q[2 * p + 1]);
                lp.resize(m - l);
                work(2 * p, l, m, lp);
                auto rp = num.mulT(q[2 * p]);
                rp.resize(r - m);
                work(2 * p + 1, m, r, rp);
            }
        };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
};

LinearRecurrence

给定一个线性递推数列 ${f_i}$ 的前 $k$ 项 $f_0\dots f_{k-1}$，和其递推式 $f_n=\sum_{i=1}^k f_{n-i}g_i$ 的各项系数 $g_i$，求 $f_n$。

定义 $F(\sum c_ix^i)=\sum c_if_i$，那么答案就是 $F(x^n)$。

由于 $f_n=\sum_{i=1}^{k}f_{n-i}a_i$，所以 $F(x^n)=F(\sum_{i=1}^{k}a_ix^{n-i})$，所以 $F(x^n-\sum_{i=1}^k a_ix^{n-i})=F(x^{n-k}(x^k-\sum_{i=0}^{k-1}a_{k-i}x^i))=0$。

设 $G(x)=x^k-\sum_{i=0}^{k-1}a_{k-i}x^i$。

那么 $F(A(x)+x^nG(x))=F(A(x))+F(x^nG(x))=F(A(x))$。

那么就可以通过多次对 $A(x)$ 加上 $x^nG(x)$ 的倍数来降低 $A(x)$ 的次数。

也就是求 $F(A(x)\bmod G(x))$。$A(x)\bmod G(x)$ 的次数不超过 $k-1$，而 $f_{0..k-1}$ 已经给出了，就可以算了。

问题转化成了快速地求 $x^n\bmod G(x)$，只要将普通快速幂中的乘法与取模换成多项式乘法与多项式取模就可以在 $O(k\log k\log n)$ 的时间复杂度内解决这个问题了。

简单来说：求 $x^{n} \pmod{G(x)}$ ，$G(x) = x^{k} - \sum_{i=0}^{k-1}a_{k-i}x^i$ ，最后累加答案 $f_{n} \equiv \sum\limits_{i} f_{i}b_{i} \pmod{P}$

template<int P>
MInt<P> linearRecurrence(vector<MInt<P>> f, vector<MInt<P>> g, long long n) {
    Poly<P> a{0, 1};
    Poly<P> b{1};
    reverse(g.begin(), g.end());
    for (auto& x : g) {
        x *= -1;
    }
    g.push_back(1);
    const Poly<P> p(g.begin(), g.end());
    while (n) {
        if (n & 1) {
            b = b * a % p;
        }
        a = a * a % p;
        n >>= 1;
    }
    MInt<P> res = 0;
    for (int i = 0; i < f.size(); i++) {
        res += f[i] * b[i];
    }
    return res;
}

相较于矩阵快速幂 $O(k^{3} \log(n))$ 更快

generating function

Ordinary generating function (OGF)

$$
F(x) = \sum_i a_{i}x^{i}
$$

系数可用于求多重组合数

卷积
$$
F(x)G(x) = \sum_n x^n \sum_{i=0}^na_ib_{n-i}
$$
封闭形式

系数 $<1,1,\cdots,1>$
$$
F(x) = 1 + x + x^{2} + \cdots + x^{n} = \frac{1}{1 - x}
$$

由 $F(x) x + 1 = F(x)$ 推出
系数 $<1,p,\cdots,p^{n}>$
$$
F(x) = 1 + px + p^{2}x^{2} + \cdots + p^{n}x^{n} = \frac{1}{1 - px}
$$
系数 $<1,2,\cdots,n>$
$$
F(x) = 1 + 2x + 3x^{2} + \cdots + x^{n - 1} = \frac{1}{(1 - x)^{2}}
$$
对第一个式子求导即可
系数 $\binom{m}{n}$
$$
F(x)=\sum_{n\ge 0}\binom{m}{n}x^n=(1 + x)^{m}
$$
系数 $\binom{n + m - 1}{n}$
$$
F(x) = \sum_{n \ge 0} \binom{n + m - 1}{n} x^{n} = \frac{1}{(1 - x)^m}
$$
该式为广义二项式定理

Exponential generating function (EGF)

$$
F(x) = \sum_{n} a_{n} \frac{x^{n}}{n!}
$$

系数可用于求多重排列数

卷积
$$
\begin{align}
F(x)G(x) &= \sum_{i \ge 0} a_{i} \frac{x^{i}}{i!} \sum_{j \ge 0} b_{j} \frac{x^{j}}{j!} \
&= \sum_{n} x^{n} \sum_{i = 0}^{n} a_{i}b_{n - i} \frac{1}{i!(n - i)!} \
&= \sum_{n} \frac{x^{n}}{n!} \sum_{i = 0}^{n} \binom{n}{i} a_{i}b_{n - i}
\end{align}
$$
封闭形式

系数 $<1, 1, \cdots, 1>$
$$
F(x) = 1 + \frac{x^{1}}{1!} + \frac{x^{2}}{2!} + \cdots + \frac{x^{n}}{n!} = \mathrm{e}^{x}
$$
系数 $<1, p, \cdots, p^{n}>$
$$
F(x) = 1 + p\frac{x^{1}}{1!} + p^{2}\frac{x^{2}}{2!} + \cdots + p^{n}\frac{x^{n}}{n!} = \mathrm{e}^{px}
$$
系数 $<1, -1, \cdots, (-1)^{n}>$
$$
F(x) = 1 - \frac{x^{1}}{1!} + \frac{x^{2}}{2!} + \cdots + (-1)^{n} \frac{x^{n}}{n!} = \mathrm{e}^{-x}
$$
系数 $<1, 0, \cdots, (n + 1) % 2>$
$$
F(x) = 1 + \frac{x^{2}}{2!} + \frac{x^{4}}{4!} + \cdots + \frac{x^{2n}}{(2n)!} = \frac{\mathrm{e}^{x} + \mathrm{e}^{-x}}{2}
$$
系数 $<0, 1, \cdots, n % 2>$
$$
F(x) = 1 + \frac{x^{1}}{1!} + \frac{x^{3}}{3!} + \cdots + \frac{x^{2n - 1}}{(2n - 1)!} = \frac{\mathrm{e}^{x} - \mathrm{e}^{-x}}{2}
$$

Dirichlet series generating function (DGF)

$$
F(x) = \sum_{n \ge 1} \frac{a_{n}}{n^{x}}
$$

卷积
$$
\begin{align}
F(x)G(x) &= \sum_{i \ge 1} \frac{a_{i}}{i^{x}} \sum_{j \ge 1} \frac{b_{j}}{j^{x}} \
&= \sum_{n \ge 1} \frac{1}{n^{x}} \sum_{d \mid n} a_{d} b_{\frac{n}{d}}
\end{align}
$$

$$
f * g = F(d)G(\frac{n}{d}) = F(\frac{n}{d})G(d)
$$

交换律： $fg=gf$。

结合律：$(fg)h=f(gh)$。

分配律：$(f+g)h=fh+g*h$。

积性函数

若函数 $f(n)$ 满足 $f(1)=1$ 且 $\forall x,y\in\mathbf{N}^*,~(x,y)=1$ 都有 $f(xy)=f(x)f(y)$，则 $f(n)$ 为 积性函数。
若函数 $f(n)$ 满足 $f(1)=1$ 且 $\forall x,y\in\mathbf{N}^*$ 都有 $f(xy)=f(x)f(y)$，则 $f(n)$ 为 完全积性函数。

欧拉函数和莫比乌斯函数都是积性函数
$$
\sum\limits_{d \mid n} \varphi(d) = n
$$

$$
\sum\limits_{d \mid n} \mu(d) = [n = 1]
$$

欧拉函数与莫比乌斯函数的联系
$$
\sum\limits_{d \mid n} \mu(d) \frac{n}{d} = \varphi(n)
$$
常用函数
$$
\begin{align}
& \text{元函数：} \quad \varepsilon(n) = [n = 1] \
& \text{常函数：} \quad 1(n) = 1 \
& \text{恒等函数：} \quad id(n) = n
\end{align}
$$
卷积关系
$$
\begin{align}
& \mu \ * \ 1 = \varepsilon \
& \varphi \ * \ 1 = id \
& \mu \ * \ id = \varphi \
& f \ * \ \varepsilon = f \
& f \ * \ 1 \neq f \
& 1 \ * \ id = \sigma \text{约数和}
\end{align}
$$