多项式二环魔法咏唱

多项式三角函数

考虑欧拉定理

\[e^{ix} = \cos{x} + i\sin{x} \]

则

\[e^{-ix} = \cos{x} - i\sin{x} \]

因此

\[\cos{x} = \frac{e^{ix} + e^{-ix}}{2} \]

\[\sin{x} = \frac{e^{ix} - e^{-ix}}{2i} \]

其中 \(i \equiv 86583718 \pmod {998244353}\)，因此套用 \(\exp\) 板子即可。

void Poly_sin(int *f, int n) {
    static int sav[N << 1];
    for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * _i % mod;
    for (int i = 0; i < n; i ++ ) sav[i] = (mod - f[i]) % mod;
    Poly_exp(f, n); Poly_exp(sav, n);
    for (int i = 0; i < n; i ++ ) f[i] = 1ll * (f[i] - sav[i] + mod) * invi % mod * inv2 % mod;
    clr(sav, n);
}

void Poly_cos(int *f, int n) {
    static int sav[N << 1];
    for (int i = 0; i < n; i ++ ) f[i] = 1ll * f[i] * _i % mod;
    for (int i = 0; i < n; i ++ ) sav[i] = (mod - f[i]) % mod;
    Poly_exp(f, n); Poly_exp(sav, n);
    for (int i = 0; i < n; i ++ ) f[i] = 1ll * (f[i] + sav[i]) * inv2 % mod;
    clr(sav, n);
}

任意模数多项式乘法

MTT

三模数 NTT 固然可以，我们来试试拆系数 FFT（MTT），考虑 \(F(x), \, G(x)\) 两个多项式改写为 \(F(x) = A(x) + c \times B(x), \, G(x) = C(x) + c \times D(x)\)，其中 \(c = 32768\)。

为计算 \(F(x)G(x)\)，我们有

\[F(x)G(x) = A(x)C(x) + c \times (A(x)D(x) + B(x)C(x)) + \times c^2 \times B(x)D(x) \]

令 \(T(x) = C(x) + i \times D(x)\)，我们可以通过 \(T(x)\) 与 \(A(x), \, B(x)\) 的卷积得到以上所有系数。上述算法使用了 \(3\) 次 DFT，\(2\) 次 IDFT，相较 \(9\) 次 NTT 的三模数 NTT 来说慢那么一点，事实上，可以优化到 \(4\) 次 FFT，这样可以比 NTT 更优。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ld long double
const int N = 2e5 + 10;
const ld PI = acos(-1);

struct Complex {
    ld x, y;
    Complex(ld xx = 0, ld yy = 0) {x = xx, y = yy;}
    Complex operator+ (const Complex &C) const {return Complex(x + C.x, y + C.y);}
    Complex operator- (const Complex &C) const {return Complex(x - C.x, y - C.y);}
    Complex operator* (const Complex &C) const {return Complex(C.x * x - C.y * y, C.y * x + C.x * y);}
    Complex operator* (const ld n) const {return Complex(x * n, y * n);}
    Complex operator/ (const ld n) const {return Complex(x / n, y / n);}
    Complex operator~ () const {return Complex(x, -y);}
}F[N << 1], G[N << 1], T[N << 1];

int rev_len, rev[N << 1];

void rev_init(int n) {
    if (n == rev_len) return;
    for (int i = 0; i < n; i ++ ) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
    rev_len = n;
}

void FFT(Complex *f, int op, int n) {
    rev_init(n);
    static Complex wk[N << 1] = {Complex(1, 0)};
    for (int i = 0; i < n; i ++ ) if (i < rev[i]) swap(f[i], f[rev[i]]);
    for (int k = 1; k < n; k <<= 1) {
        Complex w1 = Complex(cos(PI / k), op * sin(PI / k));
        for (int i = 1; i < k; i ++ ) wk[i] = wk[i - 1] * w1;
        for (int i = 0; i < n; i += k << 1) {
            for (int j = 0; j < k; j ++ ) {
                Complex x = f[i | j], y = wk[j] * f[i | j | k];
                f[i | j] = x + y, f[i | j | k] = x - y;
            }
        }
    }
    if (op == -1) for (int i = 0; i < n; i ++ ) f[i] = f[i] / n;
}

int n, m, mod;
int a[N << 1];

void solve() {
    cin >> n >> m >> mod;
    for (int i = 0; i <= n; i ++ ) cin >> a[i], F[i].x = a[i] & ((1 << 15) - 1), G[i].x = a[i] >> 15;
    for (int i = 0; i <= m; i ++ ) cin >> a[i], T[i].x = a[i] & ((1 << 15) - 1), T[i].y = a[i] >> 15;
    int len; for (len = 1; len <= n + m; len <<= 1);
    FFT(F, 1, len); FFT(G, 1, len); FFT(T, 1, len);
    for (int i = 0; i < len; i ++ ) F[i] = F[i] * T[i], G[i] = G[i] * T[i];
    FFT(F, -1, len); FFT(G, -1, len);
    for (int i = 0; i <= n + m; i ++ ) a[i] = ((ll)round(F[i].x) % mod + (((ll)(round(F[i].y) + round(G[i].x)) % mod) << 15) % mod + (((ll)round(G[i].y) % mod) << 30) % mod) % mod;
    for (int i = 0; i <= n + m; i ++ ) cout << (a[i] % mod + mod) % mod << " ";
}

分治 FFT

\[f_i=\sum_{j=1}^if_{i-j}g_j \]

类似于 CDQ 分治，由于后面的项由前面确定，每次求一个区间的时候，要保证这个区间左边所有值对这个区间的贡献都算出来了，然后用 \(f\) 在 \([l, \, mid]\) 的段乘上 \(g[0, \, r - l]\) 实现左半段到右半段的贡献计算，分治递归即可。

void CDQ(int *f, int *g, int l, int r) {
    static int b1[N << 1], b2[N << 1];
    if (l == r) {
        if (!l) f[l] = 1;
        return;
    }
    int mid = l + r >> 1;
    CDQ(f, g, l, mid);
    int n; for (n = 1; n <= r - l + 1; n <<= 1);
    cpy(b1, f + l, mid - l + 1); clr(b1 + mid - l + 1, n - (mid - l));
    cpy(b2, g, r - l + 1); clr(b2 + r - l + 1, n - (r - l));
    NTT(b1, 1, n); NTT(b2, 1, n); px(b1, b2, n); NTT(b1, -1, n);
    for (int i = mid + 1; i <= r; i ++ ) f[i] = (f[i] + b1[i - l]) % mod;
    clr(b1, n); clr(b2, n);
    CDQ(f, g, mid + 1, r);
}