MTT | 拆系数FFT

当模数抽抽又象象的时候，使用拆系数 FFT。

注意到如果是 \(10^9\) 的话不取模可以达到高贵的 \(10^{23}\)，有人说那不能用 i128 吗？不行，因为要取一个 \(10^{23}\) 向上的 ntt 模数，然鹅这个东西中间计算的过程会打到高贵的 \(10^{46}\)，i128 自愧不如。别问我怎么知道的。

于是乎我们尝试使用 FFT 解决。

但 FFT 精度炸炸炸！我们考虑把这个多项式拆一拆，然后最后再乘一乘。

具体来讲取 \(a_{0,i}=\lfloor \cfrac{a_i}{T} \rfloor,a_{1,i}=a_i \mod T\)，至于 \(b\) 也是一样的！。

于是乎我们发现 \(C(x)=T^2A_0(x)B_0(x) + TA_1(x)B_0(x) + TA_0(x)B_1(x) + A_1(x)B_1(x)\)。

然后就可以直接算啦，结果傲娇的评测姬一只脚踩在你脸上扔给了你一个 TLE。

我们求的是四个多项式乘法对不对？这个 FFT 的次数达到了恐怖的 8 次。

我们考虑偷偷懒。发现在 FFT 中我们传进去的数的虚部是 \(0\)，我们把这个虚部利用起来。

假设现在要算 \(A,B\) 的 FFT，我们可以构造 \(P=A-Bi,Q=A+Bi\)。然后对 \(P\) FFT，但 \(Q\) 共轭于 \(P\) 所以就也可以得出 \(Q\) 的点值表达，最后解一个幼儿园方程就可以啦！

好现在我们可以用两次 FFT 算出四个多项式的点值表达了。

然后我们还原回去我们可以类似考虑，我们发现一个假设以 \(i\) 为未知数的方程是可以通过两个方程解出来的！然后我们就构造，把一个乘上 \(i\) 再加上另一个，最后再解个方程回来就好啦！

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

#define int long long
#define double long double
#define endl '\n'
#define debug(x) cerr << #x << " = " << x << endl
#define rep(i, a, b) for (int i = (a); i <= (b); i++)
#define per(i, a, b) for (int i = (a); i >= (b); i--)
#define gn(u, v) for (int v : G.G[u])
#define pb emplace_back
#define mp make_pair
#define fi first
#define se second
#define sz(x) (int)(x).size()
#define pii pair<int, int>
#define vi vector<int>
#define vpii vector<pii>
#define vvi vector<vi>
#define no cout << "NO" << endl
#define yes cout << "YES" << endl
#define all(x) x.begin(), x.end()
#define rall(x) x.rbegin(), x.rend()
#define tomin(x, y) ((x) = min((x), (y)))
#define tomax(x, y) ((x) = max((x), (y)))
#define ck(mask, i) (((mask) >> (i)) & 1)
#define pq priority_queue
#define umap unordered_map
#define FLG (cerr << "Alive!" << endl);
#define popcount __builtin_popcount

int n, m, mod, len;

constexpr int MAXN = 1e6 + 5;
constexpr double EPS = 1e-8;
const double PI = acos(-1.0);
struct Complex {
    double r, i;
    Complex() { r = i = 0; }
    Complex(double _r) {
        r = _r, i = 0;
    }
    Complex(double _r, double _i) {
        r = _r, i = _i;
    }
    friend Complex operator + (Complex x, Complex y) {
        return Complex{ x.r + y.r, x.i + y.i };
    }
    friend Complex operator - (Complex x, Complex y) {
        return Complex{ x.r - y.r, x.i - y.i };
    }
    friend Complex operator * (Complex x, Complex y) {
        return Complex{ x.r * y.r - x.i * y.i, x.r * y.i + x.i * y.r };
    }
    Complex conj() {
        return {r, -i};
    }
};
Complex I = {0, 1};

int rev[MAXN];

void FFT(Complex a[], int n, int typ) {
    // if (typ == -1)
    //     rep (i, 1, n - 1)
    //         if (i < n - i)
    //             swap(a[i], a[n - i]);
    rep (i, 0, n - 1) if (i < rev[i])
        swap(a[i], a[rev[i]]);
    for (int i = 1; i < n; i <<= 1) {
        Complex omega = Complex{ cos(PI / i), typ * sin(PI / i) };
        for (int j = 0; j < n; j += i << 1) {
            Complex w = { 1, 0 };
            for (int k = 0; k < i; k++, w = w * omega) {
                Complex p = a[j + k], q = w * a[j + i + k];
                a[j + k] = p + q, a[j + i + k] = p - q;
            }
        }
    }
    if (typ == -1)
        rep (i, 0, n - 1)
            a[i] = {a[i].r / n, a[i].i / n};
}

void doubleFFT(Complex a[], Complex b[], int len, int t) {
    rep (i, 0, len - 1)
        a[i] = a[i] + I * b[i];
    FFT(a, len, t);
    // FLG
    rep (i, 0, len - 1) b[i] = a[(len - i) % len].conj();
    rep (i, 0, len - 1) {
        Complex p = a[i], q = b[i];
        a[i] = (p + q) * 0.5;
        b[i] = (q - p) * 0.5 * I;
    }
}

Complex a0[MAXN], b0[MAXN], a1[MAXN], b1[MAXN], p[MAXN], q[MAXN];
int ans[MAXN];

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);

    cin >> n >> m >> mod;
    len = 1;
    while (len < n + m + 1) len <<= 1;
    // cerr << len << endl;
    int l = __lg(len);
    rep (i, 1, len - 1)
        rev[i] = rev[i >> 1] >> 1 | ((i & 1) << l - 1);
    //     cerr << rev[i] << " ";
    // cerr << endl;
    rep (i, 0, n) {
        int x;
        cin >> x;
        x %= mod;
        a1[i] = x & (1 << 16) - 1;
        a0[i] = x >> 16;
    }
    rep (i, 0, m) {
        int x;
        cin >> x;
        x %= mod;
        b1[i] = x & (1 << 16) - 1;
        b0[i] = x >> 16;
    }

    doubleFFT(a0, a1, len, 1);
    doubleFFT(b0, b1, len, 1);
    rep (i, 0, len - 1) {
        p[i] = a0[i] * b0[i] + I * a1[i] * b0[i];
        q[i] = a0[i] * b1[i] + I * a1[i] * b1[i];
        // cerr << p[i].r << " " << p[i].i << " " << q[i].r << " " << q[i].i << endl;
    }
    // FLG

    FFT(p, len, -1);
    FFT(q, len, -1);

    int T = (1 << 16);
    rep (i, 0, n + m) {
        // cerr << p[i].r << " " << p[i].i << " " << q[i].r << " " << q[i].i << endl;
        ans[i] = T * T * (((int)round(p[i].r)) % mod) % mod + T * (((int)round(p[i].i)) % mod + ((int)round(q[i].r)) % mod) % mod + ((int)round(q[i].i)) % mod;
        ans[i] %= mod;
    }

    rep (i, 0, n + m)
        cout << ans[i] << " ";
    cout << endl;

    return 0;
}