[Luogu4245] 【模板】任意模数NTT [三模数NTT][拆系数FFT]

Link
Luogu - https://www.luogu.org/problemnew/show/P4245

---

我并没有要讲的意思 贴一个仙人板板

---

#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstring>
using namespace std;
char frBB[1<<12],*frS=frBB,*frT=frBB;
#define getchar() (frS==frT&&(frT=(frS=frBB)+fread(frBB,1,1<<12,stdin),frS==frT)?EOF:*frS++)
template<typename T>inline void read(T&x)
{
    x=0;bool w=0;char ch=getchar();
    while(!isdigit(ch))w|=(ch=='-'),ch=getchar();
    while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
    w?(x=-x):0;
}
#define R register
const int MAXN = 3e5+5;
const double tpi = acos(-1.0);
int P;
long long lP;
namespace Poly
{
    struct Complex
    {
        double r, i;
        inline Complex(const double& Re = 0, const double& Im = 0) { r = Re; i = Im;}
        inline Complex operator + (const Complex& o) const { return Complex(r + o.r, i + o.i);}
        inline Complex operator - (const Complex& o) const { return Complex(r - o.r, i - o.i);}
        inline Complex operator * (const Complex& o) const { return Complex(r * o.r - i * o.i, r * o.i + i * o.r);}
        inline Complex operator * (const double& x) const { return Complex(r * x, i * x);}
        inline Complex operator ~ () const { return Complex(r, -i);}
        inline void operator += (const Complex& o) { r += o.r; i += o.i;}
        inline void operator -= (const Complex& o) { r -= o.r; i -= o.i;}
        inline void operator *= (const Complex& o) { *this = *this * o;}
        inline void operator *= (const double& x) { r = r * x; i = i * x;}
    }A[MAXN], B[MAXN], C[MAXN], D[MAXN], Wn[MAXN];
    int Calc(const int& x)
    {
        register int t = 1;
        while (t < x) t <<= 1;
        return t;
    }
    void PreWn(const int& n)
    {
        for (R int i = 1; i < n; i <<= 1)
        {
            Wn[i] = Complex(1, 0);
            for (R int j = 1; j < i; ++j)
            {
                Wn[i + j] = ((j & 31) == 1) ? Complex(cos(tpi * j / i), sin(tpi * j / i)) : Wn[i + j - 1] * Wn[i + 1];
            }
        }
    }
    void FFT(Complex* F, const int& n)
    {
        for (R int j = 0, i = 0; i < n; ++i)
        {
            if (i > j) swap(F[i], F[j]);
            for (R int k = n >> 1; (j ^= k) < k; k >>= 1);
        }
        R Complex Y;
        for (R int Cur, Len, Mid = 1; Mid < n; Mid <<= 1)
        {
            Len = Mid << 1;
            for (R int Pos = 0; Pos < n; Pos += Len)
            {
                Cur = Pos;
                for (R int Sub = 0; Sub < Mid; ++Sub, ++Cur)
                {
                    Y = Wn[Mid + Sub] * F[Cur + Mid]; // omega_n^k\cdot F(n^{\frac{n}{2}+k})
                    F[Cur + Mid] = F[Cur] - Y;
                    F[Cur] += Y;
                }
            }
        }
    }
    void DFT(Complex* F, int n)
    {
        if (n == 1) return;
        static Complex P[MAXN>>1];
        n >>= 1;
        for (R int i = 0; i < n; ++i) P[i].r = F[i<<1].r, P[i].i = F[i<<1|1].r;
        FFT(P, n);
        R Complex Q, X, Y;
        for (R int i = 0; i < n; ++i)
        {
            Q = ~P[(n-i)&(n-1)];
            X = (P[i] + Q) * .5; // * .5 = / 2
            Y = (P[i] - Q) * Complex(0, -.5) * Wn[n + i]; // * -.5i = / 2i
            F[i] = X + Y;
            F[n + i] = X - Y; 
        }
    }
    void IDFT(Complex* F, int n)
    {
        if (n == 1) return;
        static Complex P[MAXN>>1];
        reverse(F + 1, F + n); n >>= 1;
        for (R int i = 0; i < n; ++i) P[i] = (F[i] + F[n + i]) * .5 + (F[i] - F[n + i]) * Complex(0, .5) * Wn[n + i];
        FFT(P, n);
        R double t = 1. / n;
        for (R int i = 0; i < n; ++i) F[i<<1] = Complex(P[i].r * t, 0), F[i<<1|1] = Complex(P[i].i * t, 0);
    }
    inline long long Estimate(const double& x) { return (long long)(x + .5)%lP;}
    inline long long Adjust(long long x)
    {
        x = x % lP + lP;
        return (x >= lP) ? (x - lP) : x;
    }
    void Multi(int* a, int* b, const int& n, const int& m)
    {
        R int l = Calc(n + m);
        for (R int i = 0; i <= n; ++i) A[i].r = a[i] >> 15, B[i].r = a[i] & 32767;
        for (R int i = 0; i <= m; ++i) C[i].r = b[i] >> 15, D[i].r = b[i] & 32767;
        DFT(A, l), DFT(B, l), DFT(C, l), DFT(D, l);
        R Complex t;
        for (R int i = 0; i < l; ++i)
        {
            t = A[i];
            A[i] = A[i] * D[i] + B[i] * C[i];
            B[i] *= D[i];
            C[i] *= t;
        }
        IDFT(A, l); IDFT(B, l); IDFT(C, l);
        for (R int i = 0; i <= n + m; ++i) a[i] = Adjust((Estimate(C[i].r) << 30) + (Estimate(A[i].r) << 15) + Estimate(B[i].r));
    }
}
int n, m, l, F[MAXN], G[MAXN];
int main()
{
    read(n), read(m), read(P); lP = P;;
    l = Poly::Calc(n + m);
    Poly::PreWn(l);
    for (R int i = 0; i <= n; ++i) read(F[i]);
    for (R int i = 0; i <= m; ++i) read(G[i]);
    Poly::Multi(F, G, n, m);
    for (R int i = 0; i <= n + m; ++i) printf("%d ", F[i]);
    return 0;
}