#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2171" #define ERROR 1e-8 #include <iomanip> #include <iostream> #include <vector> #include "Mylib/Graph/ShortestPath/warshall_floyd_for_matrix_graph.cpp" #include "Mylib/IO/input_vector.cpp" #include "Mylib/LinearAlgebra/simultaneous_linear_equations_float.cpp" namespace hl = haar_lib; int main() { std::cin.tie(0); std::ios::sync_with_stdio(false); int n, s, t; while (std::cin >> n >> s >> t, n) { --s, --t; auto q = hl::input_vector<int>(n); auto g = hl::input_vector<int>(n, n); auto dist = hl::warshall_floyd_for_matrix<int, 0>(g); if (not dist[s][t]) { std::cout << "impossible" << std::endl; continue; } std::vector<std::vector<double>> a(n, std::vector<double>(n)); std::vector<double> b(n); for (int i = 0; i < n; ++i) { if (i == t) { a[i][i] = 1; b[i] = 0; continue; } if (q[i]) { int d = *dist[i][t]; int k = 0; for (int j = 0; j < n; ++j) { if (g[i][j] > 0 and g[i][j] + *dist[j][t] == d) { ++k; } } for (int j = 0; j < n; ++j) { if (g[i][j] > 0 and g[i][j] + *dist[j][t] == d) { a[i][j] = -1; b[i] += g[i][j]; } } a[i][i] += k; } else { int k = 0; for (int j = 0; j < n; ++j) { if (g[i][j] > 0) { ++k; } } for (int j = 0; j < n; ++j) { if (g[i][j] > 0) { a[i][j] = -1; b[i] += g[i][j]; } } a[i][i] += k; } } auto res = hl::float_simultaneous_linear_equations(a, b, ERROR); double ans = (*res).solution[s]; std::cout << std::setprecision(12) << std::fixed << ans << std::endl; } return 0; }
#line 1 "test/aoj/2171/main.test.cpp" #define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2171" #define ERROR 1e-8 #include <iomanip> #include <iostream> #include <vector> #line 2 "Mylib/Graph/ShortestPath/warshall_floyd_for_matrix_graph.cpp" #include <optional> #line 4 "Mylib/Graph/ShortestPath/warshall_floyd_for_matrix_graph.cpp" namespace haar_lib { namespace warshall_floyd_for_matrix_impl { template <typename T> struct result { std::vector<std::vector<std::optional<T>>> dist; bool has_negative_cycle; const auto& operator[](int i) const { return dist[i]; } }; } // namespace warshall_floyd_for_matrix_impl template <typename T, T INVALID> auto warshall_floyd_for_matrix(const std::vector<std::vector<T>>& g) { const int n = g.size(); auto dist = std::vector(n, std::vector<std::optional<T>>(n)); for (int i = 0; i < n; ++i) dist[i][i] = 0; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { if (g[i][j] != INVALID) { dist[i][j] = g[i][j]; } } } for (int k = 0; k < n; ++k) { for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { if (dist[i][k] and dist[k][j]) { if (not dist[i][j]) { dist[i][j] = *dist[i][k] + *dist[k][j]; } else { dist[i][j] = std::min(*dist[i][j], *dist[i][k] + *dist[k][j]); } } } } } bool has_negative_cycle = false; for (int i = 0; i < n; ++i) if (*dist[i][i] < 0) has_negative_cycle = true; return warshall_floyd_for_matrix_impl::result<T>{dist, has_negative_cycle}; } } // namespace haar_lib #line 4 "Mylib/IO/input_vector.cpp" namespace haar_lib { template <typename T> std::vector<T> input_vector(int N) { std::vector<T> ret(N); for (int i = 0; i < N; ++i) std::cin >> ret[i]; return ret; } template <typename T> std::vector<std::vector<T>> input_vector(int N, int M) { std::vector<std::vector<T>> ret(N); for (int i = 0; i < N; ++i) ret[i] = input_vector<T>(M); return ret; } } // namespace haar_lib #line 3 "Mylib/LinearAlgebra/simultaneous_linear_equations_float.cpp" #include <utility> #line 5 "Mylib/LinearAlgebra/simultaneous_linear_equations_float.cpp" namespace haar_lib { namespace float_simultaneous_linear_equations_impl { template <typename T> struct result { int rank, dim; std::vector<T> solution; }; } // namespace float_simultaneous_linear_equations_impl template <typename T> auto float_simultaneous_linear_equations(std::vector<std::vector<T>> a, std::vector<T> b, T eps) { using result = float_simultaneous_linear_equations_impl::result<T>; std::optional<result> ret; const int n = a.size(), m = a[0].size(); int rank = 0; for (int j = 0; j < m; ++j) { int pivot = -1; double M = eps; for (int i = rank; i < n; ++i) { if (std::abs(a[i][j]) > M) { M = std::abs(a[i][j]); pivot = i; } } if (pivot == -1) continue; std::swap(a[pivot], a[rank]); std::swap(b[pivot], b[rank]); { double d = a[rank][j]; for (int k = 0; k < m; ++k) a[rank][k] /= d; b[rank] /= d; } for (int i = 0; i < n; ++i) { if (i == rank or std::abs(a[i][j]) <= eps) continue; double d = a[i][j]; for (int k = 0; k < m; ++k) { a[i][k] -= a[rank][k] * d; } b[i] -= b[rank] * d; } ++rank; } for (int i = rank; i < n; ++i) { if (std::abs(b[i]) > eps) { return ret; } } const int dim = m - rank; std::vector<T> solution(m); for (int i = 0; i < rank; ++i) solution[i] = b[i]; ret = result({rank, dim, solution}); return ret; } } // namespace haar_lib #line 10 "test/aoj/2171/main.test.cpp" namespace hl = haar_lib; int main() { std::cin.tie(0); std::ios::sync_with_stdio(false); int n, s, t; while (std::cin >> n >> s >> t, n) { --s, --t; auto q = hl::input_vector<int>(n); auto g = hl::input_vector<int>(n, n); auto dist = hl::warshall_floyd_for_matrix<int, 0>(g); if (not dist[s][t]) { std::cout << "impossible" << std::endl; continue; } std::vector<std::vector<double>> a(n, std::vector<double>(n)); std::vector<double> b(n); for (int i = 0; i < n; ++i) { if (i == t) { a[i][i] = 1; b[i] = 0; continue; } if (q[i]) { int d = *dist[i][t]; int k = 0; for (int j = 0; j < n; ++j) { if (g[i][j] > 0 and g[i][j] + *dist[j][t] == d) { ++k; } } for (int j = 0; j < n; ++j) { if (g[i][j] > 0 and g[i][j] + *dist[j][t] == d) { a[i][j] = -1; b[i] += g[i][j]; } } a[i][i] += k; } else { int k = 0; for (int j = 0; j < n; ++j) { if (g[i][j] > 0) { ++k; } } for (int j = 0; j < n; ++j) { if (g[i][j] > 0) { a[i][j] = -1; b[i] += g[i][j]; } } a[i][i] += k; } } auto res = hl::float_simultaneous_linear_equations(a, b, ERROR); double ans = (*res).solution[s]; std::cout << std::setprecision(12) << std::fixed << ans << std::endl; } return 0; }