题目地址：洛谷：【P4245】【模板】MTT – 洛谷

题目描述

多项式乘法。

输入输出样例

输入样例#1：

5 8 28
19 32 0 182 99 95
77 54 15 3 98 66 21 20 38

输出样例#1：

7 18 25 19 5 13 12 2 9 22 5 27 6 26

说明

数据范围1e5。

题解

由于NTT对模数有要求，我们需要使用任意模数卷积算法，也就是俗称的MTT。具体做法是，把一个数拆成32768x+y的形式，对拆出来的四组数分别卷积后拼起来。这样可以避免计算过程中的越界问题，一般的数据能够满足精度。
单位根最好预处理出来，这样才不会有精度问题。

代码

// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 1 << 20;
const double PI = acos(-1);
int MO;

struct Complex {
    double real, imag;
    Complex(double real = 0, double imag = 0) : real(real), imag(imag) {}
    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;
    }
    friend inline Complex operator*(double x, Complex cp) {
        return Complex(x * cp.real, x * cp.imag);
    }
    inline Complex operator/(double x) const {
        return Complex(real / x, imag / x);
    }
    inline Complex operator/=(double x) {
        return (*this) = (*this) / x;
    }
    friend inline Complex operator/(double 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 double length() {
        return sqrt(real * real + imag * imag);
    }
};

int n, m, len, rev[MAXN];
LL x[MAXN], y[MAXN], z[MAXN];
Complex a[MAXN], b[MAXN], c[MAXN], d[MAXN], e[MAXN], f[MAXN], g[MAXN], h[MAXN], wn[MAXN];

inline void fft(Complex *arr, int f) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) std::swap(arr[i], arr[rev[i]]);
    }
    for(int i = 1; i < n; i <<= 1) {
        for(int j = 0; j < n; j += i << 1) {
            for(int k = 0; k < i; k++) {
                Complex w = Complex(wn[n / i * k].real, f * wn[n / i * k].imag), 
                    x = arr[j + k], y = w * arr[j + k + i];
                arr[j + k] = x + y;
                arr[j + k + i] = x - y;
            }
        }
    }
    if(f == -1) {
        for(int i = 0; i < n; i++) {
            arr[i] /= n;
        }
    }
}

int main() {
    n = readint(); m = readint(); MO = readint();
    for(int i = 0; i <= n; i++) {
        x[i] = readint() % MO;
    }
    for(int i = 0; i <= m; i++) {
        y[i] = readint() % MO;
    }
    m += n;
    for(n = 1; n <= m; n <<= 1) len++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
    for(int i = 1; i < n; i <<= 1) {
        for(int k = 0; k < i; k++) {
            wn[n / i * k] = Complex {cos(PI * k / i), sin(PI * k / i)};
        }
    }
    for(int i = 0; i < n; i++) {
        a[i].real = x[i] >> 15; b[i].real = x[i] & 0x7fff;
        c[i].real = y[i] >> 15; d[i].real = y[i] & 0x7fff;
    }
    fft(a, 1); fft(b, 1); fft(c, 1); fft(d, 1);
    for(int i = 0; i < n; 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(e, -1); fft(f, -1); fft(g, -1); fft(h, -1);
    for(int i = 0; i < n; i++) {
        z[i] = (((LL(round(e[i].real)) % MO) << 30) % MO 
            + ((LL(round(f[i].real)) % MO) << 15) % MO
            + ((LL(round(g[i].real)) % MO) << 15) % MO 
            + LL(round(h[i].real)) % MO) % MO;
    }
    for(int i = 0; i <= m; i++) {
        printf("%lld ", z[i]);
    }
    return 0;
}