题解 P4245 【【模板】MTT】

2018-03-07 08:29:37


其实这个题卡精度来着……我也没想到要开long double才能过……因为数据是用ddd的模板做的,他用的好像是合并dft精度比较高来着QAQ
其实就是把每个数拆成 $a \cdot 32768 + b$ 的形式,然后乘起来再展开就可以了。当然你想要NTT合并也行qwq
代码如下:

#include <bits/stdc++.h>

typedef long long int64;
typedef long double D;

int MOD;
const int MAXN = 524288;
const D PI = acos(-1);

struct complex {
    D real, imag;
    complex() { real = imag = 0; }
    complex(D x): real(x), imag(0) {} 
    complex(D x, D y): real(x), imag(y) {}
    inline complex conj() { return complex(real, -imag); }
    inline complex operator+(complex rhs) const { return complex(real + rhs.real, imag + rhs.imag); }
    inline complex operator-(complex rhs) const { return complex(real - rhs.real, imag - rhs.imag); }
    inline complex operator*(complex rhs) const { return complex(real * rhs.real - imag * rhs.imag,
                                                                 imag * rhs.real + real * rhs.imag); }
    inline complex operator*=(complex rhs) { return (*this) = (*this) * rhs; }
    //(a+bi)(c+di) = (ac-bd) + (bc+ad)i
    friend inline complex operator*(D x, complex cp) { return complex(x * cp.real, x * cp.imag); }
    inline complex operator/(D x) const { return complex(real / x, imag / x); }
    inline complex operator/=(D x) { return (*this) = (*this) / x; }
    friend inline complex operator/(D x, complex cp) { return x * cp.conj() / (cp.real * cp.real - cp.imag * cp.imag); }
    inline complex operator/(complex rhs) const { 
        return (*this) * rhs.conj() / (rhs.real * rhs.real - rhs.imag * rhs.imag); 
    }
    inline complex operator/=(complex rhs) { return (*this) = (*this) / rhs; }
    inline D length() { return sqrt(real * real + imag * imag); }
};

inline complex get_omega(int len, bool inv) {
    return inv ? complex(std::cos(2*PI / len), -std::sin(2*PI / len))
                :complex(std::cos(2*PI / len),  std::sin(2*PI / len));
}

inline void fft(int len, complex* A, bool inv = false) {
    static int R[MAXN+10];
    for(int i = 0; i < len; i++)
        R[i] = ((R[i>>1]>>1) | (len >> (i&1))) & (len-1);
    for(int i = 0; i < len; i++)
        if(R[i] > i) std::swap(A[i], A[R[i]]);
    for(int step = 1; step < len; step <<= 1) 
        for(int i = 0; i < len; i += (step<<1)) {
            complex omega = 1, t = 0, rt = get_omega(step<<1, inv);
            for(int j = 0; j < step; j++, omega *= rt) {
                t = A[i+j+step] * omega;
                A[i+j+step] = A[i+j] - t;
                A[i+j] = A[i+j] + t;
            }
        }
    if(inv)
        for(int i = 0; i < len; i++)
            A[i] /= len;
}

void mtt(int deg, int *lhs, int *rhs, int *ret) {
    static complex A[MAXN+10], B[MAXN+10], C[MAXN+10], D[MAXN+10],
                   E[MAXN+10], F[MAXN+10], G[MAXN+10], H[MAXN+10];
    int len; for(len = 1; len <= 2*deg; len<<=1);
    for(int i = 0; i < len; i++) {
        lhs[i] %= MOD; rhs[i] %= MOD;
        A[i] = lhs[i] >> 15; B[i] = lhs[i] & 0x7fff;
        C[i] = rhs[i] >> 15; D[i] = rhs[i] & 0x7fff;
    }
    fft(len, A); fft(len, B); fft(len, C); fft(len, D);
    for(int i = 0; i < len; i++) {
        E[i] = A[i] * C[i]; F[i] = B[i] * C[i];
        G[i] = A[i] * D[i]; H[i] = B[i] * D[i];
    }
    fft(len, E, true); fft(len, F, true); 
    fft(len, G, true); fft(len, H, true);
    for(int i = 0; i < len; i++)
        ret[i] = (((int64(round(E[i].real)) % MOD)<<30) % MOD + ((int64(round(F[i].real)) % MOD)<<15) % MOD
            + ((int64(round(G[i].real)) % MOD)<<15) % MOD + int64(round(H[i].real)) % MOD) % MOD;
}

int deg, n, m, F[MAXN+10], G[MAXN+10], H[MAXN+10];

inline int readint() {
    int f=1, r=0; char c=getchar();
    while(!isdigit(c)) { if(c=='-')f=-1; c=getchar(); }
    while(isdigit(c)) { r=r*10+c-'0'; c=getchar(); }
    return f*r;
}

int main() {
    n = readint(); m = readint(); MOD = readint();
    for(int i = 0; i <= n; i++) 
        F[i] = readint();
    for(int i = 0; i <= m; i++)
        G[i] = readint();
    mtt(std::max(n, m), F, G, H);
    for(int i = 0; i <= n + m; i++)
        printf("%d ", H[i]);
}