|
| 1 | +--- |
| 2 | +data: |
| 3 | + _extendedDependsOn: |
| 4 | + - icon: ':heavy_check_mark:' |
| 5 | + path: linear_algebra_matrix/upper_triangular_matrix.hpp |
| 6 | + title: "Upper triangular matrix \uFF08\u5B9A\u6570\u6B21\u5143\u4E0A\u4E09\u89D2\ |
| 7 | + \u884C\u5217\uFF09" |
| 8 | + - icon: ':heavy_check_mark:' |
| 9 | + path: modint.hpp |
| 10 | + title: modint.hpp |
| 11 | + _extendedRequiredBy: [] |
| 12 | + _extendedVerifiedWith: [] |
| 13 | + _isVerificationFailed: false |
| 14 | + _pathExtension: cpp |
| 15 | + _verificationStatusIcon: ':heavy_check_mark:' |
| 16 | + attributes: |
| 17 | + '*NOT_SPECIAL_COMMENTS*': '' |
| 18 | + PROBLEM: https://yukicoder.me/problems/no/3530 |
| 19 | + links: |
| 20 | + - https://yukicoder.me/problems/no/3530 |
| 21 | + bundledCode: "#line 1 \"linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp\"\ |
| 22 | + \n#define PROBLEM \"https://yukicoder.me/problems/no/3530\"\n#line 2 \"linear_algebra_matrix/upper_triangular_matrix.hpp\"\ |
| 23 | + \n\ntemplate <class T> struct UpperTriangular3d {\n static T explicit_init_required()\ |
| 24 | + \ = delete;\n T a00 = this->explicit_init_required(), a01 = this->explicit_init_required(),\n\ |
| 25 | + \ a02 = this->explicit_init_required();\n T a11 = this->explicit_init_required(),\ |
| 26 | + \ a12 = this->explicit_init_required();\n T a22 = this->explicit_init_required();\n\ |
| 27 | + \n UpperTriangular3d operator*(const UpperTriangular3d &r) const {\n \ |
| 28 | + \ return UpperTriangular3d{\n .a00 = this->a00 * r.a00,\n \ |
| 29 | + \ .a01 = this->a00 * r.a01 + this->a01 * r.a11,\n .a02 = this->a00\ |
| 30 | + \ * r.a02 + this->a01 * r.a12 + this->a02 * r.a22,\n .a11 = this->a11\ |
| 31 | + \ * r.a11,\n .a12 = this->a11 * r.a12 + this->a12 * r.a22,\n \ |
| 32 | + \ .a22 = this->a22 * r.a22,\n };\n }\n\n UpperTriangular3d operator-()\ |
| 33 | + \ const {\n return UpperTriangular3d{\n .a00 = -this->a00,\n\ |
| 34 | + \ .a01 = -this->a01,\n .a02 = -this->a02,\n .a11\ |
| 35 | + \ = -this->a11,\n .a12 = -this->a12,\n .a22 = -this->a22,\n\ |
| 36 | + \ };\n }\n\n UpperTriangular3d operator+(const UpperTriangular3d\ |
| 37 | + \ &r) const {\n return UpperTriangular3d{\n .a00 = this->a00\ |
| 38 | + \ + r.a00,\n .a01 = this->a01 + r.a01,\n .a02 = this->a02\ |
| 39 | + \ + r.a02,\n .a11 = this->a11 + r.a11,\n .a12 = this->a12\ |
| 40 | + \ + r.a12,\n .a22 = this->a22 + r.a22,\n };\n }\n\n auto\ |
| 41 | + \ operator<=>(const UpperTriangular3d &) const = default;\n};\n#line 2 \"modint.hpp\"\ |
| 42 | + \n#include <cassert>\n#include <iostream>\n#include <set>\n#include <vector>\n\ |
| 43 | + \ntemplate <int md> struct ModInt {\n static_assert(md > 1);\n using lint\ |
| 44 | + \ = long long;\n constexpr static int mod() { return md; }\n static int\ |
| 45 | + \ get_primitive_root() {\n static int primitive_root = 0;\n if (!primitive_root)\ |
| 46 | + \ {\n primitive_root = [&]() {\n std::set<int> fac;\n\ |
| 47 | + \ int v = md - 1;\n for (lint i = 2; i * i <= v;\ |
| 48 | + \ i++)\n while (v % i == 0) fac.insert(i), v /= i;\n \ |
| 49 | + \ if (v > 1) fac.insert(v);\n for (int g = 1; g < md;\ |
| 50 | + \ g++) {\n bool ok = true;\n for (auto i\ |
| 51 | + \ : fac)\n if (ModInt(g).pow((md - 1) / i) == 1) {\n \ |
| 52 | + \ ok = false;\n break;\n \ |
| 53 | + \ }\n if (ok) return g;\n \ |
| 54 | + \ }\n return -1;\n }();\n }\n return\ |
| 55 | + \ primitive_root;\n }\n int val_;\n int val() const noexcept { return\ |
| 56 | + \ val_; }\n constexpr ModInt() : val_(0) {}\n constexpr ModInt &_setval(lint\ |
| 57 | + \ v) { return val_ = (v >= md ? v - md : v), *this; }\n constexpr ModInt(lint\ |
| 58 | + \ v) { _setval(v % md + md); }\n constexpr explicit operator bool() const {\ |
| 59 | + \ return val_ != 0; }\n constexpr ModInt operator+(const ModInt &x) const {\n\ |
| 60 | + \ return ModInt()._setval((lint)val_ + x.val_);\n }\n constexpr ModInt\ |
| 61 | + \ operator-(const ModInt &x) const {\n return ModInt()._setval((lint)val_\ |
| 62 | + \ - x.val_ + md);\n }\n constexpr ModInt operator*(const ModInt &x) const\ |
| 63 | + \ {\n return ModInt()._setval((lint)val_ * x.val_ % md);\n }\n constexpr\ |
| 64 | + \ ModInt operator/(const ModInt &x) const {\n return ModInt()._setval((lint)val_\ |
| 65 | + \ * x.inv().val() % md);\n }\n constexpr ModInt operator-() const { return\ |
| 66 | + \ ModInt()._setval(md - val_); }\n constexpr ModInt &operator+=(const ModInt\ |
| 67 | + \ &x) { return *this = *this + x; }\n constexpr ModInt &operator-=(const ModInt\ |
| 68 | + \ &x) { return *this = *this - x; }\n constexpr ModInt &operator*=(const ModInt\ |
| 69 | + \ &x) { return *this = *this * x; }\n constexpr ModInt &operator/=(const ModInt\ |
| 70 | + \ &x) { return *this = *this / x; }\n friend constexpr ModInt operator+(lint\ |
| 71 | + \ a, const ModInt &x) { return ModInt(a) + x; }\n friend constexpr ModInt operator-(lint\ |
| 72 | + \ a, const ModInt &x) { return ModInt(a) - x; }\n friend constexpr ModInt operator*(lint\ |
| 73 | + \ a, const ModInt &x) { return ModInt(a) * x; }\n friend constexpr ModInt operator/(lint\ |
| 74 | + \ a, const ModInt &x) { return ModInt(a) / x; }\n constexpr bool operator==(const\ |
| 75 | + \ ModInt &x) const { return val_ == x.val_; }\n constexpr bool operator!=(const\ |
| 76 | + \ ModInt &x) const { return val_ != x.val_; }\n constexpr bool operator<(const\ |
| 77 | + \ ModInt &x) const {\n return val_ < x.val_;\n } // To use std::map<ModInt,\ |
| 78 | + \ T>\n friend std::istream &operator>>(std::istream &is, ModInt &x) {\n \ |
| 79 | + \ lint t;\n return is >> t, x = ModInt(t), is;\n }\n constexpr\ |
| 80 | + \ friend std::ostream &operator<<(std::ostream &os, const ModInt &x) {\n \ |
| 81 | + \ return os << x.val_;\n }\n\n constexpr ModInt pow(lint n) const {\n\ |
| 82 | + \ ModInt ans = 1, tmp = *this;\n while (n) {\n if (n\ |
| 83 | + \ & 1) ans *= tmp;\n tmp *= tmp, n >>= 1;\n }\n return\ |
| 84 | + \ ans;\n }\n\n static constexpr int cache_limit = std::min(md, 1 << 21);\n\ |
| 85 | + \ static std::vector<ModInt> facs, facinvs, invs;\n\n constexpr static void\ |
| 86 | + \ _precalculation(int N) {\n const int l0 = facs.size();\n if (N\ |
| 87 | + \ > md) N = md;\n if (N <= l0) return;\n facs.resize(N), facinvs.resize(N),\ |
| 88 | + \ invs.resize(N);\n for (int i = l0; i < N; i++) facs[i] = facs[i - 1]\ |
| 89 | + \ * i;\n facinvs[N - 1] = facs.back().pow(md - 2);\n for (int i\ |
| 90 | + \ = N - 2; i >= l0; i--) facinvs[i] = facinvs[i + 1] * (i + 1);\n for (int\ |
| 91 | + \ i = N - 1; i >= l0; i--) invs[i] = facinvs[i] * facs[i - 1];\n }\n\n constexpr\ |
| 92 | + \ ModInt inv() const {\n if (this->val_ < cache_limit) {\n if\ |
| 93 | + \ (facs.empty()) facs = {1}, facinvs = {1}, invs = {0};\n while (this->val_\ |
| 94 | + \ >= int(facs.size())) _precalculation(facs.size() * 2);\n return invs[this->val_];\n\ |
| 95 | + \ } else {\n return this->pow(md - 2);\n }\n }\n\n\ |
| 96 | + \ constexpr static ModInt fac(int n) {\n assert(n >= 0);\n if\ |
| 97 | + \ (n >= md) return ModInt(0);\n while (n >= int(facs.size())) _precalculation(facs.size()\ |
| 98 | + \ * 2);\n return facs[n];\n }\n\n constexpr static ModInt facinv(int\ |
| 99 | + \ n) {\n assert(n >= 0);\n if (n >= md) return ModInt(0);\n \ |
| 100 | + \ while (n >= int(facs.size())) _precalculation(facs.size() * 2);\n \ |
| 101 | + \ return facinvs[n];\n }\n\n constexpr static ModInt doublefac(int n) {\n\ |
| 102 | + \ assert(n >= 0);\n if (n >= md) return ModInt(0);\n long\ |
| 103 | + \ long k = (n + 1) / 2;\n return (n & 1) ? ModInt::fac(k * 2) / (ModInt(2).pow(k)\ |
| 104 | + \ * ModInt::fac(k))\n : ModInt::fac(k) * ModInt(2).pow(k);\n\ |
| 105 | + \ }\n\n constexpr static ModInt nCr(int n, int r) {\n assert(n >=\ |
| 106 | + \ 0);\n if (r < 0 or n < r) return ModInt(0);\n return ModInt::fac(n)\ |
| 107 | + \ * ModInt::facinv(r) * ModInt::facinv(n - r);\n }\n\n constexpr static\ |
| 108 | + \ ModInt nPr(int n, int r) {\n assert(n >= 0);\n if (r < 0 or n\ |
| 109 | + \ < r) return ModInt(0);\n return ModInt::fac(n) * ModInt::facinv(n - r);\n\ |
| 110 | + \ }\n\n static ModInt binom(long long n, long long r) {\n static\ |
| 111 | + \ long long bruteforce_times = 0;\n\n if (r < 0 or n < r) return ModInt(0);\n\ |
| 112 | + \ if (n <= bruteforce_times or n < (int)facs.size()) return ModInt::nCr(n,\ |
| 113 | + \ r);\n\n r = std::min(r, n - r);\n assert((int)r == r);\n\n \ |
| 114 | + \ ModInt ret = ModInt::facinv(r);\n for (int i = 0; i < r; ++i) ret\ |
| 115 | + \ *= n - i;\n bruteforce_times += r;\n\n return ret;\n }\n\n\ |
| 116 | + \ // Multinomial coefficient, (k_1 + k_2 + ... + k_m)! / (k_1! k_2! ... k_m!)\n\ |
| 117 | + \ // Complexity: O(sum(ks))\n // Verify: https://yukicoder.me/problems/no/3178\n\ |
| 118 | + \ template <class Vec> static ModInt multinomial(const Vec &ks) {\n \ |
| 119 | + \ ModInt ret{1};\n int sum = 0;\n for (int k : ks) {\n \ |
| 120 | + \ assert(k >= 0);\n ret *= ModInt::facinv(k), sum += k;\n \ |
| 121 | + \ }\n return ret * ModInt::fac(sum);\n }\n template <class... Args>\ |
| 122 | + \ static ModInt multinomial(Args... args) {\n int sum = (0 + ... + args);\n\ |
| 123 | + \ ModInt result = (1 * ... * ModInt::facinv(args));\n return ModInt::fac(sum)\ |
| 124 | + \ * result;\n }\n\n // Catalan number, C_n = binom(2n, n) / (n + 1) = #\ |
| 125 | + \ of Dyck words of length 2n\n // C_0 = 1, C_1 = 1, C_2 = 2, C_3 = 5, C_4 =\ |
| 126 | + \ 14, ...\n // https://oeis.org/A000108\n // Complexity: O(n)\n static\ |
| 127 | + \ ModInt catalan(int n) {\n if (n < 0) return ModInt(0);\n return\ |
| 128 | + \ ModInt::fac(n * 2) * ModInt::facinv(n + 1) * ModInt::facinv(n);\n }\n\n \ |
| 129 | + \ ModInt sqrt() const {\n if (val_ == 0) return 0;\n if (md ==\ |
| 130 | + \ 2) return val_;\n if (pow((md - 1) / 2) != 1) return 0;\n ModInt\ |
| 131 | + \ b = 1;\n while (b.pow((md - 1) / 2) == 1) b += 1;\n int e = 0,\ |
| 132 | + \ m = md - 1;\n while (m % 2 == 0) m >>= 1, e++;\n ModInt x = pow((m\ |
| 133 | + \ - 1) / 2), y = (*this) * x * x;\n x *= (*this);\n ModInt z = b.pow(m);\n\ |
| 134 | + \ while (y != 1) {\n int j = 0;\n ModInt t = y;\n\ |
| 135 | + \ while (t != 1) j++, t *= t;\n z = z.pow(1LL << (e - j\ |
| 136 | + \ - 1));\n x *= z, z *= z, y *= z;\n e = j;\n }\n\ |
| 137 | + \ return ModInt(std::min(x.val_, md - x.val_));\n }\n};\ntemplate <int\ |
| 138 | + \ md> std::vector<ModInt<md>> ModInt<md>::facs = {1};\ntemplate <int md> std::vector<ModInt<md>>\ |
| 139 | + \ ModInt<md>::facinvs = {1};\ntemplate <int md> std::vector<ModInt<md>> ModInt<md>::invs\ |
| 140 | + \ = {0};\n\nusing ModInt998244353 = ModInt<998244353>;\n// using mint = ModInt<998244353>;\n\ |
| 141 | + // using mint = ModInt<1000000007>;\n#line 4 \"linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp\"\ |
| 142 | + \n\n#include <algorithm>\n#line 7 \"linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp\"\ |
| 143 | + \n#include <map>\n#include <utility>\n#line 10 \"linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp\"\ |
| 144 | + \n#include <atcoder/segtree>\n\nusing namespace std;\nusing S = UpperTriangular3d<ModInt998244353>;\n\ |
| 145 | + \nS op(const S &l, const S &r) { return l * r; }\nS e() { return S{1, 0, 0, 1,\ |
| 146 | + \ 0, 1}; }\n\nS GenR() { return S{ModInt998244353(3) / 4, ModInt998244353(1) /\ |
| 147 | + \ 4, 0, 1, 0, 1}; }\nS GenL() { return S{1, 0, 0, ModInt998244353(3) / 4, ModInt998244353(1)\ |
| 148 | + \ / 4, 1}; }\n\nModInt998244353 Solve(vector<pair<int, int>> ps) {\n vector<tuple<int,\ |
| 149 | + \ int, int>> yxis;\n for (int i = 0; i < (int)ps.size(); ++i) {\n auto\ |
| 150 | + \ [x, y] = ps.at(i);\n yxis.emplace_back(y, x, i);\n }\n sort(yxis.begin(),\ |
| 151 | + \ yxis.end());\n\n const auto L = GenL(), R = GenR();\n const vector<S>\ |
| 152 | + \ init(ps.size(), R);\n atcoder::segtree<S, op, e> seg(init);\n\n int last_x\ |
| 153 | + \ = -1e9;\n\n ModInt998244353 ret = 0;\n\n map<int, vector<pair<int, int>>>\ |
| 154 | + \ x2yis;\n for (int i = 0; i < (int)ps.size(); ++i) {\n auto [x, y]\ |
| 155 | + \ = ps.at(i);\n x2yis[x].emplace_back(y, i);\n }\n for (auto [x,\ |
| 156 | + \ yis] : x2yis) {\n\n const ModInt998244353 dx = x - last_x;\n ret\ |
| 157 | + \ += dx * seg.all_prod().a02;\n\n for (auto [y, i] : yis) {\n \ |
| 158 | + \ const int idx =\n lower_bound(yxis.begin(), yxis.end(), make_tuple(y,\ |
| 159 | + \ x, i)) - yxis.begin();\n seg.set(idx, L);\n }\n\n last_x\ |
| 160 | + \ = x;\n }\n\n return ret;\n}\n\nint main() {\n int N;\n cin >> N;\n\ |
| 161 | + \ vector<pair<int, int>> xy(N);\n\n for (auto &[x, y] : xy) cin >> x >>\ |
| 162 | + \ y;\n\n const ModInt998244353 coeff = ModInt998244353(4).pow(N);\n ModInt998244353\ |
| 163 | + \ ret1 = Solve(xy) * coeff;\n for (auto &[x, y] : xy) swap(x, y);\n ModInt998244353\ |
| 164 | + \ ret2 = Solve(xy) * coeff;\n cout << (ret1 + ret2) * 2 << '\\n';\n}\n" |
| 165 | + code: "#define PROBLEM \"https://yukicoder.me/problems/no/3530\"\n#include \"../upper_triangular_matrix.hpp\"\ |
| 166 | + \n#include \"../../modint.hpp\"\n\n#include <algorithm>\n#include <iostream>\n\ |
| 167 | + #include <map>\n#include <utility>\n#include <vector>\n#include <atcoder/segtree>\n\ |
| 168 | + \nusing namespace std;\nusing S = UpperTriangular3d<ModInt998244353>;\n\nS op(const\ |
| 169 | + \ S &l, const S &r) { return l * r; }\nS e() { return S{1, 0, 0, 1, 0, 1}; }\n\ |
| 170 | + \nS GenR() { return S{ModInt998244353(3) / 4, ModInt998244353(1) / 4, 0, 1, 0,\ |
| 171 | + \ 1}; }\nS GenL() { return S{1, 0, 0, ModInt998244353(3) / 4, ModInt998244353(1)\ |
| 172 | + \ / 4, 1}; }\n\nModInt998244353 Solve(vector<pair<int, int>> ps) {\n vector<tuple<int,\ |
| 173 | + \ int, int>> yxis;\n for (int i = 0; i < (int)ps.size(); ++i) {\n auto\ |
| 174 | + \ [x, y] = ps.at(i);\n yxis.emplace_back(y, x, i);\n }\n sort(yxis.begin(),\ |
| 175 | + \ yxis.end());\n\n const auto L = GenL(), R = GenR();\n const vector<S>\ |
| 176 | + \ init(ps.size(), R);\n atcoder::segtree<S, op, e> seg(init);\n\n int last_x\ |
| 177 | + \ = -1e9;\n\n ModInt998244353 ret = 0;\n\n map<int, vector<pair<int, int>>>\ |
| 178 | + \ x2yis;\n for (int i = 0; i < (int)ps.size(); ++i) {\n auto [x, y]\ |
| 179 | + \ = ps.at(i);\n x2yis[x].emplace_back(y, i);\n }\n for (auto [x,\ |
| 180 | + \ yis] : x2yis) {\n\n const ModInt998244353 dx = x - last_x;\n ret\ |
| 181 | + \ += dx * seg.all_prod().a02;\n\n for (auto [y, i] : yis) {\n \ |
| 182 | + \ const int idx =\n lower_bound(yxis.begin(), yxis.end(), make_tuple(y,\ |
| 183 | + \ x, i)) - yxis.begin();\n seg.set(idx, L);\n }\n\n last_x\ |
| 184 | + \ = x;\n }\n\n return ret;\n}\n\nint main() {\n int N;\n cin >> N;\n\ |
| 185 | + \ vector<pair<int, int>> xy(N);\n\n for (auto &[x, y] : xy) cin >> x >>\ |
| 186 | + \ y;\n\n const ModInt998244353 coeff = ModInt998244353(4).pow(N);\n ModInt998244353\ |
| 187 | + \ ret1 = Solve(xy) * coeff;\n for (auto &[x, y] : xy) swap(x, y);\n ModInt998244353\ |
| 188 | + \ ret2 = Solve(xy) * coeff;\n cout << (ret1 + ret2) * 2 << '\\n';\n}\n" |
| 189 | + dependsOn: |
| 190 | + - linear_algebra_matrix/upper_triangular_matrix.hpp |
| 191 | + - modint.hpp |
| 192 | + isVerificationFile: true |
| 193 | + path: linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp |
| 194 | + requiredBy: [] |
| 195 | + timestamp: '2026-05-06 21:05:36+09:00' |
| 196 | + verificationStatus: TEST_ACCEPTED |
| 197 | + verifiedWith: [] |
| 198 | +documentation_of: linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp |
| 199 | +layout: document |
| 200 | +redirect_from: |
| 201 | +- /verify/linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp |
| 202 | +- /verify/linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp.html |
| 203 | +title: linear_algebra_matrix/test/upper_trinaglular_matrix.yuki3530.test.cpp |
| 204 | +--- |
0 commit comments