#line 1 "test/matrix/gaussian_elimination/yuki129.test.cpp" #define PROBLEM "https://judge.yosupo.jp/problem/system_of_linear_equations" #line 2 "src/math/Modint.hpp" #define CUT #line 2 "src/Template.hpp" #define CUT #include using namespace std; #define rep(i, l, r) for (int i = (l); i < (r); ++i) #define rrep(i, l, r) for (int i = (r); i --> (l);) #define all(c) begin(c), end(c) #ifdef LOCAL #define debug(...) debug_impl(#__VA_ARGS__, __VA_ARGS__) template void debug_impl(string s, H&& h, Ts&&... ts) { cerr << '(' << s << "): (" << forward(h); ((cerr << ", " << forward(ts)), ..., (cerr << ")\n")); } #else #define debug(...) void(0) #endif template bool chmax(T& a, const T& b) { return b > a ? (a = b, true) : false; } template bool chmin(T& a, const T& b) { return b < a ? (a = b, true) : false; } template istream& operator>>(istream& in, vector& v) { for (auto& e : v) in >> e; return in; } template void read(Args&... args) { (cin >> ... >> args); } template ostream& operator<<(ostream& out, const vector& v) { int n = v.size(); rep(i, 0, n) { out << v[i]; if (i + 1 != n) out << ' '; } return out; } template void print(H&& h, Ts &&... ts) { cout << h, ((cout << ' ' << forward(ts)), ..., (cout << '\n')); } struct io_setup_ { io_setup_() { ios::sync_with_stdio(false), cin.tie(nullptr); cout << fixed << setprecision(10); } } io_setup{}; #undef CUT #define NOTE compile command: \texttt{g++ -std=gnu++17 -Wall -Wextra -g -fsanitize=address -fsanitize=undefined \$\{file\} -o \$\{fileDirname\}/\$\{fileBasenameNoExtension\}} #undef NOTE #define NOTE \texttt{-DLOCAL} を加えると \texttt{debug(...)} による出力が有効となる #undef NOTE #line 3 "src/math/ExtGCD.hpp" #define CUT constexpr long safe_mod(long long x, long long m) { return (x %= m) < 0 ? x + m : x; } // Returns `(x,g)` s.t. `g=\gcd(a,b)`, `xa\equiv g \pmod{b}`, and `0\leq x< \frac{b}{g}` constexpr pair inv_gcd(long long a, long long b) { assert(b > 0); a = safe_mod(a, b); if (a == 0) return { 0, b }; long long s = b, t = a, x = 0, y = 1, tmp = 0; while (t) { long long u = s / t; s -= t * u; x -= y * u; tmp = s, s = t, t = tmp; tmp = x, x = y, y = tmp; } if (x < 0) x += b / s; return { x, s }; } // Returns `x` s.t. `xa\equiv 1 \pmod{m}`. // Requirement: `\gcd(a, m) = 1`. constexpr long long inv_mod(long long a, long long m) { auto [x, g] = inv_gcd(a, m); assert(g == 1); return x; } // Returns `(x_0,y_0,g)` s.t. `g=\gcd(a,b)`, `ax_0 + by_0 = g`, and `0\leq x_0<\frac{b}{g}`. // 一般解は `(x,y)=(x_0+k\cdot\dfrac{b}{g}, y_0-k\cdot\dfrac{a}{g})\;(k\in\mathbb{Z})` なので、`(x_0,y_0)` は `x` が非負の下で最小の解 constexpr tuple ext_gcd(long long a, long long b) { auto [x, g] = inv_gcd(a, b); return { x, (g - x * a) / b, g }; } #undef CUT #line 4 "src/math/Modint.hpp" constexpr long long pow_mod(long long x, long long b, int m) { long long p = safe_mod(x, m), r = 1; for (; b; b >>= 1) { if (b & 1) (r *= p) %= m; (p *= p) %= m; } return r; } template struct modint { unsigned x; modint(): x(0) {} modint(long long v): x(safe_mod(v, mod)) {} int val() const { return x; } static modint raw(unsigned v) { modint res; res.x = v; return v; } modint operator-() const { return raw(x ? mod - x : 0); } modint& operator+=(modint t) { if ((x += t.x) >= mod) x -= mod; return *this; } modint& operator-=(modint t) { if ((x += mod - t.x) >= mod) x -= mod; return *this; } modint& operator*=(modint t) { x = (unsigned long long) x * t.x % mod; return *this; } modint& operator/=(modint t) { return *this *= t.inv(); } friend modint operator+(modint x, modint y) { return x += y; } friend modint operator-(modint x, modint y) { return x -= y; } friend modint operator*(modint x, modint y) { return x *= y; } friend modint operator/(modint x, modint y) { return x /= y; } friend bool operator==(modint x, modint y) { return x.x == y.x; } friend bool operator!=(modint x, modint y) { return x.x != y.x; } modint inv() const { return inv_mod(x, mod); } modint pow(long long b) const { assert(b >= 0); return pow_mod(x, b, mod); } }; #undef CUT #line 3 "src/matrix/GaussianElimination.hpp" #define CUT template struct GaussianElimination { vector sol; vector> bases; // 解が一つも存在しないかどうか bool empty = false; GaussianElimination(vector> A, const vector& b) { int n = A.size(); for (int i = 0; i < n; ++i) A[i].push_back(b[i]); solve(A); } private: static constexpr bool is_fp = is_floating_point_v; static bool is_zero(T x) { if constexpr (is_fp) { return abs(x) < 1e-9; } else { return x == 0; } } void solve(vector>& Ab) { const int n = Ab.size(), m = Ab[0].size() - 1; auto pivoting = [&](int l, int k) { int pos = m, res = n; // 浮動小数点数型のとき T max_val = 0; rep(i, k, n) { const auto& v = Ab[i]; // 浮動小数点数型のとき if constexpr (is_fp) { if (pos < m and abs(v[pos]) > max_val) { res = i, max_val = abs(v[pos]); } } rep(j, l, pos) { if (not is_zero(Ab[k][j])) { pos = j, res = i; // 浮動小数点数型のとき if constexpr (is_fp) { max_val = abs(Ab[i][j]); } break; } } } return pair{ pos, res }; }; int l = 0; rep(i, 0, n) { auto [mse, k] = pivoting(l, i); l = mse + 1; if (k == n) break; Ab[i].swap(Ab[k]); T cinv = T{ 1 } / Ab[i][mse]; rep(row, i + 1, n) if (not is_zero(Ab[row][mse])) { T c = Ab[row][mse] * cinv; rep(col, mse, m + 1) Ab[row][col] -= c * Ab[i][col]; } } int basis_num = m; vector down(m); sol.assign(m, T{0}); rrep(i, 0, n) { int mse = m + 1; rep(col, 0, m + 1) if (not is_zero(Ab[i][col])) { mse = col; break; } if (mse < m) { T cinv = T{ 1 } / Ab[i][mse]; rep(row, 0, i) if (not is_zero(Ab[row][mse])) { T c = Ab[row][mse] * cinv; rep(col, mse, m + 1) Ab[row][col] -= c * Ab[i][col]; } rep(col, mse, m + 1) Ab[i][col] *= cinv; sol[mse] = Ab[i][m]; down[mse] = true; --basis_num; } else if (mse == m) { empty = true; return; } } bases.assign(basis_num, vector(m)); int id = 0; rep(j, 0, m) if (not down[j]) { int i = 0; rep(j2, 0, m) bases[id][j2] = down[j2] ? Ab[i++][j] : 0; bases[id++][j] = -1; } } }; #undef CUT #line 5 "test/matrix/gaussian_elimination/yuki129.test.cpp" int main() { int k; read(k); vector A(k, vector(k)); vector b(k, 1); rep(i, 0, k) { rep(v, 1, 7) { if (i + v == k) { } else { int j = i + v; if (j > k) j = 0; A[i][j] -= 1. / 6.; } } A[i][i] += 1; } GaussianElimination gauss(A, b); print(gauss.sol[0]); return 0; }