算法模版-数学

发表于 更新于 分类于 阅读次数： 本文字数： 11k
数学模版，慢慢更新。

数学

取模类（模int）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
template <class T>
constexpr T power(T a, i64 b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}

constexpr i64 mul(i64 a, i64 b, i64 p) {
    i64 res = a * b - i64(1.L * a * b / p) * p;
    res %= p;
    if (res < 0) {
        res += p;
    }
    return res;
}
template <int P>
struct MInt {
    int x;
    constexpr MInt() : x{} {}
    constexpr MInt(i64 x) : x{norm(x % getMod())} {}

    static int Mod;
    constexpr static int getMod() {
        if (P > 0) {
            return P;
        } else {
            return Mod;
        }
    }
    constexpr static void setMod(int Mod_) { Mod = Mod_; }
    constexpr int norm(int x) const {
        if (x < 0) {
            x += getMod();
        }
        if (x >= getMod()) {
            x -= getMod();
        }
        return x;
    }
    constexpr int val() const { return x; }
    explicit constexpr operator int() const { return x; }
    constexpr MInt operator-() const {
        MInt res;
        res.x = norm(getMod() - x);
        return res;
    }
    constexpr MInt inv() const {
        assert(x != 0);
        return power(*this, getMod() - 2);
    }
    constexpr MInt &operator*=(MInt rhs) & {
        x = 1LL * x * rhs.x % getMod();
        return *this;
    }
    constexpr MInt &operator+=(MInt rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    constexpr MInt &operator-=(MInt rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    constexpr MInt &operator/=(MInt rhs) & { return *this *= rhs.inv(); }
    friend constexpr MInt operator*(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res *= rhs;
        return res;
    }
    friend constexpr MInt operator+(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res += rhs;
        return res;
    }
    friend constexpr MInt operator-(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res -= rhs;
        return res;
    }
    friend constexpr MInt operator/(MInt lhs, MInt rhs) {
        MInt res = lhs;
        res /= rhs;
        return res;
    }
    friend constexpr std::istream &operator>>(std::istream &is, MInt &a) {
        i64 v;
        is >> v;
        a = MInt(v);
        return is;
    }
    friend constexpr std::ostream &operator<<(std::ostream &os, const MInt &a) {
        return os << a.val();
    }
    friend constexpr bool operator==(MInt lhs, MInt rhs) {
        return lhs.val() == rhs.val();
    }
    friend constexpr bool operator!=(MInt lhs, MInt rhs) {
        return lhs.val() != rhs.val();
    }
};
template <>
int MInt<0>::Mod = 998244353;

template <int V, int P>
constexpr MInt<P> CInv = MInt<P>(V).inv();

constexpr int P = 1000000007;
using Z = MInt<P>;  // 静态模数，题目给定模数
using D = MInt<0>;  // 动态模数，可以 D::setMod(x) 改变

线性筛（欧拉筛）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// minp(x) 可以找到 x 大于 1 的最小质因子
std::vector<int> minp, primes;
void sieve(int n) {
    minp.assign(n + 1, 0);
    primes.clear();
    for (int i = 2; i <= n; i++) {
        if (minp[i] == 0) {
            minp[i] = i;
            primes.push_back(i);
        }
        for (auto p : primes) {
            if (i * p > n) {
                break;
            }
            minp[i * p] = p;
            if (p == minp[i]) {
                break;
            }
        }
    }
}

组合数类（搭配取模类）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
struct Comb {
    int n;
    std::vector<Z> _fac;
    std::vector<Z> _invfac;
    std::vector<Z> _inv;
    
    Comb() : n{0}, _fac{1}, _invfac{1}, _inv{0} {}
    Comb(int n) : Comb() {
        init(n);
    }
    
    void init(int m) {
        if (m <= n) return;
        _fac.resize(m + 1);
        _invfac.resize(m + 1);
        _inv.resize(m + 1);
        
        for (int i = n + 1; i <= m; i++) {
            _fac[i] = _fac[i - 1] * i;
        }
        _invfac[m] = _fac[m].inv();
        for (int i = m; i > n; i--) {
            _invfac[i - 1] = _invfac[i] * i;
            _inv[i] = _invfac[i] * _fac[i - 1];
        }
        n = m;
    }
    
    Z fac(int m) {
        if (m > n) init(2 * m);
        return _fac[m];
    }
    Z invfac(int m) {
        if (m > n) init(2 * m);
        return _invfac[m];
    }
    Z inv(int m) {
        if (m > n) init(2 * m);
        return _inv[m];
    }
    Z binom(int n, int m) {
        if (n < m || m < 0) return 0;
        return fac(n) * invfac(m) * invfac(n - m);
    }
} comb(200000); // Z 类型全局初始化，D 类型需要局部初始化

组合数（小范围预处理，逆元+杨辉三角）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
constexpr int P = 1000000007;
constexpr int L = 10000;

int fac[L + 1], invfac[L + 1];
int sumbinom[L + 1][7];

int binom(int n, int m) {
    if (n < m || m < 0) {
        return 0;
    }
    return 1LL * fac[n] * invfac[m] % P * invfac[n - m] % P;
}

int power(int a, int b) {
    int res = 1;
    for (; b; b /= 2, a = 1LL * a * a % P) {
        if (b % 2) {
            res = 1LL * res * a % P;
        }
    }
    return res;
}

int main() {
    fac[0] = 1;
    for (int i = 1; i <= L; i++) {
        fac[i] = 1LL * fac[i - 1] * i % P;
    }
    invfac[L] = power(fac[L], P - 2);
    for (int i = L; i; i--) {
        invfac[i - 1] = 1LL * invfac[i] * i % P;
    }

    sumbinom[0][0] = 1;
    for (int i = 1; i <= L; i++) {
        for (int j = 0; j < 7; j++) {
            sumbinom[i][j] = (sumbinom[i - 1][j] + sumbinom[i - 1][(j + 6) % 7]) % P;
        }
    }
}

矩阵（以斐波那契为例）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
constexpr int mod = 1e9 + 7;
struct Matrix {
    vector<vector<i64>> a;
    Matrix() { a.assign(3, vector<i64>(3)); }
    Matrix operator*(const Matrix b) const {
        Matrix res;
        for (int i = 1; i <= 2; i++) {
            for (int j = 1; j <= 2; j++) {
                for (int k = 1; k <= 2; k++) {
                    res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
                }
            }
        }
        return res;
    }
} ans, base;

void qpow(i64 b) {
    while (b) {
        if (b % 2 == 1) {
            ans = ans * base;
        }
        base = base * base;
        b /= 2;
    }
}

void init() {
    base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
    ans.a[1][1] = ans.a[1][2] = 1;
}

void solve() {
    i64 n;
    cin >> n;
    if (n <= 2) {
        cout << 1 << "\n";
        return;
    }
    init();
    qpow(n - 2);
    cout << ans.a[1][1] % mod << "\n";
}

扩展欧几里得 exgcd，线性同余方程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
// 返回 gcd(a, b) ax + by = gcd(a, b)
i64 exgcd(i64 a, i64 b, i64 &x, i64 &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}

// ax + b = 0 (mod m) 返回最小正整数解，通解的乘数
std::pair<i64, i64> sol(i64 a, i64 b, i64 m) {
    assert(m > 0);
    b *= -1;
    i64 x, y;
    i64 g = exgcd(a, m, x, y);
    if (g < 0) {
        g *= -1;
        x *= -1;
        y *= -1;
    }
    if (b % g != 0) {
        return {-1, -1};
    }
    x = x * (b / g) % (m / g);
    if (x < 0) {
        x += m / g;
    }
    return {x, m / g};
}

std::array<i64, 3> exgcd(i64 a, i64 b) {
    if (!b) {
        return {a, 1, 0};
    }
    auto [g, x, y] = exgcd(b, a % b);
    return {g, y, x - a / b * y};
}

线性求逆元

1
2
3
4
5
6
7
8
std::vector<int> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = 1ll * (p - p / i) * inv[p % i] % p;
    }
    for (int i = 1; i <= n; i++) {
        std::cout << inv[i] << "\n";
    }

中国剩余定理 CRT

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
// 返回 gcd(a, b) ax + by = gcd(a, b)
__int128 exgcd(__int128 a, __int128 b, __int128 &x, __int128 &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    auto g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}

__int128 crt(int n, std::vector<i64> &a, std::vector<i64> &b) {
    __int128 M = 1, ans = 0;
    for (int i = 1; i <= n; i++) {
        M = M * b[i];
    }
    for (int i = 1; i <= n; i++) {
        __int128 m = M / b[i];
        __int128 x, y;
        exgcd(m, b[i], x, y);
        ans = (ans + a[i] * m * x) % M;
    }
    ans = (ans % M + M) % M;
    return ans;
};

分解质因数

记住循环最后还有一个 \(n\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int x;
  std::cin >> x;

  for (int i = 2; i * i <= x; i++) {
      if (x % i == 0) {
          while (x % i == 0) {
              std::cout << i << "*";
              x /= i;
          }
      }
  }
  if (x != 1) {
      std::cout << x << "\n";
  }

求所有因数

1
2
3
4
5
6
7
8
9
10
11
12
13
std::vector<int> get_divisors(int n) {
    std::vector<int> res;
    for (int i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            res.push_back(i);
            if (i != n / i) {
                res.push_back(n / i);
            }
        }
    }
    std::sort(res.begin(), res.end());
    return res;
}

欧拉函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
// 求解单个数的欧拉函数
int phi(int n) {
  	int res = n;
  	for (int i = 2; i * i <= n; i++) {
    		if(n % i == 0) {
     			 while(n % i == 0) {
        			n /= i;
     			 }
     		res = res / i * (i - 1);
    }
  }
  if(n > 1) {
    res = res / n * (n - 1);
  } 
  return res;
}

// 求解全部数的欧拉函数
constexpr int N = 1e7;
constexpr int P = 1000003;

bool isprime[N + 1];
int phi[N + 1];
std::vector<int> primes;

    std::fill(isprime + 2, isprime + N + 1, true);
    phi[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (isprime[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (auto p : primes) {
            if (i * p > N) {
                break;
            }
            isprime[i * p] = false;
            if (i % p == 0) {
                phi[i * p] = phi[i] * p;
                break;
            }
            phi[i * p] = phi[i] * (p - 1);
        }
    }

欧拉定理，扩展欧拉定理

欧拉定理：若 \(\gcd(a,m) = 1, a^{\varphi(m)} = 1\pmod m,a^b = a^{b\pmod{\varphi(m)}}\pmod m\)

扩展欧拉定理：\(a^b = a^{b\pmod {\varphi(m)} + \varphi(m)}\pmod m, b\ge\varphi(m)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <bits/stdc++.h>

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;
using u128 = unsigned __int128;

// 求解单个数的欧拉函数
int phi(int n) {
    int res = n;
    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            while (n % i == 0) {
                n /= i;
            }
            res = res / i * (i - 1);
        }
    }
    if (n > 1) {
        res = res / n * (n - 1);
    }
    return res;
}

i64 qpow(i64 a, i64 b, i64 mod) {
    i64 res = 1;
    while (b) {
        if (b & 1) {
            res = res * a % mod;
        }
        a = a * a % mod;
        b /= 2;
    }
    return res;
}

int main() {
    // std::ios::sync_with_stdio(false);
    // std::cin.tie(nullptr);

    int a, m, b;
    std::cin >> a >> m;
    a %= m;
    int phim = phi(m);
    char ch;
    while ((ch = getchar()) < '0' || ch > '9');
    i64 bm = 0;
    bool flag = false;
    while (bm = bm * 10 + (ch ^ '0'), (ch = getchar()) >= '0' && ch <= '9') {
        if (bm >= phim) {
            flag = true;
            bm %= phim;
        }
    }
    if (bm >= phim) {
        flag = true;
        bm %= phim;
    }
    if (flag) {
        bm += phim;
    }
    std::cout << qpow(a, bm, m) << "\n";

    return 0;
}

第二类斯特林数（预处理阶乘及逆元）

将 \(n\) 个两两不同的元素，划分为 \(k\) 个互不区分的非空子集的方案数 \(\left\{ \begin{matrix} n\\k \end{matrix} \right\}\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
std::vector<Z> frac(n + 1), invfrac(n + 1);
  frac[0] = frac[1] = invfrac[0] = invfrac[1] = Z(1);
  for (int i = 2; i <= n; i++) {
      frac[i] = frac[i - 1] * i;
  }
  invfrac[n] = power(frac[n], P - 2);
  for (int i = n - 1; i >= 1; i--) {
      invfrac[i] = invfrac[i + 1] * (i + 1);
  }
  auto S = [&](int n, int k) {
      Z ans = 0;
      for (int i = 0; i <= k; i++) {
          ans += power(Z(-1), k - i) * power(Z(i), n) * invfrac[i] *
                 invfrac[k - i];
      }
      return ans;
  };

数论分块

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 例：求 \sum_{i=1}^n k \pmod i = nk - \sum_{i=1}^n i(k/i)
i64 ans = 1ll * n * k;
  for (int l = 1, r; l <= n; l = r + 1) {
      int q = k / l;
      if (q == 0) { 
          r = n;
      } else {
          r = std::min(k / q, n);
      }
      ans -= 1ll * (l + r) * (r - l + 1) / 2 * q;
  }
	

// 向上取整的数论分块
for (int l = 1, r; l <= n; l = r + 1) {
      int k = (m + l - 1) / l;
      r = (k == 1) ? n : (m - 1) / (k - 1);
      // 区间 [l, r] \in [1, n] 内，所有 ceil(m / i) 的值都等于 k
      ans = std::min(ans, (k - 1) * l + n - m);
  }

莫比乌斯函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
// 注意莫比乌斯函数是前缀和的形式
constexpr int N = 1e7;
std::unordered_map<int, Z> fMu;
std::vector<int> minp, primes, phi, mu;
std::vector<i64> sphi;

void sieve(int n) {
    minp.assign(n + 1, 0);
    phi.assign(n + 1, 0);
    sphi.assign(n + 1, 0);
    mu.assign(n + 1, 0);
    primes.clear();
    phi[1] = 1;
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (minp[i] == 0) {
            minp[i] = i;
            phi[i] = i - 1;
            mu[i] = -1;
            primes.push_back(i);
        }
        for (auto p : primes) {
            if (i * p > n) {
                break;
            }
            minp[i * p] = p;
            if (p == minp[i]) {
                phi[i * p] = phi[i] * p;
                break;
            }
            phi[i * p] = phi[i] * (p - 1);
            mu[i * p] = -mu[i];
        }
    }
    for (int i = 1; i <= n; i++) {
        sphi[i] = sphi[i - 1] + phi[i];
        mu[i] += mu[i - 1];
    }
}

// 快速求 \sum_{i=1}^n \mu(i) 的杜教筛
Z sumMu(int n) {
    if (n <= N) {
        return mu[n];
    }
    if (fMu.count(n)) {
        return fMu[n];
    }
    if (n == 0) {
        return 0;
    }
    Z ans = 1;
    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans -= (r - l + 1) * sumMu(n / l);
    }
    fMu[n] = ans;
    return ans;
}