数论

带模数的模板要注意是否使用了\(long ~~ long\)类。

ModInt

template <unsigned M_> struct ModInt {
    static constexpr unsigned M = M_;
    unsigned x;
    constexpr ModInt() : x(0U) {}
    constexpr ModInt(unsigned x_) : x(x_ % M) {}
    constexpr ModInt(unsigned long long x_) : x(x_ % M) {}
    constexpr ModInt(int x_) : x(((x_ %= static_cast<int>(M)) < 0) ? (x_ + static_cast<int>(M)) : x_) {}
    constexpr ModInt(long long x_) : x(((x_ %= static_cast<long long>(M)) < 0) ? (x_ + static_cast<long long>(M)) : x_) {}
    ModInt& operator +=(const ModInt& a) {
        x = ((x += a.x) >= M) ? (x - M) : x;
        return *this;
    }
    ModInt& operator -=(const ModInt& a) {
        x = ((x -= a.x) >= M) ? (x + M) : x;
        return *this;
    }
    ModInt& operator *=(const ModInt& a) {
        x = (static_cast<unsigned long long>(x) * a.x) % M;
        return *this;
    }
    ModInt& operator /=(const ModInt& a) { return (*this *= a.inv()); }
    ModInt pow(long long e) const {
        if (e < 0) {
            return inv().pow(-e);
        }
        ModInt a = *this, b = 1U;
        for ( ; e; e >>= 1) {
            if (e & 1) {
                b *= a;
            }
            a *= a;
        }
        return b;
    }
    ModInt inv() const {
        unsigned a = M, b = x;
        int y = 0, z = 1;
        while (b) {
            const unsigned q = a / b, c = a - q * b;
            a = b;
            b = c;
            const int w = y - static_cast<int>(q) * z;
            y = z;
            z = w;
        }
        assert(a == 1U);
        return ModInt(y);
    }
    ModInt operator +() const { return *this; }
    ModInt operator -() const {
        ModInt a;
        a.x = x ? (M - x) : 0U;
        return a;
    }
    ModInt operator +(const ModInt& a) const { return (ModInt(*this) += a); }
    ModInt operator -(const ModInt& a) const { return (ModInt(*this) -= a); }
    ModInt operator *(const ModInt& a) const { return (ModInt(*this) *= a); }
    ModInt operator /(const ModInt& a) const { return (ModInt(*this) /= a); }
    template <class T> friend ModInt operator +(T a, const ModInt& b) { return (ModInt(a) += b); }
    template <class T> friend ModInt operator -(T a, const ModInt& b) { return (ModInt(a) -= b); }
    template <class T> friend ModInt operator *(T a, const ModInt& b) { return (ModInt(a) *= b); }
    template <class T> friend ModInt operator /(T a, const ModInt& b) { return (ModInt(a) /= b); }
    explicit operator bool() const { return x; }
    bool operator ==(const ModInt& a) const { return (x == a.x); }
    bool operator !=(const ModInt& a) const { return (x != a.x); }
    friend std::ostream& operator <<(std::ostream& os, const ModInt& a) { return os << a.x; }
};
using mint = ModInt<998244353>;

判断素数

bool isPrime(int n) {
    if (n == 2 || n == 3) {
        return true;
    }
    if ((n % 6 != 1 && n % 6 != 5) || n == 1) {
        return false;
    }
    for (int i = 5; i * i <= n; i += 6) {
        if (n % i == 0 || n % (i + 2) == 0) { 
            return false;
        }
    }
    return true;
}

快速幂

不取模

int power(int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1) {
            res *= a;
        }
        a *= a;
        n >>= 1;
    }
    return res;
}

取模

int power(int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1) {
            res = (res * a) % MOD;
        }
        a = (a * a) % MOD;
        n >>= 1;
    }
    return res;
}

矩阵

const int R = 3;
struct Matrix {
    int mat[R + 1][R + 1];
    inline Matrix() {
        memset(mat, 0, sizeof(mat));
    }
    inline Matrix(vector<vector<int>> M) {
        for (int i = 1; i <= R; ++i) {
            for (int j = 1; j <= R; ++j) {
                mat[i][j] = M[i - 1][j - 1];
            }
        }
    }
    Matrix operator *(Matrix x) {
        Matrix res;
        for (int i = 1; i <= R; ++i) {
            for (int j = 1; j <= R; ++j) {
                for (int k = 1; k <= R; ++k) {
                    // res.mat[i][j] += mat[i][k] * x.mat[k][j];
                    res.mat[i][j] = (res.mat[i][j] + mat[i][k] * x.mat[k][j] % MOD) % MOD;
                }
            }
        }
        return res;
    }
    Matrix operator ^(int n) {
        Matrix res, tmp = *this;
        for (int i = 1; i <= R; ++i) {
            res.mat[i][i] = 1;
        }
        while (n) {
            if (n & 1) {
                res *= tmp;
            }
            tmp *= tmp;
            n >>= 1;
        }
        return res;
    }
    void operator =(Matrix x) {
        for (int i = 1; i <= R; ++i) {
            for (int j = 1; j <= R; ++j) {
                mat[i][j] = x.mat[i][j];
            }
        }
    }
    Matrix operator *=(Matrix x) {
        *this = *this * x;
        return *this;
    }
    Matrix operator ^=(int n) {
        *this = *this ^ n;
        return *this;
    }
};
ostream& operator <<(ostream& os, const Matrix& m) {
    for (int i = 1; i <= R; ++i) {
        os << '[';
        for (int j = 1; j <= R; ++j) {
            os << m.mat[i][j];
            if (j != R) {
                os << ", ";
            }
        }
        os << ']';
        if (i != R) {
            os << '\n';
        }
    }
    return os;
}

逆元

扩展欧几里得

int exgcd(int a, int b, int& x, int& y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int res = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return res;
}
int inv(int a, int b) {
    int x, y; 
    exgcd(a, b, x, y);
    return (b + x % b) % b;
}

费马小定理

int inv(int a) {
    return power(a, MOD - 2);
}

线性递推

int _inv[MAXN];
void getInvs() {
    _inv[1] = 1;
    for (int i = 2; i < MAXN; ++i) {
        _inv[i] = (MOD - MOD / i) * _inv[MOD % i] % MOD;
    }
}

组合数

调用\(C(m, n)\)可以得到\(C_m^n\)的值。

int power(int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1) {
            res = (res * a) % MOD;
        }
        a = (a * a) % MOD;
        n >>= 1;
    }
    return res;
}
int inv(int x) {
    return power(x, MOD - 2);
}
int fact[MAXN], _inv[MAXN];
void init() {
    fact[0] = 1;
    for (int i = 1; i < MAXN; ++i) {
        fact[i] = fact[i - 1] * i % MOD;
    }
    _inv[MAXN - 1] = inv(fact[MAXN - 1]);
    for (int i = MAXN - 2; i >= 0; --i) {
        _inv[i] = _inv[i + 1] * (i + 1) % MOD;
    }
}
// C(m, n) = m! / n!(m - n)!
int C(int m, int n) {
    if (n > m) return 0;
    return fact[m] * _inv[n] % MOD * _inv[m - n] % MOD;
}

中国剩余定理

int exgcd(int a, int b, int& x, int& y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int res = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return res;
}
int mul(int x, int y, int p) {
    return (__int128)x * y % p;
}
int excrt(int n, int a[], int m[]) {
    a[1] %= m[1];
    for (int i = 2; i <= n; i++) {
        a[i] %= m[i];
        int x, y, d = exgcd(m[1], m[i], x, y);
        if ((a[i] - a[1]) % d) return -1;
        int t = m[i] / d;
        x = (mul((a[i] - a[1]) / d, x, t) % t + t) % t;
        a[1] += x * m[1];
        m[1] *= m[i] / d;
        a[1] = (a[1] % m[1] + m[1]) % m[1];
    }
    return a[1];
}

求单个数所有约数

vector<int> divide(int x) {
    vector<int> d;
    for (int i = 1; i * i <= x; ++i) {
        if (x % i == 0) {
            d.push_back(i);
            d.push_back(x / i);
        }
    }
    sort(d.begin(), d.end());
    return d;
}

筛

bool isNotP[MAXN];
int mu[MAXN], phi[MAXN], p[MAXN], cnt;
void sieve() {
    isNotP[1] = phi[1] = mu[1] = 1;
    for (int i = 2; i < MAXN; ++i) {
        if (!isNotP[i]) {
            mu[i] = -1;
            phi[i] = i - 1;
            p[++cnt] = i;
        }
        for (int j = 1; j <= cnt && i * p[j] < MAXN; ++j) {
            isNotP[i * p[j]] = true;
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                phi[i * p[j]] = phi[i] * p[j];
                break;
            } else {
                mu[i * p[j]] = -mu[i];
                phi[i * p[j]] = phi[i] * (p[j] - 1);
            }
        }
    }
}

bool isNotP[MAXN];
int p[MAXN], cnt;
void sieve() {
    isNotP[1] = true;
    for (int i = 2; i < MAXN; ++i) {
        if (!isNotP[i]) {
            p[++cnt] = i;
        }
        for (int j = 1; j <= cnt && i * p[j] < MAXN; ++j) {
            isNotP[i * p[j]] = true;
            if (i % p[j] == 0) {
                break;
            }
        }
    }
}

NTT

namespace NTT {

const int MAXN = 4e6 + 5;
const int MOD = 998244353;
const int G = 3;

int power(int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1) {
            res = (res * a) % MOD;
        }
        a = (a * a) % MOD;
        n >>= 1;
    }
    return res;
}
int inv(int a) {
    return power(a, MOD - 2);
}
struct Complex {
    double x, y;
    Complex(double _x, double _y): x(_x), y(_y) {}
    Complex operator +(Complex oth) {
        return Complex(x + oth.x, y + oth.y);
    }
    Complex operator -(Complex oth) {
        return Complex(x - oth.x, y - oth.y);
    }
    Complex operator *(Complex oth) {
        return Complex(x * oth.x - y * oth.y, x * oth.y + y * oth.x);
    }
};
int R[MAXN], L, limit = 1;
void NTT(int a[], int opt) {
    for (int i = 0; i < limit; ++i) {
        if (i < R[i]) {
            swap(a[i], a[R[i]]);
        }
    }
    for (int mid = 1; mid < limit; mid <<= 1) {
        int val = power(G, (MOD - 1) / (mid * 2));
        if (opt == -1) val = inv(val);
        for (int len = mid << 1, pos = 0; pos < limit; pos += len) {
            for (int k = 0, w = 1; k < mid; ++k, w = w * val % MOD) {
                int x = a[pos + k], y = w * a[pos + mid + k] % MOD;
                a[pos + k] = (x + y) % MOD;
                a[pos + k + mid] = (x - y + MOD) % MOD;
            }
        }
    }
    if (opt == 1) return ;
    int t = inv(limit);
    for (int i = 0; i < limit; ++i) {
        a[i] = a[i] * t % MOD;
    }
}
void poly(int a[], int b[], int deg) {
    while (limit <= deg) {
        limit <<= 1, ++L;
    }
    for (int i = 0; i < limit; ++i) {
        R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    NTT(a, 1);
    NTT(b, 1);
    for (int i = 0; i < limit; ++i) {
        a[i] = a[i] * b[i] % MOD;
    }
    NTT(a, -1);
}

} // namespace NTT
using NTT::poly;

FFT

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

const double PI = acos(-1.0);
const int MAXN = 5e5 + 50;

struct Complex {
    double x, y;
    Complex(double _x = 0, double _y = 0): x(_x), y(_y) {}
    Complex operator +(Complex oth) {
        return Complex(x + oth.x, y + oth.y);
    }
    Complex operator -(Complex oth) {
        return Complex(x - oth.x, y - oth.y);
    }
    Complex operator *(Complex oth) {
        return Complex(x * oth.x - y * oth.y, x * oth.y + y * oth.x);
    }
} a[MAXN], b[MAXN];
int R[MAXN], L, limit = 1;
void FFT(Complex A[], int opt) {
    for (int i = 0; i < limit; ++i) {
        if (i < R[i]) {
            swap(a[i], a[R[i]]);
        }
    }
    for (int mid = 1; mid < limit; mid <<= 1) {
        Complex Wn(cos(PI / mid), opt * sin(PI / mid));
        for (int len = mid << 1, pos = 0; pos < limit; pos += len) {
            Complex w(1, 0);
            for (int k = 0; k < mid; ++k, w = w * Wn) {
                Complex x = A[pos + k];
                Complex y = w * A[pos + mid + k];
                A[pos + k] = x + y;
                A[pos + mid + k] = x - y;
            }
        }
    }
    if (opt == 1) return ;
    for (int i = 0; i <= limit; ++i) {
        a[i].x /= limit;
    }
}
void poly(Complex a[], Complex b[], int deg) {
    while (limit <= deg) {
        limit <<= 1;
        ++L;
    }
    for (int i = 0; i < limit; ++i) {
        R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    FFT(a, 1);
    FFT(b, 1);
    for (int i = 0; i <= limit; ++i) {
        a[i] = a[i] * b[i];
    }
    FFT(a, -1);
}

int main(int argc, char *argv[]) {
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; ++i)  {
        scanf("%lf", &a[i].x);
    }
    for (int i = 0; i <= m; ++i) {
        scanf("%lf", &b[i].x);
    }
    poly(a, b, n + m);
    for (int i = 0; i <= n + m; i++) {
        printf("%d%c", (int)(a[i].x + 0.5), " \n"[i == n + m]);
    }
    return 0;
}

FWT

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

const int MOD = 998244353;
const int inv2 = 499122177;
const int MAXN = (1 << 17) + 50;

int a[MAXN], b[MAXN], c[MAXN];
void FWT_OR(int a[], int len, int opt) {
    for (int i = 1; i < len; i <<= 1) {
        for (int j = 0, R = i << 1; j < len; j += R) {
            for (int k = 0; k < i; ++k) {
                if (opt == 1) {
                    a[i + j + k] = (a[i + j + k] + a[j + k]) % MOD;
                }
                else {
                    a[i + j + k] = (a[i + j + k] - a[j + k] + MOD) % MOD;
                }
            }
        }
    }
}
void FWT_AND(int a[], int len, int opt) {
    for (int i = 1; i < len; i <<= 1) {
        for (int j = 0, R = i << 1; j < len; j += R) {
            for (int k = 0; k < i; ++k) {
                if (opt == 1) {
                    a[j + k] = (a[j + k] + a[i + j + k]) % MOD;
                }
                else {
                    a[j + k] = (a[j + k] - a[i + j + k] + MOD) % MOD;
                }
            }
        }
    }
}
void FWT_XOR(int a[], int len, int opt) {
    for (int i = 1; i < len; i <<= 1) {
        for (int j = 0, R = i << 1; j < len; j += R) {
            for (int k = 0; k < i; ++k) {
            int x = a[j + k], y = a[i + j + k];
            if (opt == 1) {
                a[j + k] = (x + y) % MOD;
                a[i + j + k] = (x - y + MOD) % MOD;
            }
            else {
                a[j + k] = (x + y) % MOD * inv2 % MOD;
                a[i + j + k] = (x - y + MOD) % MOD * inv2 % MOD;
            }
        }
        }
    }
}
// OR, AND, XOR操作同理, len不足2的幂次用0补齐
int main(int argc, char *argv[]) {
    int n;
    scanf("%d", &n);
    int len = 1 << n;
    for (int i = 0; i < len; ++i) {
        scanf("%d", &a[i]);
    }
    for (int i = 0; i < len; ++i) {
        scanf("%d", &b[i]);
    }
    FWT_OR(a, len, 1);
    FWT_OR(b, len, 1);
    for (int i = 0; i < len; ++i) {
        c[i] = a[i] * b[i] % MOD;
    }
    FWT_OR(c, len, -1);
    for (int i = 0; i < len; ++i) {
        printf("%d%c", c[i], " \n"[i == len - 1]);
    }
    FWT_OR(a, len, -1);
    FWT_OR(b, len, -1);
    return 0;
}

杜教BM

namespace BM {

const int MAXN = 1e5 + 5;
const int MOD = 998244353;

vector<int> Md;
int res[MAXN], base[MAXN], _c[MAXN], _md[MAXN];

int power(int a, int n) {
    int res = 1;
    while (n) {
        if (n & 1) {
            res = (res * a) % MOD;
        }
        a = (a * a) % MOD;
        n >>= 1;
    }
    return res;
}
void mul(int a[], int b[], int k) {
    for (int i = 0; i < k + k; ++i) {
        _c[i] = 0;
    }
    for (int i = 0; i < k; ++i) {
        if (a[i]) {
            for (int j = 0; j < k; ++j) {
                _c[i + j] = (_c[i + j] + a[i] * b[j]) % MOD;
            }
        }
    }
    for (int i = k + k - 1; i >= k; --i) {
        if (_c[i]) {
            for (int j = 0; j < Md.size(); ++j) {
                _c[i - k + Md[j]] = (_c[i - k + Md[j]] - _c[i] * _md[Md[j]]) % MOD;
            }
        }
    }
    for (int i = 0; i < k; ++i) {
        a[i] = _c[i];
    }
}
int solve(int n, vector<int> a, vector<int> b) {
    int ans = 0, pnt = 0;
    int k = a.size();
    for (int i = 0; i < k; ++i) {
        _md[k - 1 - i] = -a[i];
    }
    _md[k] = 1;
    Md.clear();
    for (int i = 0; i < k; ++i) {
        if (_md[i]) {
            Md.push_back(i);
        }
    }
    for (int i = 0; i < k; ++i) {
        res[i] = base[i] = 0;
    }
    res[0] = 1;
    while ((1LL << pnt) <= n) {
        ++pnt;
    }
    for (int p = pnt; p >= 0; --p) {
        mul(res, res, k);
        if ((n >> p) & 1) {
            for (int i = k - 1; i >= 0; --i) {
                res[i + 1] = res[i];
            }
            res[0] = 0;
            for (int j = 0; j < Md.size(); ++j) {
                res[Md[j]] = (res[Md[j]] - res[k] * _md[Md[j]]) % MOD;
            }
        }
    }
    for (int i = 0; i < k; ++i) {
        ans = (ans + res[i] * b[i]) % MOD;
    }
    if (ans < 0) {
        ans += MOD;
    }
    return ans;
}

vector<int> calc(vector<int>& s) {
    vector<int> C(1, 1), B(1, 1);
    int L = 0, m = 1, b = 1;
    for (int n = 0; n < s.size(); ++n) {
        int d = 0;
        for (int i = 0; i <= L; ++i) {
            d = (d + C[i] * s[n - i]) % MOD;
        }
        if (d == 0) {
            ++m;
        } else if ((L << 1) <= n) {
            vector<int> T = C;
            int c = MOD - d * power(b, MOD - 2) % MOD;
            while (C.size() < B.size() + m) {
                C.push_back(0);
            }
            for (int i = 0; i < B.size(); ++i) {
                C[i + m] = (C[i + m] + c * B[i]) % MOD;
            }
            L = n + 1 - L; B = T; b = d; m = 1;
        } else {
            int c = MOD - d * power(b, MOD - 2) % MOD;
            while (C.size() < B.size() + m) {
                C.push_back(0);
            }
            for (int i = 0; i < B.size(); ++i) {
                C[i + m] = (C[i + m] + c * B[i]) % MOD;
            }
            ++m;
        }
    }
    return C;
}

int gao(vector<int>& a, int n) {
    vector<int> c = calc(a);
    c.erase(c.begin());
    for (int i = 0; i < c.size(); ++i) {
        c[i] = (MOD - c[i]) % MOD;
    }
    return solve(n, c, vector<int>(a.begin(), a.begin() + c.size()));
}

} // namespace BM
using BM::gao;