kkogoro

多项式全家桶
包含NTT模数NTT,求导,积分,多项式逆元,多项式对数函数( $\ln A$ ),多项式指数函数(exp),多项...
扫描右侧二维码阅读全文
25
2019/02

多项式全家桶

包含NTT模数NTT,求导,积分,多项式逆元,多项式对数函数( $\ln A$ ),多项式指数函数(exp),多项式开方,多项式除法,分治FFT。
文末有FFT和任意模数NTT,要注意如果保留整数输出系数需要四舍五入。

自带巨大常数

inline int FST(int base, int times) {
    int ret = 1;
    while (times) {
        if (times & 1) ret = (LL)ret * base % MOD;
        times >>= 1;
        base = (LL)base * base % MOD;
    }
    return ret;
}
inline int pls(int a, const int &b) {
    a += b;
    return a >= MOD ? a - MOD : a;
}

namespace POLY {
    int r[MAXN];
    inline void reset(int *a, int n) {
        for (int *end = a + n; a != end; ++a) *a = 0;
        return;
    }
    inline void copy(int *to, int *from, int n) {
        for (int *end = to + n; to != end; ++to, ++from) *to = *from;
        return;
    }
    inline void ButterflyChange(int n, int L) {
        for (int i = 0; i < n; ++i)
            r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
        return;
    }
    inline void NTT(int *a, int n, int type) {
        for (int i = 0; i < n; ++i)
            if (i < r[i]) std::swap(a[i], a[r[i]]);
        for (int mid = 1; mid < n; mid <<= 1) {
            int Xi = FST(type == 1 ? G : inv_G, (MOD - 1) / (mid << 1));
            for (int R = mid << 1, j = 0; j < n; j += R) {
                int w = 1;
                for (int k = 0; k < mid; ++k, w = (LL)w * Xi % MOD) {
                    int x = a[k + j], y = (LL)w * a[k + j + mid] % MOD;
                    a[k + j] = (x + y) % MOD;
                    a[k + j + mid] = ((x - y) % MOD + MOD) % MOD;
                }
            }
        }
        return;
    }
    inline void DFT(int *a, int n) {
        NTT(a, n, 1);
        return;
    }
    inline void IDFT(int *a, int n) {
        NTT(a, n, -1);
        int inv = FST(n, MOD - 2);
        for (int i = 0; i < n; ++i) a[i] = (LL)a[i] * inv % MOD;
        return;
    }
    inline int Mul(int *a, int n, int *b, int m, int *ans) {
        static int tmp1[MAXN], tmp2[MAXN];
        int length = 1, L = 0;
        while (length < (n + m)) length <<= 1, ++L;
        copy(tmp1, a,  n);
        copy(tmp2, b,  m);
        ButterflyChange(length, L);
        DFT(tmp1, length);
        DFT(tmp2, length);
        for (int i = 0; i < length; ++i) ans[i] = (LL)tmp1[i] * tmp2[i] % MOD;
        IDFT(ans, length);
        reset(ans + n + m, length - n - m);
        reset(tmp1, length);
        reset(tmp2, length);
        return length;
    }
    inline void Direv(int *a, int n) {
        for (int i = 1; i <= n; ++i) a[i - 1] = (LL)a[i] * i % MOD;
        a[n] = 0;
        return;
    }
    inline void Inter(int *a, int n) {
        for (int i = n; i >= 1; --i) a[i] = (LL)a[i - 1] * FST(i, MOD - 2) % MOD;
        a[0] = 0;
        return;
    }
    inline void Inv(int deg, int *a, int *b) {
        static int tmp[MAXN];
        if (deg == 1) return void(b[0] = FST(a[0], MOD - 2));
        int length = 1, L = 0, half = (deg + 1) >> 1;
        while (length < (deg << 1)) length <<= 1, ++L;
        Inv(half, a, b);
        reset(b + half, length - half);
        copy(tmp, a,  deg);
        reset(tmp + deg, length - deg);
        ButterflyChange(length, L);
        DFT(tmp, length);
        DFT(b, length);
        for (int i = 0; i < length; ++i) b[i] = (LL)(2 - (LL)tmp[i] * b[i] % MOD + MOD) % MOD * b[i] % MOD;
        IDFT(b, length);
        reset(b + deg, length - deg);
        return;
    }
    inline void Ln(int deg, int *a, int *b) {
        static int tmp[MAXN];
        int length = 1, L = 0;
        while (length < (deg << 1)) length <<= 1, ++L;
        Inv(deg, a, tmp);
        reset(tmp + deg, length - deg);
        copy(b, a,  deg);
        Direv(b, deg);
        ButterflyChange(length, L);
        DFT(b, length);
        DFT(tmp, length);
        for (int i = 0; i < length; ++i) b[i] = (LL)b[i] *  tmp[i] % MOD;
        IDFT(b, length);
        Inter(b, deg);
        reset(b + deg, length - deg);
        reset(tmp, length);
        return;
    }
    void Exp(int deg, int *a, int *b) {
        static int tmp[MAXN];
        if (deg == 1) return void(b[0] = 1);
        int length = 1, L = 0, half = (deg + 1) >> 1;
        while (length < (deg << 1)) length <<= 1, ++L;
        Exp(half, a, b);
        reset(b + half, length - half);
        Ln(deg, b, tmp);
        for (int i = 0; i < deg; ++i) tmp[i] = pls(a[i], MOD - tmp[i]);
        if (++tmp[0] >= MOD) tmp[0] -= MOD;
        ButterflyChange(length, L);
        DFT(tmp, length);
        DFT(b, length);
        for (int i = 0; i < length; ++i) b[i] = (LL)b[i] * tmp[i] % MOD;
        IDFT(b, length);
        reset(b + deg, length - deg);
        reset(tmp, length);
        return;
    }
    inline void Sqrt(int deg, int *a, int *b) {
        static int tmp1[MAXN], tmp2[MAXN];
        if (deg == 1) return void(b[0] = 1);
        int length = 1, L = 0, half = (deg + 1) >> 1;
        while (length < (deg << 1)) length <<= 1, ++L;
        Sqrt(half, a, b);
        reset(b + half, length - half);
        copy(tmp2, a,  deg);
        reset(tmp2 + deg, length - deg);
        Inv(deg, b, tmp1);
        ButterflyChange(length, L);
        DFT(tmp1, length);
        DFT(tmp2, length);
        for (int i = 0; i < length; ++i) tmp1[i] = (LL)tmp1[i] * tmp2[i] % MOD;
        IDFT(tmp1, length);
        int inv2 = FST(2, MOD - 2);
        for (int i = 0; i < deg; ++i) b[i] = (LL)(tmp1[i] + b[i]) % MOD * inv2 % MOD;
        reset(b + deg, length - deg);
        reset(tmp1, length);
        reset(tmp2, length);
        return;
    }
    inline void Div(int *a, int n, int *b, int m, int *c, int *r) {
        static int tmp1[MAXN], tmp2[MAXN];
        int length = 1, L = 0, t = n - m + 1;
        while (length < (t << 1)) length <<= 1, ++L;
        for (int i = 0; i < m; ++i) tmp1[m - i - 1] = b[i];
        Inv(t, tmp1, tmp2); //\pmod{x^t}
        for (int i = 0; i < n; ++i) tmp1[n - i - 1] = a[i];
        reset(tmp1 + t, length - t); //mod x^t
        ButterflyChange(length, L);
        DFT(tmp1, length);
        DFT(tmp2, length);
        for (int i = 0; i < length; ++i) tmp1[i] = (LL)tmp1[i] * tmp2[i] % MOD;
        IDFT(tmp1, length);
        for (int i = 0; i < t; ++i) c[i] = tmp1[t - i - 1];
        length = 1, L = 0;
        while (length < n) length <<= 1, ++L;
        for (int i = 0; i < t; ++i) tmp1[i] = c[i];
        for (int i = 0; i < m; ++i) tmp2[i] = b[i];
        reset(tmp1 + t, length - t);
        reset(tmp2 + m, length - m);
        ButterflyChange(length, L);
        DFT(tmp1, length);
        DFT(tmp2, length);
        for (int i = 0; i < length; ++i) tmp1[i] = (LL)tmp1[i] * tmp2[i] % MOD;
        IDFT(tmp1, length);
        for (int i = 0; i < m - 1; ++i) r[i] = ((LL)(a[i] - tmp1[i]) % MOD + MOD) % MOD;
        reset(tmp1, length);
        reset(tmp2, length);
        return;
    }
    void cdq(int l, int r, int L, int *g, int *ans) {
        static int tmp1[MAXN << 1], tmp2[MAXN << 1];
        if (L <= 0) return;
        cdq(l, (l + r) >> 1, L - 1, g, ans);
        int mid = (r - l) >> 1;
        for (int i = 0; i < mid; ++i) tmp1[i] = ans[i + l];
        for (int i = mid; i < r - l; ++i) tmp1[i] = 0;
        for (int i = 0; i < r - l; ++i) tmp2[i] = g[i];
        ButterflyChange(1 << L, L);
        DFT(tmp1, 1 << L);
        DFT(tmp2, 1 << L);
        for (int i = 0; i < r - l; ++i) tmp1[i] = (LL)tmp1[i] * tmp2[i] % MOD;
        IDFT(tmp1, 1 << L);
        for (int i = mid; i < r - l; ++i) ans[l + i] = ((LL)ans[l + i] + tmp1[i]) % MOD;
        cdq((l + r) >> 1, r, L - 1, g, ans);
        return;
    }
    void cdqNTT(int *g, int n, int *ans) {
        ans[0] = 1;
        for (int i = 1; i < n; ++i) ans[i] = 0;
        int length = 1, L = 0;
        while (length < n) length <<= 1, ++L;
        cdq(0, length, L, g, ans);
    }
}

普通FFT:

const double Pi = acos(-1.0);
struct complex {
    double re, im;
    complex (const double &real_n = 0, const double &imag_n = 0) {re = real_n, im = imag_n;}
} A[MAXN << 1], B[MAXN << 1], C[MAXN << 1];
complex operator + (const complex &a, const complex &b) {return complex(a.re + b.re, a.im + b.im);}
complex operator - (const complex &a, const complex &b) {return complex(a.re - b.re, a.im - b.im);};
complex operator * (const complex &a, const complex &b) {return complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);}
int r[MAXN << 1];
inline void ButterflyChange(int n, int L) {
    for (int i = 0; i < n; ++i) 
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
    return;
}
inline void FFT(complex *a, int n, int type) {
    for (int i = 0; i < n; ++i) 
        if (i < r[i]) std::swap(a[i], a[ r[i] ]);
    for (int mid = 1; mid < n; mid <<= 1) {
        complex Wn( cos(Pi / mid), type * sin(Pi / mid));
        for (int R = mid << 1, j = 0; j < n; j += R) {
            complex w(1, 0);
            for (int k = 0; k < mid; ++k, w = w * Wn) {
                complex x = a[j + k], y = w * a[j + mid + k];
                a[j + k] = x + y;
                a[j + mid + k] = x - y;
            }
        }
    }
    return;
}
inline void DFT(complex *a, int n) {
    FFT(a, n, 1);
    return;
}
inline void IDFT(complex *a, int n) {
    FFT(a, n, -1);
    for (int i = 0; i < n; ++i) a[i].re /= (double)n;
    return;
}
inline int multiply(complex *a, int n, complex *b, int m, complex *ans) {
    static complex tmp1[MAXN << 1], tmp2[MAXN << 1];
    int length = 1, L = 0;
    while (length < n + m) length <<= 1, ++L;
    for (int i = 0; i < n; ++i) tmp1[i] = a[i];
    for (int i = 0; i < m; ++i) tmp2[i] = b[i];
    ButterflyChange(length, L);
    DFT(tmp1, length);
    DFT(tmp2, length);
    for (int i = 0; i < length; ++i) ans[i] = tmp1[i] * tmp2[i];
    IDFT(ans, length);
    for (int i = n + m; i < length; ++i) ans[i].re = ans[i].im = 0;
    for (int i = 0; i < length; ++i) tmp1[i].re = tmp1[i].im = tmp2[i].re = tmp2[i].im = 0;
    return length;
}

任意模数NTT

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>


inline int read () {
    int x = 0, f = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9') {
        if (ch == '-') f = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9') {
        x = x * 10 + ch - '0';
        ch = getchar();
    }
    return x * f;
}


typedef long long LL;
const int MAXN = (100007 + 100007) << 1;


LL M[3] = {167772161, 469762049, 998244353}, G = 3;


inline LL FST(LL base, int times, LL mod) {
    LL ret = 1;
    while (times) {
        if (times & 1) ret = ret * base % mod;
        times >>= 1;
        base = base * base % mod;
    }
    return ret;
}
inline LL fast_mul(const LL &x, const LL &y, const LL &mod) {
    LL tmp = (x * y - (LL)((long double)x / mod * y + 1.0e-8) * mod);
    return tmp < 0 ? tmp + mod : tmp;
}

namespace POLY {
    int r[MAXN];
    inline void reset(LL *a, int n) {
        for (LL *end = a + n; a != end; ++a) *a = 0;
        return;
    }
    inline void copy(LL *to, LL *from, int n) {
        for (LL *end = to + n; to != end; ++to, ++from) *to = *from;
        return;
    }
    inline void ButterflyChange(int length, int L) {
        for (int i = 0; i < length; ++i) 
            r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
        return;
    }
    inline void NTT(LL *a, int n, int type, int mod) {
        LL MOD = M[mod];
        for (int i = 0; i < n; ++i)
            if (i < r[i]) std::swap(a[i], a[ r[i] ]);
        for (int mid = 1; mid < n; mid <<= 1) {
            LL Xi = FST(type == 1 ? G : FST(G, (int)MOD - 2, MOD), (int)(MOD - 1) / (mid << 1), MOD);
            for (int R = mid << 1, j = 0; j < n; j += R) {
                LL w = 1;
                for (int k = 0; k < mid; ++k, w = w * Xi % MOD) {
                    LL x = a[k + j], y = w * a[k + j + mid] % MOD;
                    a[k + j] = (x + y) % MOD;
                    a[k + j + mid] = ((x - y) % MOD + MOD) % MOD;
                }
            }
        }
        return;
    }
    inline void DFT(LL *a, int n, int mod) {
        NTT(a, n, 1, mod);
        return;
    }
    inline void IDFT(LL *a, int n, int mod) {
        NTT(a, n, -1, mod);
        LL inv = FST(n, (int)M[mod] - 2, M[mod]);
        for (int i = 0; i < n; ++i) a[i] = a[i] * inv % M[mod];
        return;
    }
    inline int MTT(LL *a, int n, LL *b, int m, LL p, LL *ans) {
        static LL tmp[4][MAXN];
        int length = 1, L = 0;
        while (length < n + m) length <<= 1, ++L;
        ButterflyChange(length, L);
        for (int k = 0; k < 3; ++k) {
            copy(tmp[k], a, n);
            copy(tmp[k + 1], b, m);
            DFT(tmp[k], length, k);
            DFT(tmp[k + 1], length, k);
            for (int i = 0; i < length; ++i) tmp[k][i] = tmp[k][i] * tmp[k + 1][i] % M[k];
            IDFT(tmp[k], length, k);
            reset(tmp[k] + n + m, length - n - m);
            reset(tmp[k + 1], length);
        }
        for (int i = 0; i < n + m - 1; ++i) {
            LL preM = M[0] * M[1];
            LL A = (fast_mul(tmp[0][i] * M[1] % preM, FST(M[1] % M[0], (int)M[0] - 2, M[0]), preM) +
                    fast_mul(tmp[1][i] * M[0] % preM, FST(M[0] % M[1], (int)M[1] - 2, M[1]), preM)) % preM;
            LL k = ((tmp[2][i] - A) % M[2] + M[2]) % M[2] * FST(preM % M[2], (int)M[2] - 2, M[2]) % M[2];
            ans[i] = ((k % p) * (preM % p) % p + A) % p;
        }
        for (int k = 0; k < 3; ++k) reset(tmp[k], length);
        return length;
    }
}


LL A[MAXN], B[MAXN], C[MAXN];


int main() {        
    int n = read() + 1;
    int m = read() + 1;
    int p = read();
    for (int i = 0; i < n; ++i) A[i] = (LL)read();
    for (int i = 0; i < m; ++i) B[i] = (LL)read();
    POLY::MTT(A, n, B, m, p, C);
    for (int i = 0; i < n + m - 1; ++i)
        printf("%d ", (int)C[i]);
    return 0;
}
Last modification:March 18th, 2019 at 06:12 pm

One comment

  1. CalvinJin

    神仙啊qwq

Leave a Comment