#include <cmath>
#include <functional>
#include <iostream>
#include <limits>
#include <queue>
#include <random>
#include <stack>
#include <vector>

#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"

using Rng = randutils::random_generator<pcg32>;

class Edge;

class Vertex {
public:
  unsigned index;
  double y;
  bool inZ;
  std::vector<std::reference_wrapper<Edge>> neighbors;
  std::array<unsigned, 2> coordinate;

  Vertex() : y(0), inZ(false){};

  void addEdge(Edge& e) { neighbors.push_back(e); }

  bool inMatching() const;
};

class Edge {
public:
  Vertex* s;
  Vertex* t;
  bool orientation;
  double weight;
  bool inPath;

  Edge() : orientation(false){};
  void setVertices(Vertex& red, Vertex& blue) {
    s = &red;
    t = &blue;
    red.addEdge(*this);
    blue.addEdge(*this);
  }
  Vertex& tail() {
    if (orientation) {
      return *t;
    } else {
      return *s;
    }
  }
  Vertex& head() {
    if (orientation) {
      return *s;
    } else {
      return *t;
    }
  }
  const Vertex& tail() const {
    if (orientation) {
      return *t;
    } else {
      return *s;
    }
  }
  const Vertex& head() const {
    if (orientation) {
      return *s;
    } else {
      return *t;
    }
  }
  bool isTight() const { return std::abs(s->y + t->y - weight) < 1e-13; }
};

bool Vertex::inMatching() const {
  for (const Edge& e : neighbors) {
    if (e.orientation) {
      return true;
    }
  }
  return false;
}

class Graph {
public:
  std::vector<Vertex> S;
  std::vector<Vertex> T;
  std::vector<Edge> E;

  Graph(unsigned n, Rng& r) : S(n * (n + 1)), T(n * (n + 1)), E(pow(2 * n, 2)) {
    for (unsigned i = 0; i < S.size(); i++) {
      S[i].index = i;
      S[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1};
      T[i].index = S.size() + i;
      T[i].coordinate = {2 * (i % n) + 1, 2 * (i / n)};
    }
    for (unsigned i = 0; i < E.size(); i++) {
      Vertex& redVertex = S[(1 + (i % (2 * n))) / 2 + (n + 1) * ((i / 4) / n)];
      Vertex& blueVertex = T[(i % (2 * n)) / 2 + n * (((i + 2 * n) / 4) / n)];
      E[i].setVertices(redVertex, blueVertex);
      E[i].weight = r.variate<double, std::exponential_distribution>(1);
    }
  }

  void hungarian() {
    std::queue<std::reference_wrapper<Vertex>> q;

    for (Vertex& v : T) {
      v.inZ = false;
    }

    for (Vertex& v : S) {
      v.inZ = false;

      if (!v.inMatching()) {
        q.push(v);
      }
    }

    while (!q.empty()) {
      Vertex& v = q.front();
      q.pop();

      v.inZ = true;

      for (Edge& e : v.neighbors) {
        if (e.isTight() && (e.tail().index == v.index && !e.head().inZ)) {
          q.push(e.head());
        }
      }
    }

    bool reversedPath = false;

    for (Vertex& v : T) {
      if (v.inZ && !v.inMatching()) {
        std::stack<std::tuple<std::vector<std::reference_wrapper<Edge>>::iterator,
                              std::vector<std::reference_wrapper<Edge>>::iterator, unsigned>>
            path;

        for (Edge& e : E) {
          e.inPath = false;
        }

        path.push({v.neighbors.begin(), v.neighbors.end(), v.index});

        while (!path.empty()) {
          auto [it, itEnd, i] = path.top();
          if (it == itEnd) {
            path.pop();
            if (path.empty()) {
              exit(1);
            }
            std::tie(it, itEnd, i) = path.top();
            it->get().inPath = false;
            it++;
            path.top() = {it, itEnd, i};
          } else {
            Edge& e = it->get();

            if (e.isTight() && (e.head().index == i && !e.inPath)) {
              if (!e.tail().inMatching()) {
                break;
              } else {
                e.inPath = true;
                path.push({e.tail().neighbors.begin(), e.tail().neighbors.end(), e.tail().index});
              }
            } else {
              it++;
              path.top() = {it, itEnd, i};
            }
          }
        }

        while (!path.empty()) {
          auto [it, itEnd, i] = path.top();
          Edge& e = it->get();
          path.pop();
          e.orientation = !e.orientation;
        }

        reversedPath = true;
        break;
      }
    }

    if (!reversedPath) {
      double Δ = std::numeric_limits<double>::infinity();

      for (Edge& e : E) {
        if (e.s->inZ && !e.t->inZ) {
          Δ = std::min(Δ, e.weight - e.s->y - e.t->y);
        }
      }

      for (Vertex& v : S) {
        if (v.inZ) {
          v.y += Δ;
        }
      }

      for (Vertex& v : T) {
        if (v.inZ) {
          v.y -= Δ;
        }
      }
    }
  }

  bool isPerfect() const {
    for (const Vertex& v : S) {
      if (!v.inMatching()) {
        return false;
      }
    }
    for (const Vertex& v : T) {
      if (!v.inMatching()) {
        return false;
      }
    }

    return true;
  }
};

int main(int argc, char* argv[]) {
  unsigned n = 100;
  unsigned maxSteps = 1e8;
  double beliefThreshold = 1;

  int opt;

  while ((opt = getopt(argc, argv, "n:m:t:")) != -1) {
    switch (opt) {
    case 'n':
      n = atoi(optarg);
      break;
    case 'm':
      maxSteps = (unsigned)atof(optarg);
      break;
    case 't':
      beliefThreshold = atof(optarg);
      break;
    default:
      exit(1);
    }
  }

  Rng r;
  Graph G(n, r);

  unsigned steps = 0;
  double minBelief = 0;

  while (!G.isPerfect()) {
    G.hungarian();
  }

  for (const Edge& e : G.E) {
    if (e.orientation) {
      std::cout << e.tail().coordinate[0] << " " << e.tail().coordinate[1] << " "
                << e.head().coordinate[0] << " " << e.head().coordinate[1] << std::endl;
    }
  }

  return 0;
}