<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

190. Dominoes
time limit per test: 0.25 sec.
memory limit per test: 4096 KB
input: standard input
output: standard output



There is a NxN squared chessboard (1<=N<=40). P squares were removed from the chessboard (0<=P<N*N). It is necessary to find out whether it is possible to cover the remaining part of the chessboard with dominoes dice (each die is 2x1 squares sized). Each die should lie on two neighboring squares exactly. No two dice can cover one and the same square. Your task is to find the covering, if it exists.

Input
The first line contains two integer numbers N and P separated by a space. Each of the following P lines contains a pair of numbers separated by a space - coordinates of the removed square (1<=Xi, Yi<=N). The bottom left square has the coordinates (1, 1), the bottom right square - (N, 1).

Output
If the covering exists, output "Yes" to the first line and "No" in the opposite case. If the first answer was positive, then output to the second line integer number Nh - the number of horizontally placed dice. Each of the following Nh lines should contain two integers - the coordinates of the left square covered by a corresponding die. Output to the next line Nv - the number of vertically placed dice. And the following Nv lines should contain the coordinates of the bottom square covered by a corresponding die.

Sample test(s)

Input
4 10
1 3
1 2
1 1
2 1
3 1
4 1
3 2
4 2
3 3
4 3

Output
Yes
2
1 4
3 4
1
2 2
Author:	Michael R. Mirzayanov
Resource:	ACM International Collegiate Programming Contest 2003-2004
North-Eastern European Region, Southern Subregion
Date:	2003 October, 9

<|response|>
1. Abridged problem statement
We have an N×N grid (1 ≤ N ≤ 40). P cells are removed (0 ≤ P < N²). Determine whether the remaining cells can be exactly covered by non-overlapping 2×1 dominoes placed on orthogonally adjacent cells. If a tiling exists, print "Yes", then output the lists of dominoes (vertical by the coordinates of each domino's bottom cell, horizontal by the coordinates of each domino's left cell, each list preceded by its count). Otherwise print "No".

2. Key observations needed to solve the problem
  • On a checkerboard coloring, every domino covers one black and one white cell.
  • A tiling of all free cells exists if and only if there is a perfect matching in the bipartite graph whose vertices are the free black and white cells, and edges connect orthogonally adjacent free cells.
  • We can test for a perfect matching in O(E√V) using the Hopcroft–Karp algorithm, which is efficient for up to N=40 (V≤1600, E≤4V).

3. Full solution approach based on the observations
  1. Read N and P, and mark the removed cells in a Boolean grid.
  2. Count free cells = N·N − P. If this count is odd, immediately print "No" and exit.
  3. Color each free cell black or white by parity of (x+y). Collect black cells into a list B and white cells into a list W, assigning each a unique index.
  4. Build a bipartite graph: for each black cell at (x,y), check its 4 neighbors (up/down/left/right). If a neighbor is a free white cell, add an edge from the black-index to the white-index.
  5. Run Hopcroft–Karp to find a maximum matching. Let M be the size of the matching. If 2·M ≠ (number of free cells), print "No". Otherwise proceed.
  6. Extract all matched pairs (b,w). For each pair, compare their coordinates:
     – If they share the same x, it's horizontal. Record the leftmost cell.
     – Otherwise it's vertical. Record the bottommost cell.
  7. Print "Yes", then the vertical count and list, followed by the horizontal count and list (this is the order the reference solution prints).

4. C++ implementation
```cpp
#include <bits/stdc++.h>
// #include <coding_library/graph/hopcroft_karp.hpp>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

class HopcroftKarp {
  private:
    int n, m;
    vector<int> dist;

    bool bfs() {
        queue<int> q;
        dist.assign(n, -1);
        for(int u = 0; u < n; u++) {
            if(inv_match[u] == -1) {
                dist[u] = 0;
                q.push(u);
            }
        }

        bool found = false;
        while(!q.empty()) {
            int u = q.front();
            q.pop();
            for(int v: adj[u]) {
                int m = match[v];
                if(m == -1) {
                    found = true;
                } else if(dist[m] == -1) {
                    dist[m] = dist[u] + 1;
                    q.push(m);
                }
            }
        }

        return found;
    }

    bool dfs(int u) {
        for(int v: adj[u]) {
            int m = match[v];
            if(m == -1 || (dist[m] == dist[u] + 1 && dfs(m))) {
                inv_match[u] = v;
                match[v] = u;
                return true;
            }
        }
        dist[u] = -1;
        return false;
    }

  public:
    vector<int> match, inv_match;
    vector<vector<int>> adj;

    HopcroftKarp(int _n, int _m = -1) : n(_n), m(_m == -1 ? _n : _m) {
        adj.assign(n, vector<int>());
        clear(false);
    }

    void clear(bool clear_adj = true) {
        match.assign(m, -1);
        inv_match.assign(n, -1);
        if(clear_adj) {
            adj.assign(n, vector<int>());
        }
    }

    void add_edge(int u, int v) { adj[u].push_back(v); }

    int max_matching(bool shuffle_edges = false) {
        if(shuffle_edges) {
            for(int i = 0; i < n; i++) {
                shuffle(
                    adj[i].begin(), adj[i].end(),
                    mt19937(
                        chrono::steady_clock::now().time_since_epoch().count()
                    )
                );
            }
        }

        int ans = 0;
        while(bfs()) {
            for(int u = 0; u < n; u++) {
                if(inv_match[u] == -1 && dfs(u)) {
                    ans++;
                }
            }
        }
        return ans;
    }

    vector<pair<int, int>> get_matching() {
        vector<pair<int, int>> matches;
        for(int u = 0; u < n; u++) {
            if(inv_match[u] != -1) {
                matches.emplace_back(u, inv_match[u]);
            }
        }
        return matches;
    }

    pair<vector<int>, vector<int>> minimum_vertex_cover() {
        vector<int> left_cover, right_cover;
        bfs();

        for(int u = 0; u < n; u++) {
            if(dist[u] == -1) {
                left_cover.push_back(u);
            }
        }

        for(int v = 0; v < m; v++) {
            if(match[v] != -1 && dist[match[v]] != -1) {
                right_cover.push_back(v);
            }
        }

        return {left_cover, right_cover};
    }
};

using BipartiteMatching = HopcroftKarp;

int n, p;
vector<pair<int, int>> removed;

void read() {
    cin >> n >> p;
    removed.resize(p);
    cin >> removed;
}

void solve() {
    // - Tiling the board with dominoes is a perfect matching on the grid graph
    //   where every remaining cell is a node and edges join orthogonally
    //   adjacent remaining cells. The grid is bipartite by colour (i + j)
    //   parity, so this is bipartite matching.
    //
    // - Build the bipartite graph with black cells on the left and white cells
    //   on the right, add an edge for each adjacent black/white pair, and run
    //   Hopcroft-Karp. A full tiling exists iff the total remaining cell count
    //   is even and the maximum matching covers all of them (2 * matching ==
    //   cells).
    //
    // - For the output, each matched pair is one domino: if the two cells share
    //   a row it is horizontal (report the left cell), otherwise vertical
    //   (report the bottom cell).

    vector<vector<bool>> board(n + 1, vector<bool>(n + 1, true));
    int total_cells = n * n;

    for(auto [x, y]: removed) {
        board[x][y] = false;
        total_cells--;
    }

    if(total_cells % 2 == 1) {
        cout << "No\n";
        return;
    }

    vector<pair<int, int>> black_cells, white_cells;
    map<pair<int, int>, int> black_id, white_id;

    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            if(!board[i][j]) {
                continue;
            }

            if((i + j) % 2 == 0) {
                black_id[{i, j}] = black_cells.size();
                black_cells.push_back({i, j});
            } else {
                white_id[{i, j}] = white_cells.size();
                white_cells.push_back({i, j});
            }
        }
    }

    BipartiteMatching bm(black_cells.size(), white_cells.size());

    int dx[] = {-1, 1, 0, 0};
    int dy[] = {0, 0, -1, 1};

    for(int i = 0; i < (int)black_cells.size(); i++) {
        auto [x, y] = black_cells[i];

        for(int d = 0; d < 4; d++) {
            int nx = x + dx[d];
            int ny = y + dy[d];

            if(nx >= 1 && nx <= n && ny >= 1 && ny <= n && board[nx][ny]) {
                if(white_id.count({nx, ny})) {
                    bm.add_edge(i, white_id[{nx, ny}]);
                }
            }
        }
    }

    int matching = bm.max_matching();
    if(matching * 2 != total_cells) {
        cout << "No\n";
        return;
    }

    cout << "Yes\n";
    vector<pair<int, int>> horizontal, vertical;
    for(int i = 0; i < (int)black_cells.size(); i++) {
        if(bm.inv_match[i] != -1) {
            auto [bx, by] = black_cells[i];
            auto [wx, wy] = white_cells[bm.inv_match[i]];

            if(bx == wx) {
                if(by < wy) {
                    horizontal.push_back({bx, by});
                } else {
                    horizontal.push_back({wx, wy});
                }
            } else {
                if(bx < wx) {
                    vertical.push_back({bx, by});
                } else {
                    vertical.push_back({wx, wy});
                }
            }
        }
    }

    cout << vertical.size() << "\n";
    for(auto [x, y]: vertical) {
        cout << x << " " << y << "\n";
    }

    cout << horizontal.size() << "\n";
    for(auto [x, y]: horizontal) {
        cout << x << " " << y << "\n";
    }
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

5. Python implementation with detailed comments
```python
import sys
from collections import deque

def hopcroft_karp(adj, nL, nR):
    """Return (matching_size, pairU, pairV). adj[u] = list of v's."""
    INF = 10**9
    pairU = [-1]*nL
    pairV = [-1]*nR
    dist  = [0]*nL

    def bfs():
        queue = deque()
        # Initialize free left nodes at distance 0
        for u in range(nL):
            if pairU[u] == -1:
                dist[u] = 0
                queue.append(u)
            else:
                dist[u] = INF
        found = False
        while queue:
            u = queue.popleft()
            for v in adj[u]:
                pu = pairV[v]
                if pu == -1:
                    # Found an unmatched right node
                    found = True
                elif dist[pu] == INF:
                    dist[pu] = dist[u] + 1
                    queue.append(pu)
        return found

    def dfs(u):
        for v in adj[u]:
            pu = pairV[v]
            if pu == -1 or (dist[pu] == dist[u]+1 and dfs(pu)):
                pairU[u] = v
                pairV[v] = u
                return True
        dist[u] = INF
        return False

    matching = 0
    while bfs():
        for u in range(nL):
            if pairU[u] == -1 and dfs(u):
                matching += 1
    return matching, pairU, pairV

def main():
    data = sys.stdin.read().split()
    it = iter(data)
    N, P = map(int, (next(it), next(it)))
    removed = set()
    for _ in range(P):
        x, y = map(int, (next(it), next(it)))
        removed.add((x,y))

    total_free = N*N - P
    # If odd free cells, impossible
    if total_free % 2:
        print("No")
        return

    # Partition free cells by parity
    black, white = [], []
    bidx, widx = {}, {}
    for x in range(1, N+1):
        for y in range(1, N+1):
            if (x,y) in removed: continue
            if ((x+y)&1) == 0:
                bidx[(x,y)] = len(black)
                black.append((x,y))
            else:
                widx[(x,y)] = len(white)
                white.append((x,y))

    # Build adjacency list: black → white
    adj = [[] for _ in range(len(black))]
    for i,(x,y) in enumerate(black):
        for dx,dy in [(1,0),(-1,0),(0,1),(0,-1)]:
            nx, ny = x+dx, y+dy
            if (nx,ny) in widx:
                adj[i].append(widx[(nx,ny)])

    # Run Hopcroft–Karp
    match_sz, pairU, pairV = hopcroft_karp(adj, len(black), len(white))
    if match_sz*2 != total_free:
        print("No")
        return

    print("Yes")
    horiz, vert = [], []
    # Extract domino placements
    for bi, wi in enumerate(pairU):
        if wi < 0: continue
        bx, by = black[bi]
        wx, wy = white[wi]
        if bx == wx:
            # horizontal → pick left cell
            if by < wy:
                horiz.append((bx,by))
            else:
                horiz.append((wx,wy))
        else:
            # vertical → pick bottom cell
            if bx < wx:
                vert.append((bx,by))
            else:
                vert.append((wx,wy))

    # Output vertical then horizontal (matching the reference order)
    print(len(vert))
    for x,y in vert:
        print(x, y)
    print(len(horiz))
    for x,y in horiz:
        print(x, y)

if __name__ == "__main__":
    main()
```
