This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub Luzhiled/comp-geometry
// verification-helper: PROBLEM https://onlinejudge.u-aizu.ac.jp/problems/2950 // verification-helper: ERROR 1e-4 #include <random> #include <chrono> namespace lib { using namespace std; // Rolling Hash {{{ struct RollingHash { static const uint64_t mod = (1ull << 61ull) - 1; using uint128_t = __uint128_t; const uint64_t base; vector< uint64_t > power; static inline uint64_t add(uint64_t a, uint64_t b) { if((a += b) >= mod) a -= mod; return a; } static inline uint64_t mul(uint64_t a, uint64_t b) { uint128_t c = (uint128_t) a * b; return add(c >> 61, c & mod); } static inline uint64_t generate_base() { mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count()); uniform_int_distribution< uint64_t > rand(1, RollingHash::mod - 1); return rand(mt); } inline void expand(size_t sz) { if(power.size() < sz + 1) { int pre_sz = (int) power.size(); power.resize(sz + 1); for(int i = pre_sz - 1; i < (int)sz; i++) { power[i + 1] = mul(power[i], base); } } } explicit RollingHash(uint64_t base = generate_base()) : base(base), power{1} {} vector< uint64_t > build(const string &s) const { int sz = s.size(); vector< uint64_t > hashed(sz + 1); for(int i = 0; i < sz; i++) { hashed[i + 1] = add(mul(hashed[i], base), s[i]); } return hashed; } template< typename T > vector< uint64_t > build(const vector< T > &s) const { int sz = s.size(); vector< uint64_t > hashed(sz + 1); for(int i = 0; i < sz; i++) { hashed[i + 1] = add(mul(hashed[i], base), s[i]); } return hashed; } uint64_t query(const vector< uint64_t > &s, int l, int r) { expand(r - l); return add(s[r], mod - mul(s[l], power[r - l])); } uint64_t combine(uint64_t h1, uint64_t h2, size_t h2len) { expand(h2len); return add(mul(h1, power[h2len]), h2); } int lcp(const vector< uint64_t > &a, int l1, int r1, const vector< uint64_t > &b, int l2, int r2) { int len = min(r1 - l1, r2 - l2); int low = 0, high = len + 1; while(high - low > 1) { int mid = (low + high) / 2; if(query(a, l1, l1 + mid) == query(b, l2, l2 + mid)) low = mid; else high = mid; } return low; } }; // }}} } #include "src/real-geometry/position/point-polygon-positional-relationships.hpp" #include "src/real-geometry/point-cloud/convex-hull-with-index.hpp" #include "src/real-geometry/position/intersect-ss.hpp" #include "src/real-geometry/operation/cross-product.hpp" #include "src/real-geometry/utility/sign.hpp" #include <iostream> #include <complex> #include <set> #include <unordered_set> #include <queue> #include <utility> #include <vector> #include <functional> template< typename T > using vector = std::vector<T>; void solve(int n, int k) { lib::RollingHash roll; using u64 = long long; using std::pair; using R = long double; using points = geometry::points<R>; using polygon = geometry::polygon<R>; using segment = geometry::segment<R>; points pts(n); for (auto &pt : pts) std::cin >> pt; std::unordered_set< u64 > used; using T = pair<double, int>; std::priority_queue< T, vector<T>, std::greater<T> > pq; // TODO: #37 auto calc_perimeter = [&](const vector< int > &vs) { using std::abs; double len = abs(pts[vs.front()] - pts[vs.back()]); for (int i = 1; i < (int)vs.size(); i++) { len += abs(pts[vs[i]] - pts[vs[i - 1]]); } return len; }; auto calc_hash = [&](const vector< int > &vs) { auto rh = roll.build(vs); return roll.query(rh, 0, vs.size()); }; auto insert_ptsi = [&](vector< int > vs, int i, int j) { vs.insert(vs.begin() + j + 1, i); return vs; }; vector< vector<int> > vss; { auto ds = geometry::convex_hull_with_index(pts); vector< int > vs; for (auto &v : ds.second) vs.emplace_back(v); u64 hash = calc_hash(vs); double len = calc_perimeter(vs); pq.emplace(len, vss.size()); vss.emplace_back(vs); used.emplace(hash); } for (int qi = 1; qi < k and not pq.empty(); qi++) { auto [d, idx] = pq.top(); pq.pop(); auto as = vss[idx]; int m = as.size(); std::set< int > st(as.begin(), as.end()); for (int i = 0; i < n; i++) { if (st.count(i)) continue; for (int j = 0; j < m; j++) { auto vs = insert_ptsi(as, i, j); u64 hash = calc_hash(vs); if (used.count(hash)) { continue; } used.emplace(hash); // TODO: #35 segment s1(pts[i], pts[as[j]]); segment s2(pts[i], pts[as[(j + 1) % m]]); int cnt = 0; for (int k = 0; k < m; k++) { segment s(pts[as[k]], pts[as[(k + 1) % m]]); if (intersect_ss(s, s1)) cnt++; if (intersect_ss(s, s2)) cnt++; } if (cnt != 4) continue; polygon poly; bool f = false; for (auto i : vs) poly.emplace_back(pts[i]); for (auto &p : pts) if (point_polygon_positional_relationships(p, poly) == 0) f = true; if (f) continue; double len = calc_perimeter(vs); pq.emplace(len, vss.size()); vss.emplace_back(vs); } } } if (pq.empty()) { std::cout << -1 << std::endl; return; } std::cout << pq.top().first << std::endl; } signed main() { int n, k; while (std::cin >> n >> k, n) { solve(n, k); } }
#line 1 "test/aoj/icpc/2950.test.cpp" // verification-helper: PROBLEM https://onlinejudge.u-aizu.ac.jp/problems/2950 // verification-helper: ERROR 1e-4 #include <random> #include <chrono> namespace lib { using namespace std; // Rolling Hash {{{ struct RollingHash { static const uint64_t mod = (1ull << 61ull) - 1; using uint128_t = __uint128_t; const uint64_t base; vector< uint64_t > power; static inline uint64_t add(uint64_t a, uint64_t b) { if((a += b) >= mod) a -= mod; return a; } static inline uint64_t mul(uint64_t a, uint64_t b) { uint128_t c = (uint128_t) a * b; return add(c >> 61, c & mod); } static inline uint64_t generate_base() { mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count()); uniform_int_distribution< uint64_t > rand(1, RollingHash::mod - 1); return rand(mt); } inline void expand(size_t sz) { if(power.size() < sz + 1) { int pre_sz = (int) power.size(); power.resize(sz + 1); for(int i = pre_sz - 1; i < (int)sz; i++) { power[i + 1] = mul(power[i], base); } } } explicit RollingHash(uint64_t base = generate_base()) : base(base), power{1} {} vector< uint64_t > build(const string &s) const { int sz = s.size(); vector< uint64_t > hashed(sz + 1); for(int i = 0; i < sz; i++) { hashed[i + 1] = add(mul(hashed[i], base), s[i]); } return hashed; } template< typename T > vector< uint64_t > build(const vector< T > &s) const { int sz = s.size(); vector< uint64_t > hashed(sz + 1); for(int i = 0; i < sz; i++) { hashed[i + 1] = add(mul(hashed[i], base), s[i]); } return hashed; } uint64_t query(const vector< uint64_t > &s, int l, int r) { expand(r - l); return add(s[r], mod - mul(s[l], power[r - l])); } uint64_t combine(uint64_t h1, uint64_t h2, size_t h2len) { expand(h2len); return add(mul(h1, power[h2len]), h2); } int lcp(const vector< uint64_t > &a, int l1, int r1, const vector< uint64_t > &b, int l2, int r2) { int len = min(r1 - l1, r2 - l2); int low = 0, high = len + 1; while(high - low > 1) { int mid = (low + high) / 2; if(query(a, l1, l1 + mid) == query(b, l2, l2 + mid)) low = mid; else high = mid; } return low; } }; // }}} } #line 2 "src/real-geometry/position/point-polygon-positional-relationships.hpp" #line 2 "src/real-geometry/class/point.hpp" #line 2 "src/real-geometry/class/vector.hpp" #include <complex> #include <iostream> namespace geometry { template< typename R > class vec2d : public std::complex< R > { using complex = std::complex< R >; public: using complex::complex; vec2d(const complex &c): complex::complex(c) {} const R x() const { return this->real(); } const R y() const { return this->imag(); } friend vec2d operator*(const vec2d &v, const R &k) { return vec2d(v.x() * k, v.y() * k); } friend vec2d operator*(const R &k, const vec2d &v) { return vec2d(v.x() * k, v.y() * k); } friend std::istream &operator>>(std::istream &is, vec2d &v) { R x, y; is >> x >> y; v = vec2d(x, y); return is; } }; } #line 4 "src/real-geometry/class/point.hpp" #include <vector> namespace geometry { template< typename R > using point = vec2d<R>; template< typename R > using points = std::vector< point< R > >; } #line 2 "src/real-geometry/class/polygon.hpp" #line 4 "src/real-geometry/class/polygon.hpp" #line 6 "src/real-geometry/class/polygon.hpp" namespace geometry { template< typename R > using polygon = std::vector< point<R> >; template< typename R > using polygons = std::vector< polygon<R> >; } #line 2 "src/real-geometry/numbers/posision-of-point-polygon.hpp" namespace geometry::number::point_polygon_positional_relationships { constexpr int OUT = 0; constexpr int ON_EDGE = 1; constexpr int IN = 2; } #line 2 "src/real-geometry/operation/cross-product.hpp" #line 4 "src/real-geometry/operation/cross-product.hpp" namespace geometry { template< typename R > R cross_product(const vec2d<R> &a, const vec2d<R> &b) { return a.x() * b.y() - a.y() * b.x(); } } #line 2 "src/real-geometry/operation/inner-product.hpp" #line 4 "src/real-geometry/operation/inner-product.hpp" namespace geometry { template< typename R > R inner_product(const vec2d<R> &a, const vec2d<R> &b) { return a.x() * b.x() + a.y() * b.y(); } } #line 2 "src/real-geometry/common/size-alias.hpp" #include <cstddef> namespace geometry { using isize = std::ptrdiff_t; using usize = std::size_t; } #line 2 "src/real-geometry/utility/next-idx.hpp" #line 4 "src/real-geometry/utility/next-idx.hpp" namespace geometry { inline usize next_idx(usize idx, usize size) { return idx + 1 == size ? 0 : idx + 1; } } #line 2 "src/real-geometry/utility/sign.hpp" #line 2 "src/real-geometry/common/const/eps.hpp" #line 2 "src/real-geometry/common/float-alias.hpp" namespace geometry { using f80 = long double; using f64 = double; } #line 4 "src/real-geometry/common/const/eps.hpp" namespace geometry { inline static f80 &eps() { static f80 EPS = 1e-10; return EPS; } void set_eps(f80 EPS) { eps() = EPS; } } #line 2 "src/real-geometry/numbers/sign.hpp" #line 2 "src/real-geometry/common/int-alias.hpp" namespace geometry { using i32 = int; using i64 = long long; } #line 4 "src/real-geometry/numbers/sign.hpp" namespace geometry::number::sign { constexpr i32 PLUS = +1; constexpr i32 ZERO = 0; constexpr i32 MINUS = -1; } #line 5 "src/real-geometry/utility/sign.hpp" namespace geometry { using namespace geometry::number::sign; template< typename R > inline int sign(R r) { if (r < -eps()) return MINUS; if (r > +eps()) return PLUS; return ZERO; } } #line 11 "src/real-geometry/position/point-polygon-positional-relationships.hpp" #include <algorithm> namespace geometry { // O(N) template< typename R > int point_polygon_positional_relationships(const point<R> &p, const polygon<R> &poly) { using namespace number::point_polygon_positional_relationships; usize n = poly.size(); bool in = false; for (usize i = 0; i < n; i++) { usize j = next_idx(i, n); point<R> a = poly[i] - p, b = poly[j] - p; if (a.y() > b.y()) std::swap(a, b); if (a.y() <= 0 and 0 < b.y() and cross_product(a, b) < 0) { in = not in; } if (sign(cross_product(a, b)) == 0 and sign(inner_product(a, b)) <= 0) { return ON_EDGE; } } return in ? IN : OUT; } } #line 2 "src/real-geometry/point-cloud/convex-hull-with-index.hpp" #line 2 "src/real-geometry/compare/compare-x.hpp" #line 2 "src/real-geometry/utility/equals/real-number.hpp" #line 4 "src/real-geometry/utility/equals/real-number.hpp" namespace geometry { template< typename R > bool equals(R a, R b) { return sign(a - b) == 0; } } #line 5 "src/real-geometry/compare/compare-x.hpp" namespace geometry { template< typename R > bool compare_x(const point<R> &a, const point<R> &b) { return not equals(a.x(), b.x()) ? a.x() < b.x() : a.y() < b.y(); } } #line 8 "src/real-geometry/point-cloud/convex-hull-with-index.hpp" #line 10 "src/real-geometry/point-cloud/convex-hull-with-index.hpp" #include <numeric> #include <tuple> #include <utility> #line 14 "src/real-geometry/point-cloud/convex-hull-with-index.hpp" namespace geometry { template< typename R > std::pair< polygon<R>, std::vector< usize > > convex_hull_with_index(const points<R> &pts) { usize n = pts.size(); if (n <= 2) { std::vector< usize > idxs(n); std::iota(idxs.begin(), idxs.end(), 0); return {pts, idxs}; } std::vector< std::pair< point<R>, usize > > ps; ps.reserve(n); for (usize i = 0; i < n; i++) { ps.emplace_back(pts[i], i); } auto cmp = [](const std::pair<point<R>, usize> &a, const std::pair<point<R>, usize> &b) { return compare_x(a.first, b.first); }; std::sort(ps.begin(), ps.end(), cmp); std::vector< usize > idxs(2 * n); polygon<R> poly(2 * n); usize k = 0, i = 0; auto check = [&](usize i) { return sign(cross_product<R>(poly[k - 1] - poly[k - 2], ps[i].first - poly[k - 1])) == -1; }; while (i < n) { while (k >= 2 and check(i)) k--; std::tie(poly[k], idxs[k]) = ps[i]; k++; i++; } i = n - 2; usize t = k + 1; while (true) { while (k >= t and check(i)) k--; std::tie(poly[k], idxs[k]) = ps[i]; k++; if (not i) break; i--; } poly.resize(k - 1); idxs.resize(k - 1); return {poly, idxs}; } } #line 2 "src/real-geometry/position/intersect-ss.hpp" #line 2 "src/real-geometry/class/segment.hpp" #line 2 "src/real-geometry/utility/equals/vector.hpp" #line 5 "src/real-geometry/utility/equals/vector.hpp" namespace geometry { template< typename R > bool equals(const vec2d<R> &a, const vec2d<R> &b) { return equals(a.x(), b.x()) and equals(a.y(), b.y()); } } #line 5 "src/real-geometry/class/segment.hpp" #include <cassert> #line 8 "src/real-geometry/class/segment.hpp" namespace geometry { template< typename R > class segment { public: point<R> a, b; segment() = default; segment(point<R> a, point<R> b) : a(a), b(b) { assert(not equals(a, b)); } }; template< typename R > using segments = std::vector< segment<R> >; } #line 2 "src/real-geometry/operation/ccw.hpp" #line 2 "src/real-geometry/numbers/ccw.hpp" namespace geometry::number::ccw { constexpr int COUNTER_CLOCKWISE = +1; constexpr int CLOCKWISE = -1; constexpr int ONLINE_BACK = +2; // c-a-b constexpr int ONLINE_FRONT = -2; // a-b-c constexpr int ON_SEGMENT = 0; // a-c-b } #line 8 "src/real-geometry/operation/ccw.hpp" namespace geometry { using namespace geometry::number::ccw; template< typename R > int ccw(const point<R> &a, point<R> b, point<R> c) { b = b - a, c = c - a; if (sign(cross_product(b, c)) == +1) return COUNTER_CLOCKWISE; if (sign(cross_product(b, c)) == -1) return CLOCKWISE; if (sign(inner_product(b, c)) == -1) return ONLINE_BACK; if (std::norm(b) < std::norm(c)) return ONLINE_FRONT; return ON_SEGMENT; } } #line 5 "src/real-geometry/position/intersect-ss.hpp" namespace geometry { template< typename R > bool intersect_ss(const segment<R> &s1, const segment<R> &s2) { return ccw(s1.a, s1.b, s2.a) * ccw(s1.a, s1.b, s2.b) <= 0 and ccw(s2.a, s2.b, s1.a) * ccw(s2.a, s2.b, s1.b) <= 0; } } #line 91 "test/aoj/icpc/2950.test.cpp" #line 94 "test/aoj/icpc/2950.test.cpp" #include <set> #include <unordered_set> #include <queue> #line 99 "test/aoj/icpc/2950.test.cpp" #include <functional> template< typename T > using vector = std::vector<T>; void solve(int n, int k) { lib::RollingHash roll; using u64 = long long; using std::pair; using R = long double; using points = geometry::points<R>; using polygon = geometry::polygon<R>; using segment = geometry::segment<R>; points pts(n); for (auto &pt : pts) std::cin >> pt; std::unordered_set< u64 > used; using T = pair<double, int>; std::priority_queue< T, vector<T>, std::greater<T> > pq; // TODO: #37 auto calc_perimeter = [&](const vector< int > &vs) { using std::abs; double len = abs(pts[vs.front()] - pts[vs.back()]); for (int i = 1; i < (int)vs.size(); i++) { len += abs(pts[vs[i]] - pts[vs[i - 1]]); } return len; }; auto calc_hash = [&](const vector< int > &vs) { auto rh = roll.build(vs); return roll.query(rh, 0, vs.size()); }; auto insert_ptsi = [&](vector< int > vs, int i, int j) { vs.insert(vs.begin() + j + 1, i); return vs; }; vector< vector<int> > vss; { auto ds = geometry::convex_hull_with_index(pts); vector< int > vs; for (auto &v : ds.second) vs.emplace_back(v); u64 hash = calc_hash(vs); double len = calc_perimeter(vs); pq.emplace(len, vss.size()); vss.emplace_back(vs); used.emplace(hash); } for (int qi = 1; qi < k and not pq.empty(); qi++) { auto [d, idx] = pq.top(); pq.pop(); auto as = vss[idx]; int m = as.size(); std::set< int > st(as.begin(), as.end()); for (int i = 0; i < n; i++) { if (st.count(i)) continue; for (int j = 0; j < m; j++) { auto vs = insert_ptsi(as, i, j); u64 hash = calc_hash(vs); if (used.count(hash)) { continue; } used.emplace(hash); // TODO: #35 segment s1(pts[i], pts[as[j]]); segment s2(pts[i], pts[as[(j + 1) % m]]); int cnt = 0; for (int k = 0; k < m; k++) { segment s(pts[as[k]], pts[as[(k + 1) % m]]); if (intersect_ss(s, s1)) cnt++; if (intersect_ss(s, s2)) cnt++; } if (cnt != 4) continue; polygon poly; bool f = false; for (auto i : vs) poly.emplace_back(pts[i]); for (auto &p : pts) if (point_polygon_positional_relationships(p, poly) == 0) f = true; if (f) continue; double len = calc_perimeter(vs); pq.emplace(len, vss.size()); vss.emplace_back(vs); } } } if (pq.empty()) { std::cout << -1 << std::endl; return; } std::cout << pq.top().first << std::endl; } signed main() { int n, k; while (std::cin >> n >> k, n) { solve(n, k); } }