結果
問題 | No.813 ユキちゃんの冒険 |
ユーザー |
|
提出日時 | 2023-06-07 22:37:22 |
言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 1,247 ms / 2,000 ms |
コード長 | 15,768 bytes |
コンパイル時間 | 4,372 ms |
コンパイル使用メモリ | 258,328 KB |
最終ジャッジ日時 | 2025-02-13 23:17:38 |
ジャッジサーバーID (参考情報) |
judge2 / judge5 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
other | AC * 26 |
ソースコード
#ifndef HIDDEN_IN_VS // 折りたたみ用// 警告の抑制#define _CRT_SECURE_NO_WARNINGS// ライブラリの読み込み#include <bits/stdc++.h>using namespace std;// 型名の短縮using ll = long long; // -2^63 ~ 2^63 = 9 * 10^18(int は -2^31 ~ 2^31 = 2 * 10^9)using pii = pair<int, int>; using pll = pair<ll, ll>; using pil = pair<int, ll>; using pli = pair<ll, int>;using vi = vector<int>; using vvi = vector<vi>; using vvvi = vector<vvi>;using vl = vector<ll>; using vvl = vector<vl>; using vvvl = vector<vvl>;using vb = vector<bool>; using vvb = vector<vb>; using vvvb = vector<vvb>;using vc = vector<char>; using vvc = vector<vc>; using vvvc = vector<vvc>;using vd = vector<double>; using vvd = vector<vd>; using vvvd = vector<vvd>;template <class T> using priority_queue_rev = priority_queue<T, vector<T>, greater<T>>;using Graph = vvi;// 定数の定義const double PI = acos(-1);const vi DX = { 1, 0, -1, 0 }; // 4 近傍(下,右,上,左)const vi DY = { 0, 1, 0, -1 };int INF = 1001001001; ll INFL = 4004004004004004004LL;double EPS = 1e-15;// 入出力高速化struct fast_io { fast_io() { cin.tie(nullptr); ios::sync_with_stdio(false); cout << fixed << setprecision(18); } } fastIOtmp;// 汎用マクロの定義#define all(a) (a).begin(), (a).end()#define sz(x) ((int)(x).size())#define lbpos(a, x) (int)distance((a).begin(), std::lower_bound(all(a), x))#define ubpos(a, x) (int)distance((a).begin(), std::upper_bound(all(a), x))#define Yes(b) {cout << ((b) ? "Yes\n" : "No\n");}#define YES(b) {cout << ((b) ? "YES\n" : "NO\n");}#define rep(i, n) for(int i = 0, i##_len = int(n); i < i##_len; ++i) // 0 から n-1 まで昇順#define repi(i, s, t) for(int i = int(s), i##_end = int(t); i <= i##_end; ++i) // s から t まで昇順#define repir(i, s, t) for(int i = int(s), i##_end = int(t); i >= i##_end; --i) // s から t まで降順#define repe(v, a) for(const auto& v : (a)) // a の全要素(変更不可能)#define repea(v, a) for(auto& v : (a)) // a の全要素(変更可能)#define repb(set, d) for(int set = 0; set < (1 << int(d)); ++set) // d ビット全探索(昇順)#define repp(a) sort(all(a)); for(bool a##_perm = true; a##_perm; a##_perm = next_permutation(all(a))) // a の順列全て(昇順)#define smod(n, m) ((((n) % (m)) + (m)) % (m)) // 非負mod#define uniq(a) {sort(all(a)); (a).erase(unique(all(a)), (a).end());} // 重複除去#define EXIT(a) {cout << (a) << endl; exit(0);} // 強制終了// 汎用関数の定義template <class T> inline ll pow(T n, int k) { ll v = 1; rep(i, k) v *= n; return v; }template <class T> inline bool chmax(T& M, const T& x) { if (M < x) { M = x; return true; } return false; } // 最大値を更新(更新されたら trueを返す)template <class T> inline bool chmin(T& m, const T& x) { if (m > x) { m = x; return true; } return false; } // 最小値を更新(更新されたら trueを返す)template <class T> inline T get(T set, int i) { return (set >> i) & T(1); }// 演算子オーバーロードtemplate <class T, class U> inline istream& operator>>(istream& is, pair<T, U>& p) { is >> p.first >> p.second; return is; }template <class T> inline istream& operator>>(istream& is, vector<T>& v) { repea(x, v) is >> x; return is; }template <class T> inline vector<T>& operator--(vector<T>& v) { repea(x, v) --x; return v; }template <class T> inline vector<T>& operator++(vector<T>& v) { repea(x, v) ++x; return v; }#endif // 折りたたみ用#if __has_include(<atcoder/all>)#include <atcoder/all>using namespace atcoder;#ifdef _MSC_VER#include "localACL.hpp"#endifusing mint = modint1000000007;//using mint = modint998244353;//using mint = modint; // mint::set_mod(m);namespace atcoder {inline istream& operator>>(istream& is, mint& x) { ll x_; is >> x_; x = x_; return is; }inline ostream& operator<<(ostream& os, const mint& x) { os << x.val(); return os; }}using vm = vector<mint>; using vvm = vector<vm>; using vvvm = vector<vvm>;#endif#ifdef _MSC_VER // 手元環境(Visual Studio)#include "local.hpp"#else // 提出用(gcc)inline int popcount(int n) { return __builtin_popcount(n); }inline int popcount(ll n) { return __builtin_popcountll(n); }inline int lsb(int n) { return n != 0 ? __builtin_ctz(n) : -1; }inline int lsb(ll n) { return n != 0 ? __builtin_ctzll(n) : -1; }inline int msb(int n) { return n != 0 ? (31 - __builtin_clz(n)) : -1; }inline int msb(ll n) { return n != 0 ? (63 - __builtin_clzll(n)) : -1; }#define gcd __gcd#define dump(...)#define dumpel(v)#define dump_list(v)#define dump_mat(v)#define input_from_file(f)#define output_to_file(f)#define Assert(b) { if (!(b)) while (1) cout << "OLE"; }#endifvoid WA() {int n; double p, q;cin >> n >> p >> q;double res = 0;for (int k = 1; k <= 23; k += 2) {repb(set, k) {int x = 0, d = 1;rep(i, k) {if (get(set, i)) {d *= -1;}else {x += d;}if (x == n || (i != k - 1 && x == 0)) {x = n;break;}}double prob = 1;rep(i, k) prob *= get(set, i) ? p : q;if (k < 23) {res += prob * (x == 0);}else {res += prob * (n - x) / n;}}}cout << res << endl;}//【行列】/** Matrix<T>(int n, int m) : O(n m)* n×m 零行列で初期化する.** Matrix<T>(int n) : O(n^2)* n×n 単位行列で初期化する.** Matrix<T>(vvT a) : O(n m)* 二次元配列 a[0..n)[0..m) で初期化する.** bool empty() : O(1)* 行列が空かを返す.** A + B : O(n m)* n×m 行列 A, B の和を返す.+= も使用可.** A - B : O(n m)* n×m 行列 A, B の差を返す.-= も使用可.** c * A / A * c : O(n m)* n×m 行列 A とスカラー c のスカラー積を返す.*= も使用可.** A * x : O(n m)* n×m 行列 A と n 次元列ベクトル x の積を返す.** x * A : O(n m)* m 次元行ベクトル x と n×m 行列 A の積を返す.** A * B : O(n m l)* n×m 行列 A と m×l 行列 B の積を返す.** Mat pow(ll d) : O(n^3 log d)* 自身を d 乗した行列を返す.*/template <class T>struct Matrix {int n, m; // 行列のサイズ(n 行 m 列)vector<vector<T>> v; // 行列の成分// n×m 零行列で初期化する.Matrix(int n, int m) : n(n), m(m), v(n, vector<T>(m)) {}// n×n 単位行列で初期化する.Matrix(int n) : n(n), m(n), v(n, vector<T>(n)) { rep(i, n) v[i][i] = T(1); }// 二次元配列 a[0..n)[0..m) で初期化する.Matrix(const vector<vector<T>>& a) : n(sz(a)), m(sz(a[0])), v(a) {}Matrix() : n(0), m(0) {}// 代入Matrix(const Matrix&) = default;Matrix& operator=(const Matrix&) = default;// アクセスinline vector<T> const& operator[](int i) const { return v[i]; }inline vector<T>& operator[](int i) {// verify : https://judge.yosupo.jp/problem/matrix_product// inline を付けて [] でアクセスするとなぜか v[] への直接アクセスより速くなった.return v[i];}// 入力friend istream& operator>>(istream& is, Matrix& a) {rep(i, a.n) rep(j, a.m) is >> a.v[i][j];return is;}// 空かbool empty() { return min(n, m) == 0; }// 比較bool operator==(const Matrix& b) const { return n == b.n && m == b.m && v == b.v; }bool operator!=(const Matrix& b) const { return !(*this == b); }// 加算,減算,スカラー倍Matrix& operator+=(const Matrix& b) {rep(i, n) rep(j, m) v[i][j] += b[i][j];return *this;}Matrix& operator-=(const Matrix& b) {rep(i, n) rep(j, m) v[i][j] -= b[i][j];return *this;}Matrix& operator*=(const T& c) {rep(i, n) rep(j, m) v[i][j] *= c;return *this;}Matrix operator+(const Matrix& b) const { return Matrix(*this) += b; }Matrix operator-(const Matrix& b) const { return Matrix(*this) -= b; }Matrix operator*(const T& c) const { return Matrix(*this) *= c; }friend Matrix operator*(const T& c, const Matrix<T>& a) { return a * c; }Matrix operator-() const { return Matrix(*this) *= T(-1); }// 行列ベクトル積 : O(m n)vector<T> operator*(const vector<T>& x) const {vector<T> y(n);rep(i, n) rep(j, m) y[i] += v[i][j] * x[j];return y;}// ベクトル行列積 : O(m n)friend vector<T> operator*(const vector<T>& x, const Matrix& a) {vector<T> y(a.m);rep(i, a.n) rep(j, a.m) y[j] += x[i] * a[i][j];return y;}// 積:O(n^3)Matrix operator*(const Matrix& b) const {// verify : https://judge.yosupo.jp/problem/matrix_productMatrix res(n, b.m);rep(i, res.n) rep(j, res.m) rep(k, m) res[i][j] += v[i][k] * b[k][j];return res;}Matrix& operator*=(const Matrix& b) { *this = *this * b; return *this; }// 累乗:O(n^3 log d)Matrix pow(ll d) const {Matrix res(n), pow2 = *this;while (d > 0) {if (d & 1) res *= pow2;pow2 *= pow2;d /= 2;}return res;}#ifdef _MSC_VERfriend ostream& operator<<(ostream& os, const Matrix& a) {rep(i, a.n) {os << "[";rep(j, a.m) os << a[i][j] << " ]"[j == a.m - 1];if (i < a.n - 1) os << "\n";}return os;}#endif};//【線形方程式】O(n m min(n, m))/** 与えられた n×m 行列 A と n 次元ベクトル b に対し,* 線形方程式 A x = b の特殊解 x0(m 次元ベクトル)を返す(なければ空リスト)* また同次形 A x = 0 の解空間の基底(m 次元ベクトル)のリストを xs に格納する.*/template <class T>vector<T> gauss_jordan_elimination(const Matrix<T>& A, const vector<T>& b, vector<vector<T>>* xs = nullptr) {// verify : https://judge.yosupo.jp/problem/system_of_linear_equationsint n = A.n, m = A.m;// v : 拡大係数行列 (A | b)vector<vector<T>> v(n, vector<T>(m + 1));rep(i, n) rep(j, m) v[i][j] = A[i][j];rep(i, n) v[i][m] = b[i];// pivots[i] : 第 i 行のピボットが第何列にあるかvi pivots;// 直前に見つけたピボットの位置int pi = -1, pj = -1;// 注目位置を v[i][j] とする.int i = 0, j = 0;while (i < n && j <= m) {// 同じ列の下方の行から非 0 成分を見つける.int i2 = i;while (i2 < n && v[i2][j] == 0) i2++;// 見つからなかったら注目位置を右に移す.if (i2 == n) {j++;continue;}// 見つかったら第 i 行とその行を入れ替える.pi = i; pj = j;if (i != i2) swap(v[i], v[i2]);// v[i][j] をピボットに選択する.pivots.push_back(j);// v[i][j] が 1 になるよう第 i 行全体を v[i][j] で割る.T vij_inv = T(1) / v[i][j];repi(j2, j, m) v[i][j2] *= vij_inv;// 第 i 行以外の第 j 列の成分が全て 0 になるよう第 i 行を定数倍して減じる.rep(i2, n) {if (v[i2][j] == T(0) || i2 == i) continue;T mul = v[i2][j];repi(j2, j, m) v[i2][j2] -= v[i][j2] * mul;}// 注目位置を右下に移す.i++; j++;}// 最後に見つかったピボットの位置が第 m 列ならば解なし.if (pivots.back() == m) return vector<T>();// A x = b の特殊解 x0 の構成(任意定数は全て 0 にする)vector<T> x0(m);int rnk = sz(pivots);rep(i, rnk) x0[pivots[i]] = v[i][m];// 同次形 A x = 0 の一般解 {x} の基底の構成(任意定数を 1-hot にする)if (xs != nullptr) {xs->clear();int i = 0;rep(j, m) {if (i < rnk && j == pivots[i]) {i++;continue;}vector<T> x(m, T(0));x[j] = 1;rep(i2, i) x[pivots[i2]] = -v[i2][j];xs->emplace_back(move(x));}}return x0;}//【ランダムウォーク】/** Random_walk<T>(int n) : O(1)* n 頂点 0 辺のグラフで初期化する.** add_edge(int s, int t, T prob) : O(1)* 有向辺 s→t を,選択確率 prob で追加する.* 制約:任意の s について Σs→t p[s][t] = 1** vT arrive_probability_to(int GL) : O(n^3)(バグってる)* 各頂点から出発し GL に到着する確率のリストを返す.** vT expected_turn_to(int GL) : O(n^3)* 各頂点から出発し GL に初めて到着するまでのターン数の期待値のリストを返す.* 制約:どの頂点からも GL に到達可能** vT stationary_distribution() : O(n^3)* 定常分布を返す.* 制約:どの頂点からどの頂点へも移動可能** 利用:【行列】,【線形方程式】*/template <class T>class Random_walk {int n;// 推移確率行列(p[i][j] : i から j に移動する確率)vector<vector<T>> p;public:// n 頂点 0 辺のグラフで初期化する.Random_walk(int n) : n(n), p(n, vector<T>(n)) {}Random_walk() : n(0) {}// 有向辺 s→t を,選択確率 prob で追加する.void add_edge(int s, int t, T prob) {p[s][t] += prob;}// 各頂点から出発し GL に到着する確率のリストを返す.vector<T> arrive_probability_to(int GL) {//【方法】// s から GL に到着する確率を x[s] とすると,線形方程式// x[s] = Σs→t p[s][t] x[t] (s ≠ GL)// x[GL] = 1// を得る.これを整理すると// (1 - p[s][s])x[s] - Σs→t,t≠s p[s][t] x[t] = 0// x[GL] = 1// となる.Matrix<T> mat(n); vector<T> vec(n);rep(i, n) rep(j, n) if (i != GL) mat[i][j] -= p[i][j];vec[GL] = 1;return gauss_jordan_elimination(mat, vec);}// 各頂点から出発し GL に初めて到着するまでのターン数の期待値のリストを返す.vector<T> expected_turn_to(int GL) {//【方法】// s→GL にかかるターン数の期待値を e[s] とすると,線形方程式// e[s] = 1 + Σs→t p[s][t] e[t] (s ≠ GL)// e[GL] = 0// を得る.これを整理すると// (1 - p[s][s])e[s] - Σs→t,t≠s p[s][t] e[t] = 1// e[GL] = 0// となる.Matrix<T> mat(n); vector<T> vec(n, 1);rep(i, n) rep(j, n) if (i != GL) mat[i][j] -= p[i][j];vec[GL] = 0;return gauss_jordan_elimination(mat, vec);}// 定常分布を返す.vector<T> stationary_distribution() {//【方法】// 定常分布を π[0..n) とすると,線形方程式// π[t] = Σs→t p[s][t] π[s]// Σπ[0..n) = 1// を得る.これを整理すると// (1 - p[t][t])π[t] - Σs→t,t≠s p[s][t] π[s] = 0// Σπ[0..n) = 1// となる.Matrix<T> mat(n); vector<T> vec(n);rep(i, n - 1) rep(j, n) mat[i][j] -= p[j][i];rep(j, n) mat[n - 1][j] = 1;vec[n - 1] = 1;return gauss_jordan_elimination(mat, vec);}#ifdef _MSC_VERfriend ostream& operator<<(ostream& os, const Random_walk& rw) {rep(i, rw.n) {rep(j, rw.n) os << rw.p[i][j] << " ";os << endl;}return os;}#endif};void check_Random_walk() {Random_walk<double> RW(4);RW.add_edge(0, 1, 0.5);RW.add_edge(0, 0, 0.5);RW.add_edge(1, 2, 1);RW.add_edge(2, 3, 1);RW.add_edge(3, 0, 1);dump(RW.expected_turn_to(3)); // 4 2 1 0dump(RW.arrive_probability_to(3)); // 1 1 1 1dump(RW.stationary_distribution()); // 0.4 0.2 0.2 0.2exit(0);}int main() {// input_from_file("input.txt");// output_to_file("output.txt");// check_Random_walk();int n; double p, q;cin >> n >> p >> q;int DEATH = 2 * (n + 1);Random_walk<double> g(DEATH + 1);repi(i, 0, n - 1) {g.add_edge(2 * i + 0, 2 * i + 1, p);g.add_edge(2 * i + 0, 2 * (i + 1) + 0, q);g.add_edge(2 * i + 0, DEATH, 1 - p - q);}g.add_edge(2 * n + 0, DEATH, 1);repi(i, 1, n) {g.add_edge(2 * i + 1, 2 * i + 0, p);g.add_edge(2 * i + 1, 2 * (i - 1) + 1, q);g.add_edge(2 * i + 1, DEATH, 1 - p - q);}g.add_edge(2 * 0 + 1, DEATH, 1);g.add_edge(DEATH, DEATH, 1);// dump(g);auto res = g.arrive_probability_to(2 * 0 + 1);dump(res);cout << res[2 * 0 + 0] << endl;}