|
| 1 | +#pragma once |
| 2 | + |
| 3 | +template <class mint> |
| 4 | +struct FormalDirichletSeries : vector<mint> { |
| 5 | + using vector<mint>::vector; |
| 6 | + using FDS = FormalDirichletSeries; |
| 7 | + FDS& operator+=(const mint& v) { |
| 8 | + if (1 < this->size()) (*this)[1] += v; |
| 9 | + return *this; |
| 10 | + } |
| 11 | + FDS& operator+=(const FDS& r) { |
| 12 | + if (r.size() > this->size()) this->resize(r.size()); |
| 13 | + for (int i = 1; i < (int)r.size(); i++) (*this)[i] += r[i]; |
| 14 | + return *this; |
| 15 | + } |
| 16 | + FDS& operator-=(const mint& v) { |
| 17 | + if (1 < this->size()) (*this)[1] -= v; |
| 18 | + return *this; |
| 19 | + } |
| 20 | + FDS& operator-=(const FDS& r) { |
| 21 | + if (r.size() > this->size()) this->resize(r.size()); |
| 22 | + for (int i = 1; i < (int)r.size(); i++) (*this)[i] -= r[i]; |
| 23 | + return *this; |
| 24 | + } |
| 25 | + FDS& operator*=(const mint& v) { |
| 26 | + for (int i = 1; i < (int)this->size(); i++) (*this)[i] *= v; |
| 27 | + return *this; |
| 28 | + } |
| 29 | + FDS& operator*=(const FDS& r) { |
| 30 | + FDS p(this->size(), 0); |
| 31 | + for (int i = 1; i < (int)r.size(); i++) |
| 32 | + for (int j = 1; i * j < (int)p.size(); j++) |
| 33 | + p[i * j] += r[i] * (*this)[j]; |
| 34 | + return *this = p; |
| 35 | + } |
| 36 | + FDS& operator/=(const FDS& r) { |
| 37 | + mint c = r[1].inv(); |
| 38 | + for (int i = 1; i < (int)this->size(); i++) { |
| 39 | + (*this)[i] *= c; |
| 40 | + for (int j = 2; j < (int)r.size() && i * j < (int)this->size(); j++) |
| 41 | + (*this)[i * j] -= (*this)[i] * r[j]; |
| 42 | + } |
| 43 | + return *this; |
| 44 | + } |
| 45 | + FDS operator+(const mint& v) const { return FDS(*this) += v; } |
| 46 | + FDS operator+(const FDS& r) const { return FDS(*this) += r; } |
| 47 | + FDS operator-(const mint& v) const { return FDS(*this) -= v; } |
| 48 | + FDS operator-(const FDS& r) const { return FDS(*this) -= r; } |
| 49 | + FDS operator*(const mint& v) const { return FDS(*this) *= v; } |
| 50 | + FDS operator*(const FDS& r) const { return FDS(*this) *= r; } |
| 51 | + FDS operator/(const FDS& r) const { return FDS(*this) /= r; } |
| 52 | + FDS operator-() const { |
| 53 | + FDS ret(this->size()); |
| 54 | + for (int i = 1; i < (int)this->size(); i++) ret[i] = -(*this)[i]; |
| 55 | + return ret; |
| 56 | + } |
| 57 | + |
| 58 | + static vector<mint> inv; |
| 59 | + static void init_inv() { |
| 60 | + if (inv.empty()) { |
| 61 | + inv.assign(33, 0); |
| 62 | + inv[1] = mint(1); |
| 63 | + auto mod = mint::get_mod(); |
| 64 | + for (int i = 2; i < (int)inv.size(); i++) inv[i] = (-inv[mod % i]) * (mod / i); |
| 65 | + } |
| 66 | + } |
| 67 | + FDS integral() const { |
| 68 | + int n = this->size(); |
| 69 | + if (n <= 1) return FDS(n); |
| 70 | + calc_lpf(n); |
| 71 | + FDS ret(n); |
| 72 | + init_inv(); |
| 73 | + for (int i = 2; i < ret.size(); i++) ret[i] = (*this)[i] * inv[npf[i]]; |
| 74 | + return ret; |
| 75 | + } |
| 76 | + FDS diff() const { |
| 77 | + calc_lpf(this->size()); |
| 78 | + FDS ret(this->size()); |
| 79 | + for (int i = 1; i < (int)this->size(); i++) ret[i] = (*this)[i] * npf[i]; |
| 80 | + return ret; |
| 81 | + } |
| 82 | + FDS log(int n = -1) const { |
| 83 | + if (n == -1) n = (int)this->size(); |
| 84 | + if (n <= 1) return FDS(n); |
| 85 | + assert((*this)[1] == mint(1)); |
| 86 | + FDS f = *this; |
| 87 | + f.resize(n); |
| 88 | + return (f.diff() / f).integral(); |
| 89 | + } |
| 90 | + FDS exp(int n = -1) const { |
| 91 | + if (n == -1) n = (int)this->size(); |
| 92 | + if (n <= 1) return FDS(n); |
| 93 | + assert((*this)[1] == mint(0)); |
| 94 | + auto df = this->diff(); |
| 95 | + FDS ret(n); |
| 96 | + if (1 < n) ret[1] = 1; |
| 97 | + init_inv(); |
| 98 | + for (int i = 1; i < n; i++) { |
| 99 | + if (i > 1) ret[i] *= inv[npf[i]]; |
| 100 | + for (int j = 2; j < df.size() && i * j < n; j++) ret[i * j] += df[j] * ret[i]; |
| 101 | + } |
| 102 | + return ret; |
| 103 | + } |
| 104 | + |
| 105 | + static vector<int> npf; |
| 106 | + static vector<int> lpf; |
| 107 | + static void calc_lpf(int n) { |
| 108 | + if (lpf.empty()) { |
| 109 | + lpf = {0, 1, 2, 3, 2, 5, 2, 7}; |
| 110 | + npf = {0, 0, 1, 1, 2, 1, 2, 1}; |
| 111 | + } |
| 112 | + int l = lpf.size(), r = l; |
| 113 | + while (r <= n) r <<= 1; |
| 114 | + if (l == r) return; |
| 115 | + lpf.resize(r); |
| 116 | + npf.resize(r); |
| 117 | + for (int i = 2; i < l; i++) { |
| 118 | + if (lpf[i] != i) continue; |
| 119 | + for (int j = i * 2; j < r; j += i) |
| 120 | + if (lpf[j] == 0) lpf[j] = i; |
| 121 | + } |
| 122 | + for (int i = l; i < r; i++) { |
| 123 | + if (lpf[i] != 0) continue; |
| 124 | + lpf[i] = i; |
| 125 | + for (int j = i * 2; j < r; j += i) |
| 126 | + if (lpf[j] == 0) lpf[j] = i; |
| 127 | + } |
| 128 | + for (int i = l; i < r; i++) |
| 129 | + npf[i] = npf[i / lpf[i]] + 1; |
| 130 | + } |
| 131 | +}; |
| 132 | +template <typename mint> |
| 133 | +vector<mint> FormalDirichletSeries<mint>::inv = {}; |
| 134 | +template <typename mint> |
| 135 | +vector<int> FormalDirichletSeries<mint>::npf = {}; |
| 136 | +template <typename mint> |
| 137 | +vector<int> FormalDirichletSeries<mint>::lpf = {}; |
| 138 | +/** |
| 139 | + * @brief 形式的 Dirichlet 級数 |
| 140 | + * @docs docs/number-theory/formal-dirichlet-series.md |
| 141 | + */ |
0 commit comments