diff --git a/.verify-helper/timestamps.remote.json b/.verify-helper/timestamps.remote.json index 411a4d7be..fb1f74f83 100644 --- a/.verify-helper/timestamps.remote.json +++ b/.verify-helper/timestamps.remote.json @@ -62,7 +62,7 @@ "verify/verify-aoj-other/aoj-1130-DG-bfs.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-aoj-other/aoj-1377.test.cpp": "2023-08-27 10:14:06 +0900", "verify/verify-aoj-other/aoj-1613.test.cpp": "2023-08-10 13:25:59 +0900", -"verify/verify-aoj-other/aoj-2171-bigrational.test.cpp": "2023-08-10 14:06:55 +0900", +"verify/verify-aoj-other/aoj-2171-bigrational.test.cpp": "2023-09-02 23:13:20 +0900", "verify/verify-aoj-other/aoj-2171.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-aoj-other/aoj-2891.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-aoj-other/aoj-2945-01bfs.test.cpp": "2023-08-10 13:25:59 +0900", @@ -80,6 +80,7 @@ "verify/verify-unit-test/bigint.test.cpp": "2023-08-10 14:06:55 +0900", "verify/verify-unit-test/bigint2.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-unit-test/bigint3.test.cpp": "2023-08-10 14:25:31 +0900", +"verify/verify-unit-test/bigrational.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-unit-test/bitset-find-prev.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-unit-test/complex-fft.test.cpp": "2023-08-10 14:06:55 +0900", "verify/verify-unit-test/composite-exp.test.cpp": "2023-08-31 20:44:07 +0900", @@ -127,7 +128,7 @@ "verify/verify-unit-test/primitive-root.test.cpp": "2023-08-31 20:44:07 +0900", "verify/verify-unit-test/radix-heap.test.cpp": "2023-08-10 14:06:55 +0900", "verify/verify-unit-test/radix-sort.test.cpp": "2023-08-10 14:25:31 +0900", -"verify/verify-unit-test/rational-number.test.cpp": "2023-08-10 14:06:55 +0900", +"verify/verify-unit-test/rational-number.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-unit-test/rbst-segment-tree.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-unit-test/rbst-sequence.test.cpp": "2023-08-10 14:06:55 +0900", "verify/verify-unit-test/relaxed-convolution.test.cpp": "2023-08-31 20:44:07 +0900", @@ -149,6 +150,7 @@ "verify/verify-yosupo-ds/yosupo-associative-array-rbstseg.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-ds/yosupo-associative-array-unerasable-hashmap.test.cpp": "2023-08-27 10:14:06 +0900", "verify/verify-yosupo-ds/yosupo-binary-trie.test.cpp": "2023-08-10 13:25:59 +0900", +"verify/verify-yosupo-ds/yosupo-deque-operate-all-composite.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yosupo-ds/yosupo-dynamic-li-chao-tree.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-ds/yosupo-dynamic-sequence-range-affine-range-sum-splay.test.cpp": "2023-08-10 18:41:15 +0900", "verify/verify-yosupo-ds/yosupo-dynamic-sequence-range-affine-range-sum-treap.test.cpp": "2023-08-10 14:25:31 +0900", @@ -183,6 +185,8 @@ "verify/verify-yosupo-ds/yosupo-predecessor-problem.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-yosupo-ds/yosupo-procedessor-problem-rbstseg.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-ds/yosupo-range-add-range-sum-linkcuttree.test.cpp": "2023-08-10 18:41:15 +0900", +"verify/verify-yosupo-ds/yosupo-range-affine-point-get-2.test.cpp": "2023-09-02 22:21:41 +0900", +"verify/verify-yosupo-ds/yosupo-range-affine-point-get.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yosupo-ds/yosupo-range-affine-range-sum-dynamic-segtree.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-ds/yosupo-range-affine-range-sum-rbstseg.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-ds/yosupo-range-affine-sqdec.test.cpp": "2023-08-10 13:25:59 +0900", @@ -194,7 +198,7 @@ "verify/verify-yosupo-ds/yosupo-static-range-inversion-query-2.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-yosupo-ds/yosupo-static-range-inversions-query.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-yosupo-ds/yosupo-static-rmq.test.cpp": "2023-08-10 13:25:59 +0900", -"verify/verify-yosupo-ds/yosupo-swag.test.cpp": "2023-08-10 14:25:31 +0900", +"verify/verify-yosupo-ds/yosupo-swag.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yosupo-ds/yosupo-vertex-add-path-sum-euler-tour.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-yosupo-ds/yosupo-vertex-add-path-sum.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-ds/yosupo-vertex-add-subtree-sum-dst-on-tree.test.cpp": "2023-08-10 14:25:31 +0900", @@ -267,6 +271,7 @@ "verify/verify-yosupo-math/yosupo-counting-primes-3.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-math/yosupo-counting-primes-4.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-math/yosupo-counting-primes.test.cpp": "2023-08-10 13:25:59 +0900", +"verify/verify-yosupo-math/yosupo-determinant-arbitrary-mod.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yosupo-math/yosupo-determinant-matrixlib.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-math/yosupo-determinant-of-matrix-bbla.test.cpp": "2023-08-31 20:44:07 +0900", "verify/verify-yosupo-math/yosupo-determinant-of-sparse-matrix-bbla.test.cpp": "2023-08-31 20:44:07 +0900", @@ -294,6 +299,7 @@ "verify/verify-yosupo-math/yosupo-primality-test.test.cpp": "2023-08-10 14:25:31 +0900", "verify/verify-yosupo-math/yosupo-prime-enumerate-sieve.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-math/yosupo-prime-table.test.cpp": "2023-08-10 14:25:31 +0900", +"verify/verify-yosupo-math/yosupo-primitive-root.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yosupo-math/yosupo-sparse-determinant.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-math/yosupo-stern-brocot-tree.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yosupo-math/yosupo-subset-convolution-fast.test.cpp": "2023-08-10 14:25:31 +0900", @@ -361,7 +367,7 @@ "verify/verify-yuki/yuki-0886.test.cpp": "2023-08-30 23:05:07 +0900", "verify/verify-yuki/yuki-0890.test.cpp": "2023-08-30 23:05:07 +0900", "verify/verify-yuki/yuki-0896.test.cpp": "2023-08-30 23:05:07 +0900", -"verify/verify-yuki/yuki-0952.test.cpp": "2023-08-10 14:06:55 +0900", +"verify/verify-yuki/yuki-0952.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yuki/yuki-0963-circular.test.cpp": "2023-08-31 20:44:07 +0900", "verify/verify-yuki/yuki-0963.test.cpp": "2023-08-31 20:44:07 +0900", "verify/verify-yuki/yuki-1080.test.cpp": "2023-08-31 20:44:07 +0900", @@ -396,11 +402,12 @@ "verify/verify-yuki/yuki-1786.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yuki/yuki-1789.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yuki/yuki-1875.test.cpp": "2023-08-31 20:44:07 +0900", +"verify/verify-yuki/yuki-1939-sparse-pow.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yuki/yuki-1939.test.cpp": "2023-08-31 20:44:07 +0900", "verify/verify-yuki/yuki-2012.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yuki/yuki-2231.test.cpp": "2023-08-10 13:25:59 +0900", -"verify/verify-yuki/yuki-2262.test.cpp": "2023-08-12 16:25:29 +0900", -"verify/verify-yuki/yuki-2266.test.cpp": "2023-08-12 16:25:29 +0900", +"verify/verify-yuki/yuki-2262.test.cpp": "2023-09-02 22:21:41 +0900", +"verify/verify-yuki/yuki-2266.test.cpp": "2023-09-02 22:21:41 +0900", "verify/verify-yuki/yuki-2281.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yuki/yuki-2333.test.cpp": "2023-08-10 13:25:59 +0900", "verify/verify-yuki/yuki-helloworld.test.cpp": "2023-08-10 13:25:59 +0900" diff --git a/data-structure/slide-window-aggregation-deque.hpp b/data-structure/slide-window-aggregation-deque.hpp new file mode 100644 index 000000000..341fcd9a4 --- /dev/null +++ b/data-structure/slide-window-aggregation-deque.hpp @@ -0,0 +1,58 @@ +#pragma once + +#include +using namespace std; + +template +struct SlideWindowAggregationDeque { + vector a0, a1, r0, r1; + F f; + T I; + + SlideWindowAggregationDeque(F _f, T _i) : f(_f), I(_i) {} + + private: + T get0() const { return r0.empty() ? I : r0.back(); } + T get1() const { return r1.empty() ? I : r1.back(); } + + void push0(const T &x) { + a0.push_back(x); + r0.push_back(f(x, get0())); + } + void push1(const T &x) { + a1.push_back(x); + r1.push_back(f(get1(), x)); + } + void rebalance() { + int n = a0.size() + a1.size(); + int s0 = n / 2 + (a0.empty() ? n % 2 : 0); + vector a{a0}; + reverse(begin(a), end(a)); + copy(begin(a1), end(a1), back_inserter(a)); + a0.clear(), r0.clear(); + a1.clear(), r1.clear(); + for (int i = s0 - 1; i >= 0; i--) push0(a[i]); + for (int i = s0; i < n; i++) push1(a[i]); + } + + public: + void push_front(const T &t) { push0(t); } + void push_back(const T &t) { push1(t); } + T front() const { return a0.empty() ? a1.front() : a0.back(); } + T back() const { return a1.empty() ? a0.front() : a1.back(); } + void pop_front() { + if (a0.empty()) rebalance(); + assert(!a0.empty()); + a0.pop_back(), r0.pop_back(); + } + void pop_back() { + if (a1.empty()) rebalance(); + assert(!a1.empty()); + a1.pop_back(), r1.pop_back(); + } + T query() { return f(get0(), get1()); } +}; + +/** + * @brief Slide Window Aggrigation (deque) + */ diff --git a/data-structure/slide-window-aggregation.hpp b/data-structure/slide-window-aggregation.hpp index 78442b4e6..527b143ac 100644 --- a/data-structure/slide-window-aggregation.hpp +++ b/data-structure/slide-window-aggregation.hpp @@ -1,29 +1,31 @@ #pragma once +#include +using namespace std; + template struct SlideWindowAggregation { - stack a0, a1; - stack r0, r1; + vector a0, a1, r0, r1; F f; T I, f0, f1; - SlideWindowAggregation(F f_, T I_) : f(f_), I(I_), f0(I_), f1(I_) {} + SlideWindowAggregation(F _f, T _i) : f(_f), I(_i), f0(_i), f1(_i) {} private: void push_s0(const T &x) { - a0.push(x); - r0.push(f0 = f(x, f0)); + a0.push_back(x); + r0.push_back(f0 = f(x, f0)); } void push_s1(const T &x) { - a1.push(x); - r1.push(f1 = f(f1, x)); + a1.push_back(x); + r1.push_back(f1 = f(f1, x)); } void transfer() { while (!a1.empty()) { - push_s0(a1.top()); - a1.pop(); + push_s0(a1.back()); + a1.pop_back(); } - while (!r1.empty()) r1.pop(); + while (!r1.empty()) r1.pop_back(); f1 = I; } @@ -38,9 +40,9 @@ struct SlideWindowAggregation { } void pop() { if (a0.empty()) transfer(); - a0.pop(); - r0.pop(); - f0 = r0.empty() ? I : r0.top(); + a0.pop_back(); + r0.pop_back(); + f0 = r0.empty() ? I : r0.back(); } T query() { return f(f0, f1); } }; diff --git a/dp/monge-d-edge-shortest-path-enumerate.hpp b/dp/monge-d-edge-shortest-path-enumerate.hpp index 5b2304d64..dc11f9ee1 100644 --- a/dp/monge-d-edge-shortest-path-enumerate.hpp +++ b/dp/monge-d-edge-shortest-path-enumerate.hpp @@ -16,7 +16,7 @@ vector enumerate_monge_d_edge_shortest_path( vector dp(N + 1, INF); dp[0] = 0; for (int d = 1; d <= N; d++) { - vector midx = monotone_minima(N + 1, [&](int j, int i) -> T { + vector midx = monotone_minima(N + 1, N + 1, [&](int j, int i) -> T { return i < j ? dp[i] + f(i, j) : INF; }); for (int i = N; i >= d; i--) dp[i] = dp[midx[i]] + f(midx[i], i); diff --git a/dp/monotone-minima.hpp b/dp/monotone-minima.hpp index 90d54d086..33d3ad993 100644 --- a/dp/monotone-minima.hpp +++ b/dp/monotone-minima.hpp @@ -6,27 +6,38 @@ using namespace std; // NxN 行列がある // m_i := argmin_j (A_{i,j}) が単調増加であるときに m_i を列挙する -// A(i, j) : A[i][j] を返す関数 -template -vector monotone_minima(int N, const function& A) { +// f(i, j, k) : +// A[i][j] と A[i][k] を比較 (j < k が保証されている) +// A[i][j] <= A[i][k] のとき true を返す関数を入れる (等号はどちらでもよい) +vector monotone_minima(int N, int M, + const function& f) { vector res(N); auto dfs = [&](auto rc, int is, int ie, int l, int r) -> void { if (is == ie) return; int i = (is + ie) / 2; int m = l; - T fim = A(i, m); for (int k = l + 1; k < r; k++) { - T fik = A(i, k); - if (fik < fim) fim = fik, m = k; + if (!f(i, m, k)) m = k; } res[i] = m; rc(rc, is, i, l, m + 1); rc(rc, i + 1, ie, m, r); }; - dfs(dfs, 0, N, 0, N); + dfs(dfs, 0, N, 0, M); return res; } +// NxM 行列がある +// m_i := argmin_j (A_{i,j}) が単調増加であるときに m_i を列挙する +// A(i, j) : A[i][j] を返す関数 +template +vector monotone_minima(int N, int M, const function& A) { + function f = [&](int i, int j, int k) -> bool { + return A(i, j) <= A(i, k); + }; + return monotone_minima(N, M, f); +} + /** * @brief monotone minima */ diff --git a/marathon/multi-armed-bandit.hpp b/marathon/multi-armed-bandit.hpp new file mode 100644 index 000000000..9a950dfc6 --- /dev/null +++ b/marathon/multi-armed-bandit.hpp @@ -0,0 +1,55 @@ +#include +#include +#include +#include +#include +using namespace std; + +#include "../misc/rng.hpp" + +// N 択, 報酬最大化 +struct MultiArmedBandit { + MultiArmedBandit(int n) + : N(n), last(-1), iter(0), thres(N * 5), num(N), v(N), e(N), t(1) {} + + int N, last; + long long iter, thres; + vector num; + vector v, e; + double t; + + int play() { + assert(last == -1); + iter++; + if (iter <= thres) return last = iter % N; + + double s = accumulate(begin(e), end(e), 0.0); + double x = rnd() * s; + for (int i = 0; i < N; i++) { + if ((x -= e[i]) <= 0) return last = i; + } + return last = N - 1; + } + + // 重み付け用の関数 + double f(double x) { return exp(x / t); } + + void reward(double y) { + assert(last != -1); + v[last] += y; + num[last] += 1; + e[last] = f(v[last] / num[last]); + last = -1; + + static double u = 1.0; + static double du = 0.01; + // iter % thres == 0 になったら t を再決定 + if (iter % thres == 0) { + u = max(0.7, u - du); + double average = accumulate(begin(v), end(v), 0.0) / thres; + t = average < 0.0 ? 1.0 : pow(average, u); + for (int i = 0; i < N; i++) e[i] = f(v[i] / num[i]); + } + } + int best() { return max_element(begin(e), end(e)) - begin(e); } +}; diff --git a/math/bigint-all.hpp b/math/bigint-all.hpp new file mode 100644 index 000000000..827154ca7 --- /dev/null +++ b/math/bigint-all.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include "rational-binomial.hpp" +#include "rational-fps.hpp" +#include "rational.hpp" +// +#include "bigint-rational.hpp" +#include "bigint.hpp" +// +using mint = BigRational; +using vm = vector; +using vvm = vector; +using fps = FormalPowerSeries_rational; +Binomial_rational C; diff --git a/math/primitive-root-ll.hpp b/math/primitive-root-ll.hpp new file mode 100644 index 000000000..84ceb4e46 --- /dev/null +++ b/math/primitive-root-ll.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include +#include +using namespace std; + +#include "../internal/internal-math.hpp" +#include "../prime/fast-factorize.hpp" + +long long primitive_root_ll(long long p) { + if (p == 2) return 1; + auto fs = factorize(p - 1); + sort(begin(fs), end(fs)); + fs.erase(unique(begin(fs), end(fs)), end(fs)); + for (int g = 2;; g++) { + int ok = 1; + for (auto& f : fs) { + if (internal::modpow(g, (p - 1) / f, p) == 1) { + ok = false; + break; + } + } + if (ok) return g; + } + exit(1); +} diff --git a/math/rational-binomial.hpp b/math/rational-binomial.hpp new file mode 100644 index 000000000..620c7552b --- /dev/null +++ b/math/rational-binomial.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "rational.hpp" + +template +struct Binomial_rational { + vector fc; + Binomial_rational(int = 0) { fc.emplace_back(1); } + void extend() { + int n = fc.size(); + R nxt = fc.back() * n; + fc.push_back(nxt); + } + R fac(int n) { + if (n < 0) return 0; + while ((int)fc.size() <= n) extend(); + return fc[n]; + } + R finv(int n) { + if (n < 0) return 0; + return fac(n).inverse(); + } + R inv(int n) { + if (n < 0) return -inv(-n); + return R{1, max(n, 1)}; + } + R C(int n, int r) { + if (n < 0 or r < 0 or n < r) return R{0}; + return fac(n) * finv(n - r) * finv(r); + } + R operator()(int n, int r) { return C(n, r); } + template + R multinomial(const vector& r) { + static_assert(is_integral::value == true); + int n = 0; + for (auto& x : r) { + if (x < 0) return R{0}; + n += x; + } + R res = fac(n); + for (auto& x : r) res *= finv(x); + return res; + } + + template + R operator()(const vector& r) { + return multinomial(r); + } +}; diff --git a/math/rational-fps.hpp b/math/rational-fps.hpp new file mode 100644 index 000000000..1f0d400dc --- /dev/null +++ b/math/rational-fps.hpp @@ -0,0 +1,197 @@ +#pragma once + +#include "rational-binomial.hpp" +#include "rational.hpp" + +template +struct FormalPowerSeries_rational : vector { + using vector::vector; + using fps = FormalPowerSeries_rational; + + fps &operator+=(const fps &r) { + if (r.size() > this->size()) this->resize(r.size()); + for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i]; + return *this; + } + + fps &operator+=(const R &r) { + if (this->empty()) this->resize(1); + (*this)[0] += r; + return *this; + } + + fps &operator-=(const fps &r) { + if (r.size() > this->size()) this->resize(r.size()); + for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i]; + return *this; + } + + fps &operator-=(const R &r) { + if (this->empty()) this->resize(1); + (*this)[0] -= r; + return *this; + } + + fps &operator*=(const fps &r) { + int n = this->size() + r.size() - 1; + fps f(n); + for (int i = 0; i < (int)this->size(); i++) { + for (int j = 0; j < (int)r.size(); j++) { + f[i + j] += (*this)[i] * r[j]; + } + } + return *this = f; + } + + fps &operator*=(const R &v) { + for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v; + return *this; + } + + fps &operator/=(const fps &r) { + if (this->size() < r.size()) { + this->clear(); + return *this; + } + int n = this->size() - r.size() + 1; + fps f(*this), g(r); + g.shrink(); + R coeff = g.back().inverse(); + for (auto &x : g) x *= coeff; + int deg = (int)f.size() - (int)g.size() + 1; + int gs = g.size(); + fps quo(deg); + for (int i = deg - 1; i >= 0; i--) { + quo[i] = f[i + gs - 1]; + for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j]; + } + *this = quo * coeff; + this->resize(n, R(0)); + return *this; + } + + fps &operator%=(const fps &r) { + *this -= *this / r * r; + shrink(); + return *this; + } + + fps operator+(const fps &r) const { return fps(*this) += r; } + fps operator+(const R &v) const { return fps(*this) += v; } + fps operator-(const fps &r) const { return fps(*this) -= r; } + fps operator-(const R &v) const { return fps(*this) -= v; } + fps operator*(const fps &r) const { return fps(*this) *= r; } + fps operator*(const R &v) const { return fps(*this) *= v; } + fps operator/(const fps &r) const { return fps(*this) /= r; } + fps operator%(const fps &r) const { return fps(*this) %= r; } + fps operator-() const { + fps ret(this->size()); + for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i]; + return ret; + } + + void shrink() { + while (this->size() && this->back() == R(0)) this->pop_back(); + } + + fps rev() const { + fps ret(*this); + reverse(begin(ret), end(ret)); + return ret; + } + + fps dot(fps r) const { + fps ret(min(this->size(), r.size())); + for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i]; + return ret; + } + + // 前 sz 項を取ってくる。sz に足りない項は 0 埋めする + fps pre(int sz) const { + fps ret(begin(*this), begin(*this) + min((int)this->size(), sz)); + if ((int)ret.size() < sz) ret.resize(sz); + return ret; + } + + fps operator>>(int sz) const { + if ((int)this->size() <= sz) return {}; + fps ret(*this); + ret.erase(ret.begin(), ret.begin() + sz); + return ret; + } + + fps operator<<(int sz) const { + fps ret(*this); + ret.insert(ret.begin(), sz, R(0)); + return ret; + } + + fps diff() const { + const int n = (int)this->size(); + fps ret(max(0, n - 1)); + R one(1), coeff(1); + for (int i = 1; i < n; i++) { + ret[i - 1] = (*this)[i] * coeff; + coeff += one; + } + return ret; + } + + fps integral() const { + const int n = (int)this->size(); + fps ret(n + 1); + for (int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / (i + 1); + return ret; + } + + R eval(R x) const { + R r = 0, w = 1; + for (auto &v : *this) r += w * v, w *= x; + return r; + } + + fps inv(int deg = -1) const { + assert((*this)[0] != R(0)); + if (deg == -1) deg = (*this).size(); + fps ret{R(1) / (*this)[0]}; + for (int i = 1; i < deg; i <<= 1) { + ret = (ret + ret - ret * ret * (*this).pre(i << 1)).pre(i << 1); + } + return ret.pre(deg); + } + fps log(int deg = -1) const { + assert(!(*this).empty() && (*this)[0] == R(1)); + if (deg == -1) deg = (int)this->size(); + return (this->diff() * this->inv(deg)).pre(deg - 1).integral(); + } + fps exp(int deg = -1) const { + assert((*this).size() == 0 || (*this)[0] == R(0)); + if (deg == -1) deg = (int)this->size(); + fps ret{R(1)}; + for (int i = 1; i < deg; i <<= 1) { + ret = (ret * (pre(i << 1) + R(1) - ret.log(i << 1))).pre(i << 1); + } + return ret.pre(deg); + } + fps pow(int64_t k, int deg = -1) const { + const int n = (int)this->size(); + if (deg == -1) deg = n; + if (k == 0) { + fps ret(deg); + if (deg) ret[0] = 1; + return ret; + } + for (int i = 0; i < n; i++) { + if ((*this)[i] != R(0)) { + R rev = R(1) / (*this)[i]; + fps ret = (((*this * rev) >> i).log(deg) * k).exp(deg); + ret *= (*this)[i].pow(k); + ret = (ret << (i * k)).pre(deg); + if ((int)ret.size() < deg) ret.resize(deg, R(0)); + return ret; + } + if (__int128_t(i + 1) * k >= deg) return fps(deg, R(0)); + } + return fps(deg, R(0)); + } +}; diff --git a/math/rational.hpp b/math/rational.hpp index 454fb9aa7..37b7c855c 100644 --- a/math/rational.hpp +++ b/math/rational.hpp @@ -92,6 +92,7 @@ struct RationalBase { return os; } + // T にキャストされるので T が bigint の場合は to_ll も要る T to_mint(T mod) const { assert(mod != 0); T a = y, b = mod, u = 1, v = 0, t; @@ -105,49 +106,3 @@ struct RationalBase { }; using Rational = RationalBase; - -template -struct Binomial { - vector fc; - Binomial(int = 0) { fc.emplace_back(1); } - void extend() { - int n = fc.size(); - R nxt = fc.back() * n; - fc.push_back(nxt); - } - R fac(int n) { - if (n < 0) return 0; - while ((int)fc.size() <= n) extend(); - return fc[n]; - } - R finv(int n) { - if (n < 0) return 0; - return fac(n).inverse(); - } - R inv(int n) { - if (n < 0) return -inv(-n); - return R{1, max(n, 1)}; - } - R C(int n, int r) { - if (n < 0 or r < 0 or n < r) return R{0}; - return fac(n) * finv(n - r) * finv(r); - } - R operator()(int n, int r) { return C(n, r); } - template - R multinomial(const vector& r) { - static_assert(is_integral::value == true); - int n = 0; - for (auto& x : r) { - if (x < 0) return R{0}; - n += x; - } - R res = fac(n); - for (auto& x : r) res *= finv(x); - return res; - } - - template - R operator()(const vector& r) { - return multinomial(r); - } -}; diff --git a/matrix/determinant-arbitrary-mod.hpp b/matrix/determinant-arbitrary-mod.hpp new file mode 100644 index 000000000..a830cdc1f --- /dev/null +++ b/matrix/determinant-arbitrary-mod.hpp @@ -0,0 +1,29 @@ +#pragma once + +#include +using namespace std; + +template +mint determinant_arbitrary_mod(vector> a) { + int n = a.size(); + mint det = 1; + for (int j = 0; j < n; j++) { + for (int i = j; i < n; i++) { + if (a[i][j] == 0) continue; + if (i != j) swap(a[i], a[j]), det = -det; + break; + } + if (a[j][j] == 0) return 0; + for (int i = j + 1; i < n; i++) { + while (a[i][j] != 0) { + long long q = a[j][j].get() / a[i][j].get(); + mint c = -q; + for (int k = j; k < n; k++) a[j][k] += a[i][k] * c; + swap(a[i], a[j]), det = -det; + } + } + } + + for (int i = 0; i < n; i++) det *= a[i][i]; + return det; +} diff --git a/verify/verify-aoj-other/aoj-2171-bigrational.test.cpp b/verify/verify-aoj-other/aoj-2171-bigrational.test.cpp index 8c9b322f4..e968734c6 100644 --- a/verify/verify-aoj-other/aoj-2171-bigrational.test.cpp +++ b/verify/verify-aoj-other/aoj-2171-bigrational.test.cpp @@ -41,7 +41,9 @@ double calc(int N, int s, int t, vl q, vvi A) { } m[i][i] -= 1; } - auto ans = LinearEquation(m, b)[0]; + auto anss = LinearEquation(m, b); + assert(sz(anss) >= 1); + auto ans = anss[0]; if constexpr (is_same_v) { return ans[s]; } else { @@ -50,7 +52,7 @@ double calc(int N, int s, int t, vl q, vvi A) { } void test() { - rep(t, 100) { + rep(t, 500) { int N = rng(2, 30); int S = -1, T = -1; do { @@ -60,7 +62,7 @@ void test() { vl q(N); each(x, q) x = rng(0, 1); vvi A(N, vi(N)); - rep(i, N) rep(j, i) A[j][i] = A[i][j] = rng(0, 1) ? rng(1, 10) : 0; + rep(i, N) rep(j, i) A[j][i] = A[i][j] = rng(0, 10) ? rng(1, 10) : 0; double a1 = calc(N, S, T, q, A); double a2 = calc(N, S, T, q, A); double error = abs(a1 - a2) / max(1.0, a2); diff --git a/verify/verify-unit-test/bigrational.test.cpp b/verify/verify-unit-test/bigrational.test.cpp new file mode 100644 index 000000000..72a309a2c --- /dev/null +++ b/verify/verify-unit-test/bigrational.test.cpp @@ -0,0 +1,111 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/aplusb" + +#include "../../template/template.hpp" +// +#include "../../math/rational-binomial.hpp" +#include "../../math/rational-fps.hpp" +#include "../../math/rational.hpp" +// +#include "../../math/bigint-rational.hpp" +// +#include "../../modint/montgomery-modint.hpp" +#include "../../modulo/binomial.hpp" +using mint = LazyMontgomeryModInt<998244353>; +// +#include "../../misc/rng.hpp" + +using namespace Nyaan; + +void Nyaan::solve() { + { + BigRational a{4, 3}, b{2, 3}; + assert(a + b == 2); + assert(a - b == BigRational(2, 3)); + assert(a * b == BigRational(8, 9)); + assert(a / b == 2); + assert(a.inverse() == BigRational(3, 4)); + assert(a.pow(3) == BigRational(64, 27)); + + assert((a > b) == true); + assert((a >= b) == true); + assert((a < b) == false); + assert((a <= b) == false); + } + + { + Binomial_rational C; + assert(C.fac(3) == 6); + assert(C.finv(3) == BigRational(1, 6)); + assert(C(4, 2) == 6); + assert(C(vi{3, 2}) == 10); + } + + { + using fps = FormalPowerSeries_rational; + { + fps f{1, 2, {3, 2}}, g{{1, 4}, 5}; + fps h{{5, 4}, 7, {3, 2}}; + assert(f + g == h); + h = fps{{3, 4}, -3, {3, 2}}; + assert(f - g == h); + assert(f * g % g == fps{}); + assert(f * g % f == fps{}); + } + + { + fps e{1, 1, {1, 2}, {1, 6}, {1, 24}, {1, 120}}; + fps f = e.pow(TEN(10)); + trc(f); + rep(i, sz(e)) { assert(e[i] * BigRational{TEN(10)}.pow(i) == f[i]); } + } + } + + // mint と挙動の比較 + { + auto comp = [&](int i, int j, int k, int l) { + rep(b, 16) { + int ii = (b >> 0) % 2 ? -i : +i; + int jj = (b >> 1) % 2 ? -j : +j; + int kk = (b >> 2) % 2 ? -k : +k; + int ll = (b >> 3) % 2 ? -l : +l; + BigRational x{ii, jj}, y{kk, ll}; + mint X = mint{ii} / jj; + mint Y = mint{kk} / ll; + assert(X + Y == (x + y).to_mint(998244353).to_ll()); + assert(X - Y == (x - y).to_mint(998244353).to_ll()); + assert(X * Y == (x * y).to_mint(998244353).to_ll()); + if (Y != 0) { + assert(X / Y == (x / y).to_mint(998244353).to_ll()); + } + } + }; + rep(i, 10) rep1(j, 10) rep(k, 10) rep1(l, 10) comp(i, j, k, l); + rep(t, 1000) { + ll lower = t % 2 ? 1 : TEN(17); + ll i = rng(lower, TEN(18)); + ll j = rng(lower, TEN(18)); + ll k = rng(lower, TEN(18)); + ll l = rng(lower, TEN(18)); + comp(i, j, k, l); + } + } + + // binom, mint と挙動の比較 + { + Binomial_rational C1; + Binomial C2; + reg(i, -15, 15) { + assert(C2.fac(i) == C1.fac(i).to_mint(998244353).to_ll()); + assert(C2.finv(i) == C1.finv(i).to_mint(998244353).to_ll()); + assert(C2.inv(i) == C1.inv(i).to_mint(998244353).to_ll()); + reg(j, -15, 15) assert(C2(i, j) == C1(i, j).to_mint(998244353).to_ll()); + } + } + + cerr << "OK" << endl; + { + int s, t; + cin >> s >> t; + cout << s + t << "\n"; + } +} diff --git a/verify/verify-unit-test/rational-number.test.cpp b/verify/verify-unit-test/rational-number.test.cpp index 5c56b11a6..ff1ac1a50 100644 --- a/verify/verify-unit-test/rational-number.test.cpp +++ b/verify/verify-unit-test/rational-number.test.cpp @@ -1,87 +1,109 @@ #define PROBLEM "https://judge.yosupo.jp/problem/aplusb" #include "../../template/template.hpp" - // +#include "../../math/rational-binomial.hpp" +#include "../../math/rational-fps.hpp" #include "../../math/rational.hpp" -#include "../../misc/rng.hpp" +// #include "../../modint/montgomery-modint.hpp" -using mint = LazyMontgomeryModInt<998244353>; - -namespace mint_binom { - #include "../../modulo/binomial.hpp" - -} +using mint = LazyMontgomeryModInt<998244353>; +// +#include "../../misc/rng.hpp" using namespace Nyaan; void Nyaan::solve() { - Rational a{4, 3}, b{2, 3}; + { + Rational a{4, 3}, b{2, 3}; + assert(a + b == 2); + assert(a - b == Rational(2, 3)); + assert(a * b == Rational(8, 9)); + assert(a / b == 2); + assert(a.inverse() == Rational(3, 4)); + assert(a.pow(3) == Rational(64, 27)); + + assert((a > b) == true); + assert((a >= b) == true); + assert((a < b) == false); + assert((a <= b) == false); + } - trc(a + b); - assert(a + b == 2); - trc(a - b); - assert(a - b == Rational(2, 3)); - trc(a * b); - assert(a * b == Rational(8, 9)); - trc(a / b); - assert(a / b == 2); - trc(a.inverse()); - assert(a.inverse() == Rational(3, 4)); - trc(a.pow(3)); - assert(a.pow(3) == Rational(64, 27)); + { + Binomial_rational C; + assert(C.fac(3) == 6); + assert(C.finv(3) == Rational(1, 6)); + assert(C(4, 2) == 6); + assert(C(vi{3, 2}) == 10); + } - trc(a > b); - assert(a > b == true); - trc(a >= b); - assert(a >= b == true); - trc(a < b); - assert(a < b == false); - trc(a <= b); - assert(a <= b == false); + { + using fps = FormalPowerSeries_rational; - Binomial C; - assert(C.fac(3) == 6); - assert(C.finv(3) == Rational(1, 6)); - assert(C(4, 2) == 6); - assert(C(vi{3, 2}) == 10); + { + fps f{1, 2, {3, 2}}, g{{1, 4}, 5}; + fps h{{5, 4}, 7, {3, 2}}; + assert(f + g == h); + h = fps{{3, 4}, -3, {3, 2}}; + assert(f - g == h); + assert(f * g % g == fps{}); + assert(f * g % f == fps{}); + } - auto comp = [&](int i, int j, int k, int l) { - rep(b, 16) { - int ii = (b >> 0) % 2 ? -i : +i; - int jj = (b >> 1) % 2 ? -j : +j; - int kk = (b >> 2) % 2 ? -k : +k; - int ll = (b >> 3) % 2 ? -l : +l; - Rational x{ii, jj}, y{kk, ll}; - mint X = mint{ii} / jj; - mint Y = mint{kk} / ll; - assert(X + Y == (x + y).to_mint(998244353)); - assert(X - Y == (x - y).to_mint(998244353)); - assert(X * Y == (x * y).to_mint(998244353)); - if (Y != 0) { - assert(X / Y == (x / y).to_mint(998244353)); + + { + fps e{1, 1, {1, 2}, {1, 6}, {1, 24}, {1, 120}}; + fps f = e.pow(10); + trc(f); + rep(i, sz(e)) { + assert(e[i] * Rational{10}.pow(i) == f[i]); } } - }; - rep(i, 20) rep1(j, 20) rep(k, 20) rep1(l, 20) comp(i, j, k, l); - rep(t, 10000) { - int lower = t % 2 ? 1 : 32000; - ll i = rng(lower, 35000); - ll j = rng(lower, 35000); - ll k = rng(lower, 35000); - ll l = rng(lower, 35000); - comp(i, j, k, l); } - Binomial C1; - mint_binom::Binomial C2; - reg(i, -15, 15) { - assert(C2.fac(i) == C1.fac(i).to_mint(998244353)); - assert(C2.finv(i) == C1.finv(i).to_mint(998244353)); - assert(C2.inv(i) == C1.inv(i).to_mint(998244353)); - reg(j, -15, 15) assert(C2(i, j) == C1(i, j).to_mint(998244353)); + // mint と挙動の比較 + { + auto comp = [&](int i, int j, int k, int l) { + rep(b, 16) { + int ii = (b >> 0) % 2 ? -i : +i; + int jj = (b >> 1) % 2 ? -j : +j; + int kk = (b >> 2) % 2 ? -k : +k; + int ll = (b >> 3) % 2 ? -l : +l; + Rational x{ii, jj}, y{kk, ll}; + mint X = mint{ii} / jj; + mint Y = mint{kk} / ll; + assert(X + Y == (x + y).to_mint(998244353)); + assert(X - Y == (x - y).to_mint(998244353)); + assert(X * Y == (x * y).to_mint(998244353)); + if (Y != 0) { + assert(X / Y == (x / y).to_mint(998244353)); + } + } + }; + rep(i, 20) rep1(j, 20) rep(k, 20) rep1(l, 20) comp(i, j, k, l); + rep(t, 10000) { + int lower = t % 2 ? 1 : 32000; + ll i = rng(lower, 35000); + ll j = rng(lower, 35000); + ll k = rng(lower, 35000); + ll l = rng(lower, 35000); + comp(i, j, k, l); + } + } + + // binom, mint と挙動の比較 + { + Binomial_rational C1; + Binomial C2; + reg(i, -15, 15) { + assert(C2.fac(i) == C1.fac(i).to_mint(998244353)); + assert(C2.finv(i) == C1.finv(i).to_mint(998244353)); + assert(C2.inv(i) == C1.inv(i).to_mint(998244353)); + reg(j, -15, 15) assert(C2(i, j) == C1(i, j).to_mint(998244353)); + } } + cerr << "OK" << endl; { int s, t; diff --git a/verify/verify-yosupo-ds/yosupo-deque-operate-all-composite.test.cpp b/verify/verify-yosupo-ds/yosupo-deque-operate-all-composite.test.cpp new file mode 100644 index 000000000..3e35f49ec --- /dev/null +++ b/verify/verify-yosupo-ds/yosupo-deque-operate-all-composite.test.cpp @@ -0,0 +1,54 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/deque_operate_all_composite" +// +#include "../../template/template.hpp" +// +#include "../../data-structure/slide-window-aggregation-deque.hpp" +#include "../../math/affine-transformation.hpp" +// +#include "../../modint/montgomery-modint.hpp" +// +using namespace Nyaan; +using mint = LazyMontgomeryModInt<998244353>; +// using mint = LazyMontgomeryModInt<1000000007>; +using vm = vector; +using vvm = vector; +using namespace Nyaan; + +void q() { + using A = Affine; + + SlideWindowAggregationDeque swag([](const A &l, const A &r) { return l * r; }, + A{}); + deque dq; + + ini(Q); + while (Q--) { + ini(cmd); + if (cmd == 0) { + ini(a, b); + swag.push_front({a, b}); + dq.push_front({a, b}); + } else if (cmd == 1) { + ini(a, b); + swag.push_back({a, b}); + dq.push_back({a, b}); + } else if (cmd == 2) { + if (dq.front() != swag.front()) exit(1); + swag.pop_front(); + dq.pop_front(); + } else if (cmd == 3) { + if (dq.back() != swag.back()) exit(1); + swag.pop_back(); + dq.pop_back(); + } else { + ini(x); + out(swag.query()(x)); + } + } +} + +void Nyaan::solve() { + int t = 1; + // in(t); + while (t--) q(); +} diff --git a/verify/verify-yosupo-ds/yosupo-range-affine-point-get-2.test.cpp b/verify/verify-yosupo-ds/yosupo-range-affine-point-get-2.test.cpp new file mode 100644 index 000000000..045e0aa67 --- /dev/null +++ b/verify/verify-yosupo-ds/yosupo-range-affine-point-get-2.test.cpp @@ -0,0 +1,38 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/range_affine_point_get" +// +#include "../../template/template.hpp" +// +#include "../../modint/montgomery-modint.hpp" +#include "../../segment-tree/lazy-segment-tree-utility.hpp" + +using mint = LazyMontgomeryModInt<998244353>; +using Pair = pair; +Pair f(Pair a, Pair b) { return Pair{a.first + b.first, a.second + b.second}; }; +Pair g(Pair a, Pair b) { + return Pair{a.first * b.first + a.second * b.second, a.second}; +}; +Pair h(Pair a, Pair b) { + return Pair{a.first * b.first, a.second * b.first + b.second}; +}; +Pair ti() { return {0, 0}; } +Pair ei() { return {1, 0}; } + +using namespace Nyaan; +void Nyaan::solve() { + ini(N, Q); + V a(N, {0, 1}); + rep(i, N) in(a[i].first); + + LazySegmentTreeBase seg(a); + + rep(_, Q) { + ini(cmd); + if (!cmd) { + inl(l, r, b, c); + seg.update(l, r, {b, c}); + } else { + inl(i); + out(seg.get_val(i).first); + } + } +} diff --git a/verify/verify-yosupo-ds/yosupo-range-affine-point-get.test.cpp b/verify/verify-yosupo-ds/yosupo-range-affine-point-get.test.cpp new file mode 100644 index 000000000..bf73c7f26 --- /dev/null +++ b/verify/verify-yosupo-ds/yosupo-range-affine-point-get.test.cpp @@ -0,0 +1,38 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/range_affine_point_get" +// +#include "../../template/template.hpp" +// +#include "../../modint/montgomery-modint.hpp" +#include "../../segment-tree/lazy-segment-tree.hpp" + +using mint = LazyMontgomeryModInt<998244353>; + +using namespace Nyaan; +void Nyaan::solve() { + using P = pair; + auto f = [](P a, P b) { return P{a.first + b.first, a.second + b.second}; }; + auto g = [](P a, P b) { + return P{a.first * b.first + a.second * b.second, a.second}; + }; + auto h = [](P a, P b) { + return P{a.first * b.first, a.second * b.first + b.second}; + }; + P ti = {0, 0}; + P ei = {1, 0}; + ini(N, Q); + V

a(N, {0, 1}); + rep(i, N) in(a[i].first); + + LazySegmentTree seg(a, f, g, h, ti, ei); + + rep(_, Q) { + ini(cmd); + if (!cmd) { + inl(l, r, b, c); + seg.update(l, r, {b, c}); + } else { + inl(i); + out(seg.get_val(i).first); + } + } +} diff --git a/verify/verify-yosupo-math/yosupo-determinant-arbitrary-mod.test.cpp b/verify/verify-yosupo-math/yosupo-determinant-arbitrary-mod.test.cpp new file mode 100644 index 000000000..075ac904e --- /dev/null +++ b/verify/verify-yosupo-math/yosupo-determinant-arbitrary-mod.test.cpp @@ -0,0 +1,23 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/matrix_det_arbitrary_mod" +// +#include "../../template/template.hpp" +// +#include "../../matrix/determinant-arbitrary-mod.hpp" +// +#include "../../modint/arbitrary-modint.hpp" +using mint = ArbitraryModInt; +using namespace Nyaan; + +void q() { + inl(N, mod); + mint::set_mod(mod); + vector a(N, vector(N, mint{})); + in(a); + out(determinant_arbitrary_mod(a)); +} + +void Nyaan::solve() { + int t = 1; + // in(t); + while (t--) q(); +} diff --git a/verify/verify-yosupo-math/yosupo-primitive-root.test.cpp b/verify/verify-yosupo-math/yosupo-primitive-root.test.cpp new file mode 100644 index 000000000..b94fef9de --- /dev/null +++ b/verify/verify-yosupo-math/yosupo-primitive-root.test.cpp @@ -0,0 +1,21 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/primitive_root" +// +#include "../../template/template.hpp" +// +#include "../../math/primitive-root-ll.hpp" + +using namespace Nyaan; + +void q() { + ini(Q); + while (Q--) { + inl(p); + out(primitive_root_ll(p)); + } +} + +void Nyaan::solve() { + int t = 1; + // in(t); + while (t--) q(); +} diff --git a/verify/verify-yuki/yuki-1939-sparse-pow.test.cpp b/verify/verify-yuki/yuki-1939-sparse-pow.test.cpp new file mode 100644 index 000000000..d05b301e8 --- /dev/null +++ b/verify/verify-yuki/yuki-1939-sparse-pow.test.cpp @@ -0,0 +1,34 @@ +#define PROBLEM "https://yukicoder.me/problems/no/1939" +// +#include "../../template/template.hpp" +// +#include "../../fps/sparse-fps.hpp" +// +#include "../../fps/ntt-friendly-fps.hpp" +#include "../../modint/montgomery-modint.hpp" +#include "../../modulo/binomial.hpp" +// +using namespace Nyaan; +using mint = LazyMontgomeryModInt<998244353>; +using vm = vector; +using vvm = vector; +Binomial C; +using fps = FormalPowerSeries; +// +using namespace Nyaan; + +void q() { + inl(N, M); + vl L(M); + in(L); + fps g(N + 1); + g[0] = 1; + each(n, L) g[n] = 1; + out(sparse_pow(g, N + 1)[N] / (N + 1)); +} + +void Nyaan::solve() { + int t = 1; + // in(t); + while (t--) q(); +}