<|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.

286. Ancient decoration
Time limit per test: 0.25 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard



The ancient race ofhas built N cities connected with two-way roads. They believed in magic of an even number K and built the roads in such a way that each city had exactly K roads going from it.

The Hollars decided to decorate some roads because of their religious holiday. Because they also believe in magic of an even number 2, each city must have exactly 2 decorated roads going from it.

You have to find the roads needing to be decorated.

Input
The first line of the input contains integers N and K (2≤ KNK/2 lines contains description of one road, being the numbers of the cities connected by this road. The cities are numbered starting from 1. There is no road from a city to itself; each pair of cities is connected by at most one road.

Output
If it is impossible to decorate the roads, the only line of the output must contain 'NO' (without quotes). Otherwise the first line of the output must contain 'YES' (without quotes); the rest of the output must contain N lines, each containing one number of a road to be decorated. The roads are numbered starting from 1 in the same order as they appear in the input.

Example(s)
sample input
sample output
9 4
1 2
1 3
1 4
1 5
2 3
3 4
4 5
5 2
2 6
3 7
4 8
5 9
6 7
7 8
8 9
9 6
7 9
8 6
YES
4
9
5
3
12
13
10
11
15



Novosibirsk SU Contest #2, by Novosibirsk Team #1

<|response|>
1. Abridged problem statement  
You are given an undirected, K-regular graph on N vertices (K is even). There are M = N·K/2 edges, each indexed from 1 to M. Find a subset of exactly N edges so that every vertex is incident to exactly 2 chosen edges. If such a selection exists, output "YES" and any valid list of N edge indices (one per line); otherwise output "NO".

2. Key observations  
• A subgraph where every vertex has degree exactly 2 is called a 2-factor (a disjoint union of cycles covering all vertices).  
• In an even-regular graph one can first orient edges so that each vertex has equal out- and in-degree (i.e.\ out-degree = in-degree = K/2).  
• Once edges are oriented, consider a bipartite graph with left part L = {vertices as sources} and right part R = {vertices as sinks}. Each directed edge u→v becomes a bipartite edge (u in L)–(v in R).  
• This bipartite graph is regular of degree K/2 on both sides. By König's edge-coloring theorem it decomposes into exactly K/2 perfect matchings. Any one matching gives each vertex exactly one outgoing and one incoming matched edge ⇒ total degree 2 in the original undirected graph. Such a decomposition always exists, so the answer is always "YES".

3. Full solution approach
Step A: Read input and build adjacency lists with edge indices.
Step B: Orient edges with a single Hierholzer-style DFS (`orient_graph`).
  - Maintain for each vertex its list of incident (unused) edges.
  - Walk along unused edges marking them used; when a recursive call returns, record the edge oriented in the direction it was first traversed.
  - In the end, we obtain a list of directed edges (u→v) each annotated by its original index (1…M). Each vertex has out-degree = in-degree = K/2.
Step C: Build the K/2-regular bipartite graph on 2N nodes: left side L = {0…N−1}, right side R = {0…N−1}; each directed edge (u→v, idx) becomes a bipartite edge from u∈L to v∈R.
Step D: Edge-color that bipartite graph into K/2 perfect matchings (`BipartiteColoring::euler_colour`). The maximum degree is halved recursively by splitting each Euler tour's edges alternately into two subgraphs; whenever the current maximum degree is odd, one matching is first peeled off with Hopcroft–Karp (`slow_one_colour`) to make the degree even. Each resulting color class is a perfect matching of size N.
Step E: Take the first color class `ans[0]`, collect the original indices of its edges, sort them ascending, print "YES" followed by those N indices.

Time complexity:
  - Building adjacency and orienting edges takes O(M + N).
  - The coloring has O(log K) recursion levels; Hopcroft–Karp (O(E·√V)) is called only O(log K) times, which is efficient for typical contest limits.

4. C++ Solution
```cpp
#include <bits/stdc++.h>

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;
}

template<class T>
void make_larger_if_needed(vector<T>& v, int size) {
    if(v.size() < size) {
        v.resize(size);
    }
}

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;
    }
};

struct Edge {
    int u, v, idx;
    Edge(int _u = 0, int _v = 0, int _idx = 0) : u(_u), v(_v), idx(_idx) {}
};

class BipartiteColoring {
  private:
    vector<vector<pair<Edge, int>>> edges_for_ver;
    vector<bool> used;
    vector<vector<pair<int, int>>> adj;
    vector<int> memory_m;

    vector<Edge> slow_one_colour(
        const vector<Edge>& edges, int n, vector<vector<Edge>>& answer
    ) {
        make_larger_if_needed(edges_for_ver, n);
        make_larger_if_needed(used, (int)edges.size());

        for(int i = 0; i < n; i++) {
            edges_for_ver[i].clear();
        }
        for(int i = 0; i < (int)edges.size(); i++) {
            used[i] = false;
        }

        HopcroftKarp bm(n, n);
        for(int i = 0; i < (int)edges.size(); i++) {
            auto e = edges[i];
            bm.add_edge(e.u, e.v);
            edges_for_ver[e.u].push_back({e, i});
        }

        int max_match = bm.max_matching();
        assert(max_match == n);

        vector<Edge> assigned;
        vector<pair<int, int>> matches = bm.get_matching();
        for(auto [u, v]: matches) {
            for(auto [e, ei]: edges_for_ver[u]) {
                if(e.v == v && !used[ei]) {
                    used[ei] = true;
                    assigned.push_back(e);
                    break;
                }
            }
        }

        vector<Edge> new_edges;
        for(int i = 0; i < (int)edges.size(); i++) {
            if(!used[i]) {
                new_edges.push_back(edges[i]);
            }
        }

        answer.push_back(assigned);
        return new_edges;
    }

    pair<vector<Edge>, vector<Edge>> partition_edges_euler(
        const vector<Edge>& edges, const vector<int>& vers, int n, int m
    ) {
        make_larger_if_needed(adj, 2 * n);
        make_larger_if_needed(memory_m, m);

        for(int v: vers) {
            adj[v].clear();
            adj[v + n].clear();
        }

        for(int ei = 0; ei < (int)edges.size(); ei++) {
            auto e = edges[ei];
            adj[e.u].push_back({e.v + n, ei});
            adj[e.v + n].push_back({e.u, ei});
            memory_m[ei] = 0;
        }

        function<void(int, vector<Edge>&, vector<Edge>&)> dfs =
            [&](int v, vector<Edge>& subgraph_0, vector<Edge>& subgraph_1) {
                while(!adj[v].empty()) {
                    auto [u, ei] = adj[v].back();
                    adj[v].pop_back();
                    if(memory_m[ei] == 1) {
                        continue;
                    }
                    memory_m[ei] = 1;
                    dfs(u, subgraph_0, subgraph_1);
                    if(v < n) {
                        subgraph_0.push_back(edges[ei]);
                    } else {
                        subgraph_1.push_back(edges[ei]);
                    }
                    break;
                }
            };

        vector<Edge> subgraph_0, subgraph_1;
        for(int v: vers) {
            while(!adj[v].empty()) {
                dfs(v, subgraph_0, subgraph_1);
            }
        }

        return {subgraph_0, subgraph_1};
    }

  public:
    int euler_colour(
        const vector<Edge>& edges, int n, int m, vector<vector<Edge>>& answer
    ) {
        static vector<int> memory;
        make_larger_if_needed(memory, n);

        vector<int> vers, _vers;
        for(auto e: edges) {
            _vers.push_back(e.u);
            _vers.push_back(e.v);
        }

        int max_degree = 0;
        for(int v: _vers) {
            memory[v] = -1;
        }
        for(int v: _vers) {
            if(memory[v] == -1) {
                vers.push_back(v);
                memory[v] = 0;
            }
        }
        for(auto e: edges) {
            memory[e.u]++;
            max_degree = max(max_degree, memory[e.u]);
        }

        if(max_degree == 0) {
            return 0;
        }
        if(max_degree == 1) {
            answer.push_back({});
            for(auto e: edges) {
                answer.back().push_back(e);
            }
            return 1;
        }

        if(max_degree % 2 == 1) {
            auto subgraph = slow_one_colour(edges, n, answer);
            return 1 + euler_colour(subgraph, n, m, answer);
        }

        auto [subgraph_0, subgraph_1] =
            partition_edges_euler(edges, vers, n, m);
        int colour_num_subgraph_0 = euler_colour(subgraph_0, n, m, answer);

        int d = max_degree, q = 0;
        while((1 << q) < (max_degree / 2)) {
            q++;
        }
        int to_remove_count = (1 << q) - (max_degree / 2);
        if(to_remove_count > 0 && colour_num_subgraph_0 >= to_remove_count) {
            for(int i = answer.size() - 1; i >= answer.size() - to_remove_count;
                i--) {
                for(auto& e: answer[i]) {
                    subgraph_1.push_back(e);
                }
            }
            answer.erase(answer.end() - to_remove_count, answer.end());
        }

        int colour_num_subgraph_1 = euler_colour(subgraph_1, n, m, answer);
        return colour_num_subgraph_0 + colour_num_subgraph_1;
    }
};

vector<Edge> orient_graph(
    int n, int m, const vector<int>& a, const vector<int>& b
) {
    vector<vector<pair<int, int>>> adj(n);
    vector<Edge> edges;
    for(int i = 0; i < m; i++) {
        int u = a[i], v = b[i];
        adj[u].push_back({v, i});
        adj[v].push_back({u, i});
        edges.emplace_back(u, v, i + 1);
    }

    vector<bool> used(m, false);
    vector<Edge> directed_edges;
    function<void(int)> dfs = [&](int v) {
        while(!adj[v].empty()) {
            auto [u, ei] = adj[v].back();
            adj[v].pop_back();
            if(used[ei]) {
                continue;
            }
            used[ei] = true;
            dfs(u);
            directed_edges.push_back({v, u, ei + 1});
            break;
        }
    };

    for(int v = 0; v < n; v++) {
        while(!adj[v].empty()) {
            dfs(v);
            break;
        }
    }

    return directed_edges;
}

int n, k, m;
vector<int> a, b;

void read() {
    cin >> n >> k;
    m = n * k / 2;
    a.resize(m);
    b.resize(m);
    for(int i = 0; i < m; i++) {
        cin >> a[i] >> b[i];
        a[i]--;
        b[i]--;
    }
}

void solve() {
    // The K-regular graph (K even) must be split so each city keeps exactly 2
    // decorated roads, i.e. extract a 2-regular spanning subgraph (one perfect
    // matching's worth in the oriented bipartite model).
    //
    // - Eulerian-orient the multigraph so every vertex has out-degree K/2; each
    //   directed edge u->v becomes a bipartite edge between left u and right v,
    //   giving a K/2-regular bipartite (multi)graph.
    //
    // - Edge-colour that bipartite graph into K/2 perfect matchings via Euler
    //   partitioning: repeatedly split the edge set along Euler tours to halve
    //   the degree, peeling one matching with Hopcroft-Karp whenever the max
    //   degree is odd. Each resulting colour class is a perfect matching of size
    //   n, so picking any one class gives every city exactly 2 incident roads.
    //
    // - Output the original road indices of the first colour class, sorted. The
    //   construction always succeeds for an even-regular graph, so the answer is
    //   always YES.

    vector<Edge> edges = orient_graph(n, m, a, b);

    vector<vector<Edge>> ans;
    BipartiteColoring bc;
    bc.euler_colour(edges, n, m, ans);

    cout << "YES" << "\n";
    vector<int> roads;
    for(auto& e: ans[0]) {
        roads.push_back(e.idx);
    }

    sort(roads.begin(), roads.end());
    for(int i = 0; i < n; i++) {
        cout << roads[i] << "\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
This Python reference implements the simpler single-perfect-matching variant (orient, then take one matching of the bipartite graph as the 2-factor); the C++ above instead computes the full K/2-color decomposition and emits its first class. Both produce a valid 2-factor.
```python
import sys
from collections import deque

sys.setrecursionlimit(10**7)
input = sys.stdin.readline

# Hopcroft–Karp on bipartite graph
class HopcroftKarp:
    def __init__(self, n, m):
        self.n, self.m = n, m
        self.adj = [[] for _ in range(n)]
        self.dist = [0]*n
        self.matchL = [-1]*n
        self.matchR = [-1]*m

    def add_edge(self, u, v):
        self.adj[u].append(v)

    def bfs(self):
        q = deque()
        for u in range(self.n):
            if self.matchL[u] < 0:
                self.dist[u] = 0
                q.append(u)
            else:
                self.dist[u] = -1
        foundAug = False
        while q:
            u = q.popleft()
            for v in self.adj[u]:
                u2 = self.matchR[v]
                if u2 < 0:
                    foundAug = True
                elif self.dist[u2] < 0:
                    self.dist[u2] = self.dist[u] + 1
                    q.append(u2)
        return foundAug

    def dfs(self, u):
        for v in self.adj[u]:
            u2 = self.matchR[v]
            if u2 < 0 or (self.dist[u2] == self.dist[u] + 1 and self.dfs(u2)):
                self.matchL[u] = v
                self.matchR[v] = u
                return True
        self.dist[u] = -1
        return False

    def max_matching(self):
        res = 0
        while self.bfs():
            for u in range(self.n):
                if self.matchL[u] < 0 and self.dfs(u):
                    res += 1
        return res

def main():
    N, K = map(int, input().split())
    M = N * K // 2
    A = [0]*M
    B = [0]*M
    for i in range(M):
        a, b = map(int, input().split())
        A[i], B[i] = a-1, b-1

    # Build adjacency with edge indices
    adj = [[] for _ in range(N)]
    for i in range(M):
        u, v = A[i], B[i]
        adj[u].append((v, i))
        adj[v].append((u, i))

    # Hierholzer to orient edges
    used = [False]*M
    directed = []  # will store tuples (u, v, edge_index)
    for start in range(N):
        if not adj[start]:
            continue
        stackV = [start]
        stackEi = []
        while stackV:
            v = stackV[-1]
            # skip used edges
            while adj[v] and used[adj[v][-1][1]]:
                adj[v].pop()
            if not adj[v]:
                stackV.pop()
                if stackEi:
                    ei = stackEi.pop()
                    u = stackV[-1]
                    directed.append((u, v, ei))
            else:
                u, ei = adj[v].pop()
                if used[ei]: continue
                used[ei] = True
                stackV.append(u)
                stackEi.append(ei)

    # Build bipartite graph on N + N
    hk = HopcroftKarp(N, N)
    outEdges = [[] for _ in range(N)]
    for u, v, ei in directed:
        hk.add_edge(u, v)
        outEdges[u].append((v, ei+1))  # store 1-based index

    # Find matching
    if hk.max_matching() < N:
        print("NO")
        return

    print("YES")
    answer = [0]*N
    for u in range(N):
        v = hk.matchL[u]
        # find stored (v, idx)
        for vv, idx in outEdges[u]:
            if vv == v:
                answer[u] = idx
                break

    # Print one per line
    print("\n".join(map(str, answer)))

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