#783098 (C++17(clang)) No.5007 Steiner Space Travel

• 問題: No.5007 Steiner Space Travel
• ユーザー: nagiss
• 提出日時: 2022-07-30 16:22:47
• 言語: C++17(clang) (17.0.6 + boost 1.83.0)
• 結果: AC
• 実行時間: 906 ms / 1,000 ms
• コード長: 10,598 bytes
• コンパイル時間: 30,416 ms
• 実行使用メモリ: 3,696 KB
• スコア: 8,241,906
• 最終ジャッジ日時: 2022-07-31 03:12:54
• 合計ジャッジ時間: 31,101 ms
• ジャッジサーバーID （参考情報）: judge15 / judge13

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	903 ms; 3,604 KB
testcase_01	AC	904 ms; 3,660 KB
testcase_02	AC	903 ms; 3,584 KB
testcase_03	AC	903 ms; 3,668 KB
testcase_04	AC	904 ms; 3,596 KB
testcase_05	AC	905 ms; 3,632 KB
testcase_06	AC	903 ms; 3,696 KB
testcase_07	AC	903 ms; 3,620 KB
testcase_08	AC	904 ms; 3,584 KB
testcase_09	AC	902 ms; 3,588 KB
testcase_10	AC	906 ms; 3,664 KB
testcase_11	AC	904 ms; 3,592 KB
testcase_12	AC	903 ms; 3,632 KB
testcase_13	AC	902 ms; 3,616 KB
testcase_14	AC	902 ms; 3,576 KB
testcase_15	AC	904 ms; 3,588 KB
testcase_16	AC	905 ms; 3,592 KB
testcase_17	AC	902 ms; 3,600 KB
testcase_18	AC	904 ms; 3,608 KB
testcase_19	AC	906 ms; 3,520 KB
testcase_20	AC	904 ms; 3,576 KB
testcase_21	AC	903 ms; 3,592 KB
testcase_22	AC	905 ms; 3,588 KB
testcase_23	AC	905 ms; 3,612 KB
testcase_24	AC	905 ms; 3,564 KB
testcase_25	AC	903 ms; 3,588 KB
testcase_26	AC	906 ms; 3,592 KB
testcase_27	AC	902 ms; 3,664 KB
testcase_28	AC	906 ms; 3,592 KB
testcase_29	AC	904 ms; 3,584 KB

ソースコード

#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <climits>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

using namespace std;
using ll = long long;
constexpr double PI = 3.1415926535897932;

template <class T, class S> inline bool chmin(T& m, S q) {
    if (m > q) {
        m = q;
        return true;
    } else
        return false;
}

template <class T, class S> inline bool chmax(T& m, const S q) {
    if (m < q) {
        m = q;
        return true;
    } else
        return false;
}

// 2 次元ベクトル
template <typename T> struct Vec2 {
    /*
    y 軸正は下方向
    x 軸正は右方向
    回転は時計回りが正（y 軸正を上と考えると反時計回りになる）
    */
    using value_type = T;
    T y, x;
    constexpr inline Vec2() = default;
    constexpr Vec2(const T& arg_y, const T& arg_x) : y(arg_y), x(arg_x) {}
    inline Vec2(const Vec2&) = default;            // コピー
    inline Vec2(Vec2&&) = default;                 // ムーブ
    inline Vec2& operator=(const Vec2&) = default; // 代入
    inline Vec2& operator=(Vec2&&) = default;      // ムーブ代入
    template <typename S>
    constexpr inline Vec2(const Vec2<S>& v) : y((T)v.y), x((T)v.x) {}
    inline Vec2 operator+(const Vec2& rhs) const {
        return Vec2(y + rhs.y, x + rhs.x);
    }
    inline Vec2 operator+(const T& rhs) const { return Vec2(y + rhs, x + rhs); }
    inline Vec2 operator-(const Vec2& rhs) const {
        return Vec2(y - rhs.y, x - rhs.x);
    }
    template <typename S> inline Vec2 operator*(const S& rhs) const {
        return Vec2(y * rhs, x * rhs);
    }
    inline Vec2 operator*(const Vec2& rhs) const { // x + yj とみなす
        return Vec2(x * rhs.y + y * rhs.x, x * rhs.x - y * rhs.y);
    }
    template <typename S> inline Vec2 operator/(const S& rhs) const {
        assert(rhs != 0.0);
        return Vec2(y / rhs, x / rhs);
    }
    inline Vec2 operator/(const Vec2& rhs) const { // x + yj とみなす
        return (*this) * rhs.inv();
    }
    inline Vec2& operator+=(const Vec2& rhs) {
        y += rhs.y;
        x += rhs.x;
        return *this;
    }
    inline Vec2& operator-=(const Vec2& rhs) {
        y -= rhs.y;
        x -= rhs.x;
        return *this;
    }
    template <typename S> inline Vec2& operator*=(const S& rhs) const {
        y *= rhs;
        x *= rhs;
        return *this;
    }
    inline Vec2& operator*=(const Vec2& rhs) { return *this = (*this) * rhs; }
    inline Vec2& operator/=(const Vec2& rhs) { return *this = (*this) / rhs; }
    inline bool operator!=(const Vec2& rhs) const {
        return x != rhs.x || y != rhs.y;
    }
    inline bool operator==(const Vec2& rhs) const {
        return x == rhs.x && y == rhs.y;
    }
    inline void rotate(const double& rad) { *this = rotated(rad); }
    inline Vec2<double> rotated(const double& rad) const {
        return (*this) * rotation(rad);
    }
    static inline Vec2<double> rotation(const double& rad) {
        return Vec2(sin(rad), cos(rad));
    }
    static inline Vec2<double> rotation_deg(const double& deg) {
        return rotation(PI * deg / 180.0);
    }
    inline Vec2<double> rounded() const {
        return Vec2<double>(round(y), round(x));
    }
    inline Vec2<double> inv() const { // x + yj とみなす
        const double norm_sq = l2_norm_square();
        assert(norm_sq != 0.0);
        return Vec2(-y / norm_sq, x / norm_sq);
    }
    inline double l2_norm() const { return sqrt(x * x + y * y); }
    inline double l2_norm_square() const { return x * x + y * y; }
    inline T l1_norm() const { return std::abs(x) + std::abs(y); }
    inline double abs() const { return l2_norm(); }
    inline double phase() const { // [-PI, PI) のはず
        return atan2(y, x);
    }
    inline double phase_deg() const { // [-180, 180) のはず
        return phase() * (180.0 / PI);
    }
};
template <typename T, typename S>
inline Vec2<T> operator*(const S& lhs, const Vec2<T>& rhs) {
    return rhs * lhs;
}
template <typename T> ostream& operator<<(ostream& os, const Vec2<T>& vec) {
    os << vec.y << ' ' << vec.x;
    return os;
}

// 2 次元配列
template <class T, int height, int width> struct Board {
    array<T, height * width> data;
    template <class Int> constexpr inline auto& operator[](const Vec2<Int>& p) {
        return data[width * p.y + p.x];
    }
    template <class Int>
    constexpr inline const auto& operator[](const Vec2<Int>& p) const {
        return data[width * p.y + p.x];
    }
    template <class Int>
    constexpr inline auto& operator[](const initializer_list<Int>& p) {
        return data[width * *p.begin() + *(p.begin() + 1)];
    }
    template <class Int>
    constexpr inline const auto&
    operator[](const initializer_list<Int>& p) const {
        return data[width * *p.begin() + *(p.begin() + 1)];
    }
    constexpr inline void Fill(const T& fill_value) {
        fill(data.begin(), data.end(), fill_value);
    }
    void Print() const {
        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                cout << data[width * y + x] << " \n"[x == width - 1];
            }
        }
    }
};

inline double Time() {
    return static_cast<double>(
               chrono::duration_cast<chrono::nanoseconds>(
                   chrono::steady_clock::now().time_since_epoch())
                   .count()) *
           1e-9;
}

// まず、全点間距離を求める
// 次に、 2-opt する

static constexpr auto N = 100;
static constexpr auto M = 8;
static constexpr auto NM = N + M;

using Point = Vec2<float>;

// struct Point {
//     enum PointType {
//         kPlanet = 1,
//         kStation = 2,
//     };
//     Vec2<float> yx;
//     PointType type;

//     short id;

//     array<float>
// };

auto rng = mt19937();

struct State {
    array<Point, NM> points;
    double score;
    array<short, NM + 1> order;

    double entire_best_score;
    array<Vec2<short>, M> entire_best_points;
    vector<short> solution;

    void Read() {

        auto dummy = 0;
        cin >> dummy >> dummy;
        for (auto i = 0; i < N; i++) {
            cin >> points[i].y >> points[i].x;
        }
        for (auto i = 0; i < NM; i++) {
            order[i] = i;
        }
        order.back() = 0;
        entire_best_score = 1e30;
    }

    void Randomize() {
        for (auto i = 0; i < M; i++) {
            points[N + i].y = uniform_int_distribution<>(100, 900)(rng);
            points[N + i].x = uniform_int_distribution<>(100, 900)(rng);
        }
        shuffle(order.begin() + 1, order.end() - 1, rng);
    };

    void PrintBest() const {
        for (auto i = 0; i < M; i++) {
            cout << entire_best_points[i] << endl;
        }
        cout << solution.size() << endl;
        for (const auto p : solution) {
            if (p < 100) {
                cout << "1 " << p + 1 << endl;
            } else {
                cout << "2 " << p - 99 << endl;
            }
        }
    }
};

void TwoOpt(State& state) {

    auto& points = state.points;
    auto& order = state.order;

    auto distances = Board<float, NM, NM>();
    auto warshall_floyd_next = Board<short, NM, NM>();

    for (auto i = 0; i < NM; i++) {
        for (auto j = 0; j < NM; j++) {
            distances[{i, j}] = (points[i] - points[j]).l2_norm_square();
            if (i < N)
                distances[{i, j}] *= 5.0f;
            if (j < N)
                distances[{i, j}] *= 5.0f;
            warshall_floyd_next[{i, j}] = j;
        }
    }

    for (int k = 0; k < NM; k++)
        for (int i = 0; i < NM; i++)
            for (int j = 0; j < NM; j++)
                if (distances[{i, k}] + distances[{k, j}] < distances[{i, j}]) {
                    distances[{i, j}] = distances[{i, k}] + distances[{k, j}];
                    warshall_floyd_next[{i, j}] = warshall_floyd_next[{i, k}];
                }

    // スコア初期化
    auto& current_score = state.score;
    current_score = 0.0;
    for (auto i = 0; i < NM; i++) {
        current_score += distances[{order[i], order[i + 1]}];
    }

    for (auto iteration = 0; iteration < 100000; iteration++) {
        // 0 -> a -> b ...
        // 0 <- d <- c ...
        //        vvv
        // 0 -> a -> c ...
        // 0 <- d <- b ...
        auto idx_a = uniform_int_distribution<>(0, NM - 1)(rng);
        auto idx_c = uniform_int_distribution<>(0, NM - 1)(rng);
        if (abs(idx_a - idx_c) <= 1)
            continue;
        if (idx_a > idx_c) {
            swap(idx_a, idx_c);
        }
        const auto a = order[idx_a];
        const auto b = order[idx_a + 1];
        const auto c = order[idx_c];
        const auto d = order[idx_c + 1];

        const auto distance_before = distances[{a, b}] + distances[{c, d}];
        const auto distance_after = distances[{a, c}] + distances[{b, d}];
        const auto delta = distance_after - distance_before;
        if (delta < 0.0f) {
            reverse(order.begin() + idx_a + 1, order.begin() + idx_c + 1);
            current_score += delta;
        }
    }

    if (current_score >= state.entire_best_score)
        return;

    cerr << current_score << endl;

    state.entire_best_score = current_score;
    for (auto i = 0; i < M; i++) {
        state.entire_best_points[i] = {(short)points[N + i].y,
                                       (short)points[N + i].x};
    }

    // 経路復元
    // from https://zeosutt.hatenablog.com/entry/2015/05/05/045943
    state.solution.clear();
    for (auto idx_order = 0; idx_order < N + M; idx_order++) {
        const auto start = order[idx_order];
        const auto goal = order[idx_order + 1];
        for (short cur = start; cur != goal;
             cur = warshall_floyd_next[{cur, goal}])
            state.solution.push_back(cur);
    }
    state.solution.push_back(0);
}

void Solve() {
    auto state = State();
    state.Read();
    auto t0 = Time();
    auto outer_iteraiton = 0;
    while (Time() - t0 < 0.9) {
        state.Randomize();
        TwoOpt(state);
        outer_iteraiton++;
    }
    cerr << "outer_iteraiton=" << outer_iteraiton << endl;
    state.PrintBest();
}

int main() { Solve(); }