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

525. Revolutionary Roads
Time limit per test: 1 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard

Governments of different countries like to boast about their achievements. For instance, the President of Flatland has announced that his country has the most advanced road system. He said the degree of a country road system development is equal to the amount of cities in the largest  subset of cities. A subset of cities is called  if it is possible to get from any city of the subset to all other cities of the subset.

Not to lag behind the neighbors Berland's President decided to undertake a reform and modernize roads in his country. All the roads in Berland are one-way, each of them connects a pair of cities in one direction. There is at most one road in each direction between any two given cities.

Since there is little money in the budget, President's plans aren't very ambitious. He can turn at most one of all given one-way roads into a two-way road. And he wants to do it in such a way that the resulting road system degree of development in Berland becomes as high as possible. Let's say the maximum degree of development, which can be achieved by this action, is equal to w.

A road is called  if, after it is changed from one-way to two-way, the degree of road system development becomes equal to w. Your task is to find all  roads.

Input
The first line of input contains a pair of numbers n, m (1 ≤ n ≤ 1000, 0 ≤ m ≤ 20000), where n — the number cities, m — the number of roads. The following m lines contain descriptions of the roads. Each line contains a pair of integers ai, bi (1 ≤ ai, bi ≤ n, ai ≠ bi), representing a one-way road from city ai to city bi. Cities are numbered from 1 to n.

Output
Write w to the first line of output. To the second line write t — number of roads in the required subset. To the third line write indices of the roads in this subset. Roads are numbered from 1 to m according to their order in the input file.

Example(s)
sample input
sample output
5 4
1 2
2 3
1 3
4 1
3
1
3

sample input
sample output
3 4
1 2
2 1
1 3
3 1
3
4
1 2 3 4

<|response|>
1. Abridged Problem Statement
You are given a directed graph of n cities (nodes) and m one-way roads (edges). You may choose at most one existing directed edge u→v and add the reverse edge v→u, effectively making that road two-way. After this change, the "degree of development" is the size (number of cities) of the largest strongly connected component (SCC) in the new graph.
• Compute the maximum possible SCC size w achievable by upgrading at most one road.
• List all edges (by their input indices) whose upgrading attains this maximum w.

2. Key Observations
- In the original graph, find its SCC decomposition. Contract each SCC into a single node to form a condensation DAG.
- Adding the reverse of an edge u→v can only merge SCCs that lie on some path cu→…→cv in the condensation (where cu=comp[u], cv=comp[v]).
- A component i merges into the new big SCC if and only if cu can reach i (closure[cu][i] is set) AND i can reach cv (closure[i][cv] is set) — using only the forward reachability closure.
- We need to know, for every component x, which components it can reach (forward closure), and we can then check both conditions using only the single forward closure array.
- For each edge, merging gives total size = sum of comp_size over the merged set. The resulting largest SCC is max(merged_size, largest comp_size among components not in that merged set).
- Edge cases: if cu==cv, upgrading does nothing and the best SCC size is comp_size[cu]; if m=0, no edges exist, best w=1 (any single city).

3. Full Solution Approach
Step 1: Read n, m and the list of directed edges (0-based) in a vector edges.
Step 2: Run SCC (Kosaraju) on the original graph to get for each vertex v its component id comp[v], and count C = number of components. Also compute comp_size[c] = number of original vertices in component c.
Step 3: Build the condensation DAG of size C: for each original edge u→v, if comp[u]≠comp[v], add a single DAG edge cu→cv (deduplicated), and record indegrees.
Step 4: Topologically sort the DAG (Kahn's algorithm), obtaining an order topo[].
Step 5: Compute forward reachability closure: create a bitset<1024> closure[c] of size C, initialize closure[c][c]=1. Process components in reverse topological order; for each u, for each v in adj[u], do closure[u] |= closure[v].
Step 6: Now iterate over each original edge i: let u,v be its endpoints and cu=comp[u], cv=comp[v].
   - If cu==cv: current_max = comp_size[cu].
   - Else: for each component w, it merges if closure[cu][w] && closure[w][cv].
        merged_size = sum(comp_size[w] for all such w)
        other_max = maximum comp_size[w] for the rest
        current_max = max(merged_size, other_max)
   Track the global maximum w and maintain a list of all edge indices i (1-based) for which current_max == w.
Step 7: If m==0, set w=1 and answer-set empty.
Step 8: Output w, the count of chosen edges, then the sorted list of their indices.

4. C++ Solution
```cpp
#include <bits/stdc++.h>
using namespace std;

class StronglyConnectedComponents {
  private:
    vector<bool> visited;

    void dfs1(int u) {
        visited[u] = true;
        for(int v: adj[u]) {
            if(!visited[v]) {
                dfs1(v);
            }
        }
        top_sort.push_back(u);
    }

    void dfs2(int u) {
        for(int v: radj[u]) {
            if(comp[v] == -1) {
                comp[v] = comp[u];
                dfs2(v);
            }
        }
    }

  public:
    int n;
    vector<vector<int>> adj, radj;
    vector<int> comp, comp_ids, top_sort;

    StronglyConnectedComponents() {}
    StronglyConnectedComponents(int _n) { init(_n); }

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

    void init(int _n) {
        n = _n;
        comp_ids.clear();
        top_sort.clear();
        adj.assign(n, {});
        radj.assign(n, {});
    }

    void find_components() {
        comp.assign(n, -1);
        visited.assign(n, false);

        for(int i = 0; i < n; i++) {
            if(!visited[i]) {
                dfs1(i);
            }
        }

        reverse(top_sort.begin(), top_sort.end());
        for(int u: top_sort) {
            if(comp[u] == -1) {
                comp[u] = (int)comp_ids.size();
                comp_ids.push_back(comp[u]);
                dfs2(u);
            }
        }
    }
};

int n, m;
vector<pair<int, int>> edges;
StronglyConnectedComponents G;

void read() {
    cin >> n >> m;
    edges.resize(m);
    G.init(n);
    for(int i = 0; i < m; i++) {
        int u, v;
        cin >> u >> v;
        u--, v--;
        edges[i] = {u, v};
        G.add_edge(u, v);
    }
}

void solve() {
    // Both N and M aren't huge, and the naive solution is O(M^2) with trying to
    // add every edge and then finding the largest strongly connected component.
    // However, to have something acceptable we would have to implement a O(NM)
    // approach. This is possible by a slight modification - we will only
    // calculate the strongly connected components decomposition in the
    // beginning in O(M), and then calculate the transitive closure between the
    // components in O(NM), or even O(NM/w) with bitset (word size w = 32 or
    // 64). We will then try every edge (u, v) and if u and v are in different
    // components, we will add an edge v -> u (note u -> v are already an edge),
    // we will merge comp[u], comp[v], and all w such that C[comp[u], w] and
    // C[w, comp[v]] (here C[x, y] means y is reachable from x or essentially
    // this is the transitive closure). This check can be done in O(N) by simply
    // iterating through all w, which is quicker than O(M).

    G.find_components();

    int num_components = G.comp_ids.size();

    vector<vector<bool>> used(
        num_components, vector<bool>(num_components, false)
    );
    vector<vector<int>> comp_adj(num_components);
    vector<int> indegree(num_components, 0);

    for(int i = 0; i < m; i++) {
        int u = edges[i].first, v = edges[i].second;
        int cu = G.comp[u], cv = G.comp[v];
        if(cu != cv && !used[cu][cv]) {
            used[cu][cv] = true;
            comp_adj[cu].push_back(cv);
            indegree[cv]++;
        }
    }

    vector<int> comp_size(num_components, 0);
    for(int i = 0; i < n; i++) {
        comp_size[G.comp[i]]++;
    }

    queue<int> q;
    vector<int> topo_order;
    for(int i = 0; i < num_components; i++) {
        if(indegree[i] == 0) {
            q.push(i);
        }
    }

    while(!q.empty()) {
        int u = q.front();
        q.pop();
        topo_order.push_back(u);

        for(int v: comp_adj[u]) {
            indegree[v]--;
            if(indegree[v] == 0) {
                q.push(v);
            }
        }
    }

    vector<bitset<1024>> closure(num_components);
    for(int i = 0; i < num_components; i++) {
        closure[i][i] = 1;
    }
    for(int i = num_components - 1; i >= 0; i--) {
        int u = topo_order[i];
        for(int v: comp_adj[u]) {
            closure[u] |= closure[v];
        }
    }

    int max_size = 0;
    vector<int> good_edges;

    for(int edge_idx = 0; edge_idx < m; edge_idx++) {
        int u = edges[edge_idx].first, v = edges[edge_idx].second;
        int cu = G.comp[u], cv = G.comp[v];

        vector<bool> will_merge(num_components, false);
        for(int i = 0; i < num_components; i++) {
            if(closure[cu][i] && closure[i][cv]) {
                will_merge[i] = true;
            }
        }

        vector<int> new_comp_size(num_components, 0);
        int merged_size = 0;

        for(int i = 0; i < n; i++) {
            int comp_id = G.comp[i];
            if(will_merge[comp_id]) {
                merged_size++;
            } else {
                new_comp_size[comp_id]++;
            }
        }

        int current_max = merged_size;
        for(int i = 0; i < num_components; i++) {
            current_max = max(current_max, new_comp_size[i]);
        }

        if(current_max > max_size) {
            max_size = current_max;
            good_edges = {edge_idx + 1};
        } else if(current_max == max_size) {
            good_edges.push_back(edge_idx + 1);
        }
    }

    if(m == 0) {
        max_size = 1;
        good_edges = {};
    }

    cout << max_size << endl;
    cout << good_edges.size() << endl;
    for(int i = 0; i < (int)good_edges.size(); i++) {
        if(i > 0) {
            cout << " ";
        }
        cout << good_edges[i];
    }
    if(!good_edges.empty()) {
        cout << endl;
    }
}

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

    int T = 1;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}
```

5. Python Implementation with Detailed Comments
```python
import sys
sys.setrecursionlimit(10**7)

def readints():
    return map(int, sys.stdin.readline().split())

# Step 1: read input
n, m = map(int, sys.stdin.readline().split())
edges = []
g = [[] for _ in range(n)]
rg = [[] for _ in range(n)]
for _ in range(m):
    u, v = map(int, sys.stdin.readline().split())
    u -= 1; v -= 1
    edges.append((u, v))
    g[u].append(v)
    rg[v].append(u)

# Step 2: Kosaraju to find SCCs
visited = [False]*n
order = []

def dfs1(u):
    visited[u] = True
    for w in g[u]:
        if not visited[w]:
            dfs1(w)
    order.append(u)

for i in range(n):
    if not visited[i]:
        dfs1(i)

comp = [-1]*n
cid = 0
def dfs2(u, cid):
    comp[u] = cid
    for w in rg[u]:
        if comp[w] == -1:
            dfs2(w, cid)

for u in reversed(order):
    if comp[u] == -1:
        dfs2(u, cid)
        cid += 1

C = cid

# Step 2B: component sizes
comp_size = [0]*C
for v in range(n):
    comp_size[comp[v]] += 1

# Step 3: build condensation DAG and its reverse
cadj = [set() for _ in range(C)]
indeg = [0]*C
for u, v in edges:
    cu, cv = comp[u], comp[v]
    if cu != cv and cv not in cadj[cu]:
        cadj[cu].add(cv)
        indeg[cv] += 1

# Step 4: topological sort (Kahn's algorithm)
from collections import deque
q = deque([i for i in range(C) if indeg[i] == 0])
topo = []
while q:
    u = q.popleft()
    topo.append(u)
    for v in cadj[u]:
        indeg[v] -= 1
        if indeg[v] == 0:
            q.append(v)

# Step 5: compute forward closure using ints as bitsets
closure = [0]*C
for i in range(C):
    closure[i] = 1 << i

# forward: process reverse topo
for u in reversed(topo):
    for v in cadj[u]:
        closure[u] |= closure[v]

# Step 6: evaluate each edge
best = 0
ans = []
if m == 0:
    best = 1
else:
    for idx, (u, v) in enumerate(edges, start=1):
        cu, cv = comp[u], comp[v]
        if cu == cv:
            cur = comp_size[cu]
        else:
            merged_sz = 0
            other_max = 0
            # component i merges if closure[cu] has bit i AND closure[i] has bit cv
            for c in range(C):
                if (closure[cu] >> c) & 1 and (closure[c] >> cv) & 1:
                    merged_sz += comp_size[c]
                else:
                    other_max = max(other_max, comp_size[c])
            cur = max(merged_sz, other_max)
        if cur > best:
            best = cur
            ans = [idx]
        elif cur == best:
            ans.append(idx)

# Step 7: output
print(best)
print(len(ans))
if ans:
    print(*ans)
```

Explanation Comments in Code
- We first decompose the original graph into SCCs in O(n+m).
- We build a smaller DAG of components and topologically sort it in O(C + #dag_edges).
- We compute the forward reachability closure in O(C²/word_size) using bitsets.
- Finally, for each of the m edges we compute the effect of adding its reverse in O(C) by checking both closure[cu][i] and closure[i][cv] for every component i.
This runs comfortably within the constraints (n≤1000, m≤20000) under 1 second.
