1. Abridged problem statement
You are given an undirected K-regular graph on N vertices (K is even) with M = N·K/2 edges. You must select exactly N edges so that every vertex is incident to exactly 2 selected edges. Print “YES” and list the chosen edge indices (1-based). For an even-regular graph such a selection always exists, so the answer is always “YES.”

2. Detailed editorial
We need a 2-factor (every vertex of degree 2) in an even-regular graph. The solution proceeds in two steps.

Step A: Orient each undirected edge so that every vertex has equal in- and out-degree (hence out-degree = in-degree = K/2).
Because K is even, every component of the graph is Eulerian. We can find, separately in each component, an Euler circuit and orient its edges along the traversal direction. In an Euler circuit each time you enter a vertex you later leave it, so in- and out-counts match. After processing all components we get a directed graph where every vertex has out-degree K/2. The code does this with a recursive Hierholzer-style DFS (`orient_graph`) that pops unused incident edges and records each edge oriented in the direction it was first traversed.

Step B: From this directed K/2-regular digraph on N vertices, build a bipartite graph of size N+N:
  – Left part L = {0..N−1} represents “sources” (the tail of each directed edge).
  – Right part R = {0..N−1} represents “sinks” (the head).
  – For every directed edge u→v we add a bipartite edge (u in L)–(v in R).
This is a K/2-regular balanced bipartite graph. By König's edge-coloring theorem a bipartite graph of maximum degree Δ decomposes into exactly Δ perfect matchings, so here the edges split into K/2 perfect matchings (color classes). Any single perfect matching gives every L-vertex exactly one outgoing matched edge and every R-vertex exactly one incoming matched edge, i.e. degree 2 in the original undirected sense — exactly a 2-factor. So we only need one color class, and since a decomposition always exists, the answer is always “YES.”

How the edge coloring is produced (`BipartiteColoring::euler_colour`): rather than running K/2 separate matchings, the code halves the maximum degree recursively. If the current maximum degree is 1, the whole edge set is already a matching. If it is even, an Euler partition splits every Euler tour's edges alternately into two subgraphs of half the degree, and each half is colored recursively; a small power-of-two bookkeeping step (`to_remove_count`) moves a few finished color classes between the two halves so the recursion balances out. If the maximum degree is odd, one perfect matching is peeled off first with Hopcroft–Karp (`slow_one_colour`) to make the remaining degree even, then the even case applies. The result is K/2 color classes, each a perfect matching of size N.

Implementation details:
– `orient_graph` orients all M edges via a single Hierholzer DFS, marking edges used as it goes.
– `slow_one_colour` extracts one perfect matching using Hopcroft–Karp (O(E√V)); it is invoked only when the current degree is odd, so it runs O(log K) times overall.
– The Euler partition routine runs in linear time per level, and there are O(log K) recursion levels, so the whole coloring is efficient.
– Finally we output the original edge indices of the first color class `ans[0]`, sorted ascending.

3. 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;
}
```

4. Python solution with detailed comments
This Python reference implements the simpler single-perfect-matching variant from Step B (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
sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()

# Hopcroft–Karp for bipartite perfect matching
from collections import deque
class HopcroftKarp:
    def __init__(self,n,m):
        self.n, self.m = n, m
        self.adj = [[] for _ in range(n)]
    def add_edge(self,u,v):
        self.adj[u].append(v)
    def max_matching(self):
        # -1 means unmatched
        matchL = [-1]*self.n
        matchR = [-1]*self.m
        dist = [0]*self.n
        INF = 10**9

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

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

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

def main():
    data = sys.stdin.read().split()
    it = iter(data)
    N = int(next(it))
    K = int(next(it))
    M = N*K//2
    A = [0]*M
    B = [0]*M
    for i in range(M):
        A[i] = int(next(it))-1
        B[i] = int(next(it))-1

    # Build adjacency for undirected graph
    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))

    # Prepare to find Euler tours and orient edges
    used = [False]*M
    directed = []  # list of (u,v,index+1)
    stack = []
    # We'll reuse a stack for each component
    for start in range(N):
        if not adj[start]:
            continue
        stack = [start]
        path = []    # will hold the Euler circuit in vertex order
        edge_stack = []  # parallel stack to record which edge we took
        while stack:
            v = stack[-1]
            # remove used edges
            while adj[v] and used[adj[v][-1][1]]:
                adj[v].pop()
            if not adj[v]:
                stack.pop()
                path.append(v)
                if edge_stack:
                    ei, pu = edge_stack.pop()
                    # orient from pu->v
                    directed.append((pu,v,ei+1))
            else:
                u,ei = adj[v].pop()
                if used[ei]:
                    continue
                used[ei] = True
                stack.append(u)
                edge_stack.append((ei,v))

    # Now build bipartite graph for the directed digraph
    hk = HopcroftKarp(N,N)
    # We also need to remember indices for each (u->v)
    edge_map = [[] for _ in range(N)]
    for u,v,ei in directed:
        hk.add_edge(u,v)
        edge_map[u].append((v,ei))

    if hk.max_matching() < N:
        print("NO")
        return
    print("YES")
    # extract one matched edge per u
    res = [0]*N
    for u in range(N):
        v = hk.matchL[u]
        # pick the stored edge-index that goes to v
        for vv,ei in edge_map[u]:
            if vv==v:
                res[u] = ei
                break
    # print each on its own line
    out = "\n".join(str(x) for x in res)
    print(out)

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

5. Compressed editorial
Find a 2-factor in an even-regular graph by:
a) Euler-orient all edges so each vertex has out-degree = in-degree = K/2.
b) Form a bipartite graph L,R of size N each, where each directed edge u→v becomes (u∈L)–(v∈R), a K/2-regular bipartite graph.
c) Edge-color that bipartite graph into K/2 perfect matchings by recursively halving the maximum degree along Euler partitions (peeling one Hopcroft–Karp matching whenever the degree is odd). Output any one color class: each vertex then has exactly one incoming and one outgoing matched edge ⇒ degree 2. Such a decomposition always exists, so the answer is always YES.