namespace atcoder {}

#ifdef LOCAL
#define dbg(x) cerr << __LINE__ << " : " << #x << " = " << (x) << endl;
#else
#define NDEBUG
#define dbg(x) true;
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#endif

#ifdef GTEST
#include <gtest/gtest.h>
#endif

#include <math.h>

#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#ifdef PERF
#include <gperftools/profiler.h>
#endif

using namespace std;
using namespace atcoder;
#define fast_io                     \
  ios_base::sync_with_stdio(false); \
  cin.tie(0);                       \
  cout.tie(0);
#define ll long long int
#define rep(i, n) for (int i = 0; i < (int)(n); i++)
#define reps(i, n) for (int i = 1; i <= (int)(n); i++)
#define REP(i, n) for (int i = n - 1; i >= 0; i--)
#define REPS(i, n) for (int i = n; i > 0; i--)
#define MOD (long long int)(1e9 + 7)
#define INF (int)(1e9)
#define LINF (long long int)(1e18)
#define chmax(a, b) a = (((a) < (b)) ? (b) : (a))
#define chmin(a, b) a = (((a) > (b)) ? (b) : (a))
#define all(v) v.begin(), v.end()
typedef pair<int, int> Pii;
typedef pair<ll, ll> Pll;
constexpr double PI = acos(-1);

#ifdef NDEBUG
#define CHECK(v1, op, v2)
#else
#define CHECK(v1, op, v2)                            \
  if (!((v1)op(v2))) {                               \
    cerr << "ERROR:" << (v1) << " " << (v2) << endl; \
    assert((v1)op(v2));                              \
  }
#endif

long double nCr(const int n, const int r) {
  long double ret = 1;
  rep(t, r) {
    ret *= (n - t);
    ret /= (r - t);
  }
  return ret;
}

template <typename T>
string to_string(const vector<T>& vec) {
  string ret = "";
  rep(i, vec.size()) {
    ret += vec[i].to_string();
    if (i + 1 != vec.size()) {
      ret += ",";
    }
  }
  return ret;
}

template <typename T>
ostream& operator<<(ostream& os, const vector<T>& vec) {
  os << to_string(vec);
  return os;
}

uint32_t xorshift() {
  static uint32_t x = 12345789;
  static uint32_t y = 362436069;
  static uint32_t z = 521288629;
  static uint32_t w = 88675123;
  uint32_t t;
  t = x ^ (x << 11);
  x = y;
  y = z;
  z = w;
  w ^= t ^ (t >> 8) ^ (w >> 19);

  return w;
}

int rand(const uint32_t l, const uint32_t r) {
  return xorshift() % (r - l) + l;
}

uint32_t rand_other_than(const uint32_t l, const uint32_t r,
                         const uint32_t other) {
  const uint32_t num = rand(l, r - 1);
  return num + (num >= other);
}

template <typename T>
const T& rand_vec(const vector<T>& vec) {
  assert(vec.size() > 0);
  return vec[rand(0, vec.size())];
}

template <typename T>
void shuffle(vector<T>& vec) {
  rep(l, (int)vec.size() - 1) {
    const int idx = rand(l, vec.size());
    swap(vec[idx], vec[l]);
  }
}




/* start */

vector<double> PARAMS = {};





class Timer {
  chrono::system_clock::time_point _start, _end;
  ll _sum = 0, _count = 0;

 public:
  void start() { _start = chrono::system_clock::now(); }

  void stop() { _end = chrono::system_clock::now(); }

  void add() {
    const chrono::system_clock::time_point now = chrono::system_clock::now();
    _sum += static_cast<double>(
        chrono::duration_cast<chrono::nanoseconds>(now - _start).count());
    _count++;
  }

  ll sum() const { return _sum / 1000; }

  int count() const { return _count; }

  string average() const {
    if (_count == 0) {
      return "NaN";
    }
    return to_string(_sum / 1000 / _count);
  }

  void reset() {
    _start = chrono::system_clock::now();
    _sum = 0;
    _count = 0;
  }

  inline int ms() const {
    const chrono::system_clock::time_point now = chrono::system_clock::now();
    return static_cast<double>(
        chrono::duration_cast<chrono::microseconds>(now - _start).count() /
        1000);
  }

  inline int ns() const {
    const chrono::system_clock::time_point now = chrono::system_clock::now();
    return static_cast<double>(
        chrono::duration_cast<chrono::microseconds>(now - _start).count());
  }
};

#ifdef LOCAL
struct Timers : unordered_map<string, Timer> {
  friend ostream& operator<<(ostream& os, const Timers& timers) {
    for (const auto& pa : timers) {
      os << pa.first << " time: " << pa.second.sum() / 1000
         << " count: " << pa.second.count() << endl;
    }
    return os;
  }
};
#else
struct Timers {
  struct Dummy {
    void start() const {}
    void add() const {}
  };
  Dummy dummy;
  const Dummy& operator[](const std::string& str) { return dummy; }
  friend ostream& operator<<(ostream& os, const Timers& timers) { return os; }
};
#endif

Timers global_timers;


// NBD = neighborhood
template <class Solution, class NBDGenerator>
struct SimulatedAnnealing {
 public:
  SimulatedAnnealing(double end_time, double start_temp, double end_temp)
      : end_time_(end_time * 1000),  // ns
        inv_time_(1.0 / end_time_),
        cur_time_(0.0),
        start_temp_(start_temp),
        end_temp_(end_temp),
        diff_temp_(end_temp - start_temp),
        cur_temp_(start_temp) {
    assert(start_temp >= end_temp);
    timer_.start();
  }

  Solution run(const Solution& initial_sol, NBDGenerator& nbd_generator) {
    Solution sol(initial_sol);

    int acc = 0;
    int g = 0;
    for (;; ++g) {
      if ((g & 0xf) == 0) {
#ifdef SA_MAX_G
        if (g >= SA_MAX_G) {
          break;
        }
#endif
        UpdateTime();
        UpdateTemp();
        if (cur_time_ > end_time_) {
          break;
        }
      }

      const auto nbd = nbd_generator.Generate(&sol, g);
      constexpr int mod = 1'000'000'000;
      const int r = rand(0, mod);
      if (static_cast<double>(r) / mod < exp(-nbd->Diff() / cur_temp_)) {
        acc++;
        nbd->Update(&sol);
      } else {
        nbd->Restore(&sol);
      }
    }
    cerr << "g,acc: " << g << " " << acc << endl;
    return sol;
  }

 private:
  double end_time_, inv_time_, cur_time_;
  double start_temp_, end_temp_, diff_temp_, cur_temp_;
  Timer timer_;

  void UpdateTime() { cur_time_ = timer_.ns(); }
  void UpdateTemp() {
    cur_temp_ = max(0.0, start_temp_ - cur_time_ * diff_temp_ * inv_time_);
  }
};




/* start */

template <typename Point, typename Dist_T>
class Graph_T {
 public:
  struct Edge {
    Point to;
    Dist_T dist;

    Edge(const Point _to, const Dist_T _dist) : to(_to), dist(_dist) {}
  };
  Graph_T() {}
  explicit Graph_T(const int n)
      : n_(n), adjs_(n), dists_(n, vector<Dist_T>(n)) {}
  int N() const { return n_; }
  Dist_T Dist(const Point& from, const Point& to) const {
    return dists_[from][to];
  }
  void AddEdge(const Point a, const Point b, const Dist_T dist) {
    adjs_[a].emplace_back(b, dist);
    dists_[a][b] = dist;
  }
  const vector<Edge>& Adj(const Point a) const { return adjs_[a]; }
  void RemoveEdge(const Point a, const Point b) {
    {
      const auto itr = std::find_if(
          all(adjs_[a]), [&](const Edge& edge) { return edge.to == b; });
      assert(itr != adjs_[a].end());
      adjs_[a].erase(itr);
      dists_[a][b] = INF;
    }
  }

 private:
  int n_;
  vector<vector<Edge>> adjs_;
  vector<vector<Dist_T>> dists_;
};



enum Dir {
  // y-1, y+1, x-1, x+1
  kU,
  kD,
  kL,
  kR
};

struct Pos {
  int idx_;
  Pos() {}
  explicit Pos(const int _idx) : idx_(_idx) {}
  Pos(int _x, int _y) : idx_(_y * N + _x) { assert(N > 0); }
  int X() const { return pos_2_x[*this]; }
  int Y() const { return pos_2_y[*this]; }
  int Idx() const { return idx_; }
  operator int() const { return Idx(); }
  operator size_t() const { return Idx(); }
  const vector<Pos>& Adj() const { return adj_poses[*this]; }
  const vector<Dir>& AdjDirs() const { return adj_dirs[*this]; }
  int Manhattan(const Pos& other) const {
    return abs(X() - other.X()) + abs(Y() - other.Y());
  }
  int Euclid2(const Pos& other) const {
    const int dx = X() - other.X();
    const int dy = Y() - other.Y();
    return dx * dx + dy * dy;
  }
  bool Move(const Dir dir) {
    if (move_to[dir][*this].IsDummy()) {
      return false;
    } else {
      *this = move_to[dir][*this];
      return true;
    }
  }
  bool operator<(const Pos& other) const { return this->Idx() < other.Idx(); }
  bool operator==(const Pos& other) const { return this->Idx() == other.Idx(); }
  bool operator!=(const Pos& other) const { return this->Idx() != other.Idx(); }
  friend ostream& operator<<(ostream& os, const Pos& pos) {
    os << pos.X() << " " << pos.Y();
    return os;
  }
  bool IsDummy() const { return this->idx_ < 0; }
  static Pos Dummy() {
    Pos p;
    p.idx_ = -1;
    return p;
  }

  static Pos TryCreate(const int x, const int y, const int z) {
    if (y < 0 || y >= N || x < 0 || x >= N) {
      return Pos::Dummy();
    } else {
      return Pos(x, y);
    }
  }

  static void StaticInit(const int n) {
    N = n;
    N2 = N * N;
    N4 = N2 * N2;
    adj_poses = vector<vector<Pos>>(N2);
    adj_dirs = vector<vector<Dir>>(N2);
    move_to = vector<vector<Pos>>(4, vector<Pos>(N2, Pos::Dummy()));
    pos_2_y.resize(N2);
    pos_2_x.resize(N2);

    rep(y, N) {
      rep(x, N) {
        const Pos p{x, y};
        pos_2_y[p] = y;
        pos_2_x[p] = x;
        for (int dy = -1; dy <= 1; ++dy) {
          for (int dx = -1; dx <= 1; ++dx) {
            if (abs(dy) + abs(dx) != 1) continue;
            const int adj_y = y + dy;
            const int adj_x = x + dx;
            if (adj_y < 0 || adj_y >= N) continue;
            if (adj_x < 0 || adj_x >= N) continue;
            adj_poses[p].emplace_back(adj_x, adj_y);

            Dir dir;
            if (dy == -1) {
              dir = kU;
            } else if (dy == 1) {
              dir = kD;
            } else if (dx == -1) {
              dir = kL;
            } else if (dx == 1) {
              dir = kR;
            } else {
              assert(false);
            }

            move_to[dir][p] = {adj_x, adj_y};
            adj_dirs[p].emplace_back(dir);
          }
        }
      }
    }
  }
  static int N, N2, N4;
  static vector<vector<Pos>> adj_poses;
  static vector<vector<Dir>> adj_dirs;
  static vector<vector<Pos>> move_to;
  static vector<int> pos_2_x, pos_2_y;
};

vector<vector<Pos>> Pos::adj_poses;
vector<vector<Dir>> Pos::adj_dirs;
vector<vector<Pos>> Pos::move_to;
vector<int> Pos::pos_2_x, Pos::pos_2_y;
int Pos::N, Pos::N2, Pos::N4;








/* start */

constexpr int kPlanetType = 1;
constexpr int kStationType = 2;

struct Solution {
  using Graph = Graph_T<int, int>;
  Solution() {
    orders.resize(N);
    std::iota(all(orders), 0);
    stations.resize(M);
    rep(i, M) { stations[i] = Pos{rand(0, Pos::N + 1), rand(0, Pos::N + 1)}; }
  }

  static void StaticInit(const int _N, const int _M,
                         const vector<Pos>& _planets, const Graph& _graph) {
    N = _N;
    M = _M;
    planets = _planets;
    graph = _graph;
  }

  // [l, r)のreverse
  int DiffSwap(const int l, const int r) const {
    // l-1, l, r-1, r
    // l-1, r-1, l, r
    const int i0 = (l - 1 + N) % N;
    const int i1 = l % N;
    const int i2 = (r - 1 + N) % N;
    const int i3 = r % N;
    const int prev_dist =
        graph.Dist(orders[i0], orders[i1]) + graph.Dist(orders[i2], orders[i3]);
    const int next_dist =
        graph.Dist(orders[i0], orders[i2]) + graph.Dist(orders[i1], orders[i3]);
    return next_dist - prev_dist;
  }

  void Swap(int l, int r) {
    if (l > r) {
      swap(l, r);
    }

    reverse(orders.begin() + l, orders.begin() + r);
  }

  vector<pair<int, int>> CalcAns() const {
    int start_i = 0;
    for (start_i = 0; start_i < N; start_i++) {
      if (this->orders[start_i] == 0) {
        break;
      }
    }

    vector<pair<int, int>> ans;
    ans.emplace_back(kPlanetType, this->orders[start_i]);
    rep(add_i, N) {
      const int i = (start_i + add_i) % N;
      const int next_i = (i + 1) % N;
      const Pos& from_planet = planets[this->orders[i]];
      const Pos& to_planet = planets[this->orders[next_i]];

      // stationを検討
      int min_dist =
          Solution::graph.Dist(this->orders[i], this->orders[next_i]);
      vector<int> best_stations;

      // i,jの順に経由
      rep(i, M) {
        const Pos& station_i = this->stations[i];
        rep(j, M) {
          if (i == j) continue;
          const Pos& station_j = this->stations[j];
          const int dist = from_planet.Euclid2(station_i) * alpha +
                           station_i.Euclid2(station_j) +
                           station_j.Euclid2(to_planet) * alpha;
          if (dist < min_dist) {
            min_dist = dist;
            best_stations = {i, j};
          }
        }
      }

      // iのみ経由
      rep(i, M) {
        const Pos& station = this->stations[i];
        const int dist = from_planet.Euclid2(station) * alpha +
                         station.Euclid2(to_planet) * alpha;
        if (dist < min_dist) {
          min_dist = dist;
          best_stations = {i};
        }
      }

      if (best_stations.size() == 0) {
        // 惑星to惑星
        int from = this->orders[i];
        int to = this->orders[next_i];
        vector<int> rev_planets;
        while (from != to) {
          rev_planets.emplace_back(to);
          rep(i, N) {
            if (i == to) continue;
            if (graph.Dist(from, i) + graph.Dist(i, to) ==
                graph.Dist(from, to)) {
              to = i;
              break;
            }
          }
        }
        REP(i, rev_planets.size()) {
          ans.emplace_back(kPlanetType, rev_planets[i]);
        }
      } else {
        for (const auto station : best_stations) {
          ans.emplace_back(kStationType, station);
        }
        ans.emplace_back(kPlanetType, this->orders[next_i]);
      }
    }

    return ans;
  }

  friend ostream& operator<<(ostream& os, const Solution& sol) {
    const auto ans = sol.CalcAns();

    rep(i, M) { os << sol.stations[i] << endl; }

    os << ans.size() << endl;
    rep(i, ans.size()) {
      const auto& [t, r] = ans[i];
      os << t << " " << r + 1;
      if (i + 1 != (int)ans.size()) {
        os << endl;
      }
    }

    int dist_sum = 0;
    rep(i, (int)ans.size() - 1) {
      const auto& [t0, r0] = ans[i];
      const auto& [t1, r1] = ans[(i + 1) % ans.size()];
      const Pos& from = (t0 == kPlanetType ? planets[r0] : sol.stations[r0]);
      const Pos& to = (t1 == kPlanetType ? planets[r1] : sol.stations[r1]);
      dist_sum += from.Euclid2(to) * (t0 == kPlanetType ? alpha : 1) *
                  (t1 == kPlanetType ? alpha : 1);
    }
    cerr << "score: " << std::round(1e9 / (1000 + sqrt(dist_sum))) << endl;

    return os;
  }

  static vector<Pos> planets;
  static int N, M;
  static Graph graph;
  static constexpr int alpha = 5;
  static constexpr int alpha2 = 25;

  vector<int> orders;
  vector<Pos> stations;
};

int Solution::N;
int Solution::M;
vector<Pos> Solution::planets;
Graph_T<int, int> Solution::graph;







/* start */

struct TSPSolution {
  using Graph = Graph_T<int, int>;
  TSPSolution() {
    orders.resize(N);
    std::iota(all(orders), 0);
  }

  TSPSolution(const vector<int>& _orders) : orders(_orders) {}

  static void StaticUpdateGraph(const Graph& _graph) {
    graph = _graph;
    N = _graph.N();
  }

  // [l, r)のreverse
  int DiffSwap(const int l, const int r) const {
    // l-1, l, r-1, r
    // l-1, r-1, l, r
    const int i0 = (l - 1 + N) % N;
    const int i1 = l % N;
    const int i2 = (r - 1 + N) % N;
    const int i3 = r % N;
    const int prev_dist =
        graph.Dist(orders[i0], orders[i1]) + graph.Dist(orders[i2], orders[i3]);
    const int next_dist =
        graph.Dist(orders[i0], orders[i2]) + graph.Dist(orders[i1], orders[i3]);
    return next_dist - prev_dist;
  }

  void Swap(int l, int r) {
    if (l > r) {
      swap(l, r);
    }

    reverse(orders.begin() + l, orders.begin() + r);
  }

  // toの前にtargetを入れる
  int DiffInsert(const int to, const int target) const {
    // to-1, to, target-1, target, target, target+1
    // to-1, target, target, to, target-1, target+1
    const int i0 = (to - 1 + N) % N;
    const int i1 = to;
    const int i2 = (target - 1 + N) % N;
    const int i3 = target;
    const int i4 = target;
    const int i5 = (target + 1) % N;

    const int prev_dist = graph.Dist(orders[i0], orders[i1]) +
                          graph.Dist(orders[i2], orders[i3]) +
                          graph.Dist(orders[i4], orders[i5]);
    const int next_dist = graph.Dist(orders[i0], orders[i3]) +
                          graph.Dist(orders[i4], orders[i1]) +
                          graph.Dist(orders[i2], orders[i5]);
    return next_dist - prev_dist;
  }

  void Insert(const int to, const int target) {
    orders.insert(orders.begin() + to, orders[target]);
    orders.erase(orders.begin() + target + (to <= target ? 1 : 0));
  }

  static Graph graph;
  static int N;

  vector<int> orders;
};

int TSPSolution::N;
Graph_T<int, int> TSPSolution::graph;


/* start */

using Graph = Graph_T<int, int>;

struct NBD {
  int Diff() const { return diff_; }
  virtual void Update(TSPSolution*) const = 0;
  void Restore(TSPSolution*) const {}

  int diff_;
};

struct NBD_swap : public NBD {
  NBD_swap() {}
  NBD_swap(TSPSolution* const sol, const int l, const int r) {
    l_ = l;
    r_ = r;
    diff_ = sol->DiffSwap(l, r);
  }
  void Update(TSPSolution* sol) const override { sol->Swap(l_, r_); }

  int l_;
  int r_;
};

struct NBD_insert : public NBD {
  NBD_insert() {}
  NBD_insert(TSPSolution* const sol, const int to, const int target) {
    to_ = to;
    target_ = target;
    diff_ = sol->DiffInsert(to, target);
  }
  void Update(TSPSolution* sol) const override { sol->Insert(to_, target_); }

  int to_;
  int target_;
};

struct NBDGenerator {
  NBDGenerator() {}
  NBD* Generate(TSPSolution* const sol, const int g) {
    if ((g & 0x4) == 0) {
      // 1.5-opt
      const int to = rand(0, TSPSolution::N);
      const int target = rand(0, TSPSolution::N);
      nbd_insert0 = NBD_insert(sol, to, target);
      return &nbd_insert0;
    } else {
      // 2-opt
      const int l = rand(0, TSPSolution::N);
      const int r = rand_other_than(0, TSPSolution::N, l);
      nbd_swap0 = NBD_swap(sol, l, r);
      return &nbd_swap0;
    }
  }

  NBD_swap nbd_swap0;      // 2-opt
  NBD_insert nbd_insert0;  // 1.5-opt
};

struct State {
  static State Dummy() {
    State state;
    state.sum_dist = INF;
    return state;
  }

  static void StaticInit(const int _N, const int _M,
                         const vector<Pos>& _planets) {
    N = _N;
    M = _M;
    planets = _planets;
  }

  bool IsDummy() const { return sum_dist == INF; }

  void UpdateWarshall() {
    if (stations.size() == 0) return;
    const int total_N = N + stations.size();
    const int k = total_N - 1;
    const Pos last_station = stations.back();

    // kにまっすぐ向かった時の距離
    vector<int> dist_to_k_straight(total_N);
    rep(i, N) {
      dist_to_k_straight[i] = planets[i].Euclid2(last_station) * alpha;
    }
    for (int i = N; i < total_N; ++i) {
      dist_to_k_straight[i] = stations[i - N].Euclid2(last_station);
    }

    // iからkへの最短距離をまずは求める
    rep(i, total_N) {
      int min_dist = INF;
      rep(j, total_N - 1) {
        const int dist = warshall[i][j] + dist_to_k_straight[j];
        chmin(min_dist, dist);
      }
      assert(min_dist < INF);
      warshall[i][k] = warshall[k][i] = min_dist;
    }

    // warshallを更新
    warshall[k][k] = 0;
    rep(i, total_N) {
      rep(j, i) {
        const int next_dist = warshall[i][k] + warshall[k][j];
        if (next_dist < warshall[i][j]) {
          warshall[i][j] = warshall[j][i] = next_dist;
        }
      }
    }
  }

  int CalcSumDist() const {
    int ret = 0;
    rep(i, orders.size()) {
      const int next_i = (i + 1) % orders.size();
      ret += warshall[orders[i]][orders[next_i]];
    }
    return ret;
  }

  void TwoOpt(const double time_limit) {
    SimulatedAnnealing<TSPSolution, NBDGenerator> sa(time_limit, 1.0, 0.03);
    NBDGenerator nbd_generator;
    Graph graph(N);
    rep(i, N) {
      rep(j, N) { graph.AddEdge(i, j, warshall[i][j]); }
    }
    TSPSolution::StaticUpdateGraph(graph);
    assert(orders.size() == N);
    TSPSolution initial_sol(orders);
    const auto tsp_sol = sa.run(initial_sol, nbd_generator);
    orders = tsp_sol.orders;

    sum_dist = CalcSumDist();
  }

  vector<vector<int>> warshall;
  vector<int> orders;
  vector<Pos> stations;
  int sum_dist;
  int prev_j;
  int prev_k;

  static constexpr int alpha = 5;
  static constexpr int alpha2 = 25;

  static vector<Pos> planets;
  static int N, M;
};

int State::N;
int State::M;
vector<Pos> State::planets;

class Solver {
 public:
  Solver(istream& is) {
    is >> N >> M;
    Pos::StaticInit(1001);
    planets.reserve(N);
    rep(i, N) {
      int a, b;
      is >> a >> b;
      planets.emplace_back(a, b);
    }

    warshall_without_stations = CalcWarshallWithoutStations();
    const auto graph = InitGraph(warshall_without_stations);
    TSPSolution::StaticUpdateGraph(graph);
    Solution::StaticInit(N, M, planets, graph);
    State::StaticInit(N, M, planets);
  }

  vector<vector<int>> CalcWarshallWithoutStations() const {
    // stationの分も確保しておく
    vector<vector<int>> dists(N + M, vector<int>(N + M, INF));
    rep(i, N + M) { dists[i][i] = 0; }
    rep(i, N) {
      rep(j, N) { dists[i][j] = alpha2 * planets[i].Euclid2(planets[j]); }
    }

    rep(k, N) {
      rep(i, N) {
        rep(j, N) { chmin(dists[i][j], dists[i][k] + dists[k][j]); }
      }
    }
    return dists;
  }

  Graph InitGraph(const vector<vector<int>>& dists) const {
    // dists.size() = N+M
    // graph.size() = N
    Graph graph(N);
    rep(i, N) {
      rep(j, N) {
        const int dist = dists[i][j];
        graph.AddEdge(i, j, dist);
      }
    }

    return graph;
  }

  using BS = bitset<100>;

  vector<pair<Pos, int>> GridSearchNextStation(const State& state,
                                               const int delta) const {
    vector<pair<Pos, int>> best_next_stations(N + 1, {Pos::Dummy(), INF});
    assert(state.orders.size() == N);
    for (int y = 0; y <= 1000; y += delta) {
      for (int x = 0; x <= 1000; x += delta) {
        const Pos station(x, y);
        // このstationの経由を可能としたときの、各辺のdistを求める
        int dist_sum = 0;
        int used_station_count = 0;
        rep(order_i, N) {
          const int i = state.orders[order_i];
          const int j = state.orders[(order_i + 1) % N];
          const auto pos_i = planets[i];
          const auto pos_j = planets[j];
          const int used_dist =
              (pos_i.Euclid2(station) + pos_j.Euclid2(station)) * alpha;
          const int dist = min(state.warshall[i][j], used_dist);
          dist_sum += dist;

          used_station_count += (dist < warshall_without_stations[i][j]);
        }

        if (best_next_stations[used_station_count].second > dist_sum) {
          best_next_stations[used_station_count] = {station, dist_sum};
        }
      }
    }

    return best_next_stations;
  }

  Solution BeamSearch(const vector<int>& initial_orders) const {
    assert(initial_orders.size() == N);
    assert(warshall_without_stations.size() == N + M);
    vector dp(M + 1, vector(N + 1, State::Dummy()));
    dp[0][0].warshall = warshall_without_stations;
    dp[0][0].orders = initial_orders;
    dp[0][0].sum_dist = dp[0][0].CalcSumDist();

    rep(j, M + 1) {
      rep(k, N + 1) {
        auto& state = dp[j][k];
        if (state.IsDummy()) continue;

        // TSPを追加で1ms行い、既存の駅を有効活用できている状態にする
        if (j != 0) {
          state.warshall = dp[state.prev_j][state.prev_k].warshall;
        }
        state.UpdateWarshall();
        state.TwoOpt(0.8);

        if (j == M) continue;

        // 訪問順は固定した状態で、次の駅を置く場所をグリッドサーチして、
        //「何らかの駅を経由している辺の本数」ごとに「total_distが最小」となるような
        // 駅の置き場所を求める
        const auto best_next_stations = GridSearchNextStation(state, 33);

        // それらを用いて遷移
        rep(next_k, N + 1) {
          const auto& [station, dist] = best_next_stations[next_k];
          if (dist == INF) continue;

          auto& next_state = dp[j + 1][next_k];
          if (next_state.sum_dist <= dist) continue;
          next_state.orders = state.orders;
          next_state.stations = state.stations;
          next_state.stations.emplace_back(station);
          next_state.sum_dist = dist;
          next_state.prev_j = j;
          next_state.prev_k = k;
        }
      }
    }

    using IDX = tuple<int, int>;
    IDX best_idx;
    int min_dist = INF;
    const int j = M;
    rep(k, N + 1) {
      auto& state = dp[j][k];
      if (state.IsDummy()) continue;

      if (state.sum_dist < min_dist) {
        min_dist = state.sum_dist;
        best_idx = {j, k};
      }
    }

#if 0
    IDX idx = best_idx;
    while (true) {
      const auto [j, k] = idx;
      const auto& state = dp[j][k];

      Solution sol;
      sol.orders = state.orders;
      sol.stations = state.stations;
      sol.stations.resize(M, Pos{0, 0});
      cout << endl << endl << sol << endl;

      if (j == 0) {
        break;
      }

      idx = {state.prev_j, state.prev_k};
    }
#endif

    cerr << "min_dist: " << min_dist << endl;

    const auto& [best_j, best_k] = best_idx;
    Solution best_sol;
    best_sol.orders = dp[best_j][best_k].orders;
    best_sol.stations = dp[best_j][best_k].stations;
    return best_sol;
  }

  Solution OptimizeStations(const Solution& base_sol) const {
    auto sol = base_sol;
    rep(_t, 3) {
      const auto ans = base_sol.CalcAns();
      // 各stationから伸びている辺のdistの総和を最小化する位置に各stationを移動する

      // 重心に移動すればよい
      // station->stationについては、重みを1/5にするのが自然だが、
      // もう一つのstationもずれることが予想されるので、重みは1/10にする
      vector<float> weight_sums(M);
      vector<float> Xs(M), Ys(M);

      rep(i, (int)ans.size() - 1) {
        auto [t0, r0] = ans[i];
        auto [t1, r1] = ans[i + 1];
        Pos p0 = (t0 == kPlanetType ? planets[r0] : sol.stations[r0]);
        Pos p1 = (t1 == kPlanetType ? planets[r1] : sol.stations[r1]);
        rep(_t, 2) {
          swap(t0, t1);
          swap(r0, r1);
          swap(p0, p1);

          // 0について考える
          if (t0 != kStationType) continue;
          // station->planetなら1.0
          // station->stationなら0.1
          const float weight = (t1 == kPlanetType ? 1.0 : 0.1);
          weight_sums[r0] += weight;
          Xs[r0] += p1.X() * weight;
          Ys[r0] += p1.Y() * weight;
        }
      }

      rep(i, M) {
        const auto weight = weight_sums[i];
        assert(weight > 0.0);
        float float_x = Xs[i] / weight;
        float float_y = Ys[i] / weight;
        const int x = std::clamp((int)std::round(float_x), 0, 1000);
        const int y = std::clamp((int)std::round(float_y), 0, 1000);
        sol.stations[i] = Pos{x, y};
      }

      // warshallを0から更新
      vector<vector<int>> dists(N + M, vector<int>(N + M, INF));
      vector<Pos> planets_and_stations(N + M);
      rep(i, N) { planets_and_stations[i] = planets[i]; }
      for (int i = N; i < N + M; i++) {
        planets_and_stations[i] = sol.stations[i - N];
      }
      rep(i, N + M) { dists[i][i] = 0; }
      rep(i, N + M) {
        rep(j, i) {
          dists[i][j] = dists[j][i] =
              planets_and_stations[i].Euclid2(planets_and_stations[j]) *
              (i < N ? alpha : 1) * (j < N ? alpha : 1);
        }
      }

      rep(k, N + M) {
        rep(i, N + M) {
          rep(j, N + M) { chmin(dists[i][j], dists[i][k] + dists[k][j]); }
        }
      }

      const auto graph = InitGraph(dists);
      TSPSolution::StaticUpdateGraph(graph);
      SimulatedAnnealing<TSPSolution, NBDGenerator> sa(10, 1.0, 0.03);
      NBDGenerator nbd_generator;
      TSPSolution initial_sol(sol.orders);
      const auto tsp_sol = sa.run(initial_sol, nbd_generator);
      sol.orders = tsp_sol.orders;
    }

    return sol;
  }

  Solution Solve(const int time_limit) const {
    const TSPSolution initial_sol;
    SimulatedAnnealing<TSPSolution, NBDGenerator> sa(100, 1.0, 0.03);
    NBDGenerator nbd_generator;
    const auto tsp_sol = sa.run(initial_sol, nbd_generator);
    const auto beam_sol = BeamSearch(tsp_sol.orders);
    const auto optimized_sol = OptimizeStations(beam_sol);

    return optimized_sol;
  }

 private:
  int N, M;
  vector<vector<int>> warshall_without_stations;
  vector<Pos> planets;
  static constexpr int alpha = 5;
  static constexpr int alpha2 = 25;
};

int main(int argc, char* argv[]) {
  fast_io;

  if (argc >= 2) {
    int idx = 0;
    for (int i = 1; i < argc; ++i) {
      PARAMS[idx++] = std::stod(argv[i]);
    }
  }

  Timer timer;
  timer.start();
  Solver solver(cin);
  auto sol = solver.Solve(900 - timer.ms());
  cerr << "ms: " << timer.ms() << endl;
  cout << sol << endl;
  return 0;
}