Problem: There are 2N balls in the xy-plane. The coordinates of the i-th of them is (x_i, y_i). Here, x_i and y_i are integers between 1 and N (inclusive) for all i, and no two balls occupy the same coordinates.  In order to collect these balls, Snuke prepared 2N robots, N of type A and N of type B. Then, he placed the type-A robots at coordinates (1, 0), (2, 0), ..., (N, 0), and the type-B robots at coordinates (0, 1), (0, 2), ..., (0, N), one at each position.  When activated, each type of robot will operate as follows.    * When a type-A robot is activated at coordinates (a, 0), it will move to the position of the ball with the lowest y-coordinate among the balls on the line x = a, collect the ball and deactivate itself. If there is no such ball, it will just deactivate itself without doing anything.    * When a type-B robot is activated at coordinates (0, b), it will move to the position of the ball with the lowest x-coordinate among the balls on the line y = b, collect the ball and deactivate itself. If there is no such ball, it will just deactivate itself without doing anything.  Once deactivated, a robot cannot be activated again. Also, while a robot is operating, no new robot can be activated until the operating robot is deactivated.  When Snuke was about to activate a robot, he noticed that he may fail to collect all the balls, depending on the order of activating the robots.  Among the (2N)! possible orders of activating the robots, find the number of the ones such that all the balls can be collected, modulo 1 000 000 007.

Submission attempts (in order):
import sys
input = sys.stdin.readline

from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import connected_components

MOD = 10**9 + 7

N = int(input())
ball = (tuple(int(x) for x in row.split()) for row in sys.stdin.readlines())




graph = [set() for _ in range(N+N+1)]
for x,y in ball:
    graph[x].add(y+N)
    graph[y+N].add(x)

visited = [False] * (N+N+1)
components = []
for x in range(1,N+N+1):
    if visited[x]:
        continue
    V = set([x])
    E = []
    q = [x]
    visited[x] = True
    while q:
        y = q.pop()
        for z in graph[y]:
            if y < z:
                E.append((y,z))
            if visited[z]:
                continue
            V.add(z)
            visited[z] = True
            q.append(z)
    components.append((V,E))


import sys
input = sys.stdin.readline

MOD = 10**9 + 7

N = int(input())
ball = (tuple(int(x) for x in row.split()) for row in sys.stdin.readlines())




graph = [set() for _ in range(N+N+1)]
for x,y in ball:
    graph[x].add(y+N)
    graph[y+N].add(x)

visited = [False] * (N+N+1)
components = []
for x in range(1,N+N+1):
    if visited[x]:
        continue
    V = set([x])
    E = []
    q = [x]
    visited[x] = True
    while q:
        y = q.pop()
        for z in graph[y]:
            if y < z:
                E.append((y,z))
            if visited[z]:
                continue
            V.add(z)
            visited[z] = True
            q.append(z)
    components.append((V,E))

def make_get_pattern(V):
    deg1 = [x for x in V if len(graph[x]) == 1]
    get = {}
    while deg1:
        x = deg1.pop()
        if not graph[x]:
            continue
        y = graph[x].pop()
        se = graph[y]; se.remove(x)
        if len(se) == 1: deg1.append(y)
        if x < y:
            get[(x,y)] = 0
        else:
            get[(y,x)] = 1
    for x in V:
        if graph[x]:
            y = graph[x].pop()
            break
    
    graph[y].remove(x)
    if x > y: x,y = y,x
    get[(x,y)] = 2
    while graph[x]:
        y = graph[x].pop()
        graph[y].remove(x)
        if x < y:
            get[(x,y)] = 3
        else:
            get[(y,x)] = 2
        x = y
    return get

def F(V,E):
    # V is connected
    if len(E) != len(V):
        return 0
    ret = 0
    E.sort()
    get = make_get_pattern(V)
    den1,den2 = 1,1
    dp1 = {x:0 for x in V}
    dp2 = {x:0 for x in V}
    for x,y in E:
        if get[(x,y)] == 0:
            k1 = dp1[x] + 1; k2 = dp2[x] + 1
        elif get[(x,y)] == 1:
            k1 = dp1[y] + 1; k2 = dp2[y] + 1
        elif get[(x,y)] == 2:
            k1 = dp1[x] + 1; k2 = dp2[y] + 1
        else:
            k1 = dp1[y] + 1; k2 = dp2[x] + 1
        dp1[x] += k1; dp1[y] += k1
        dp2[x] += k2; dp2[y] += k2
        den1 *= k1; den2 *= k2
        den1 %= MOD; den2 %= MOD
    return sum(pow(x,MOD-2,MOD) for x in (den1,den2))

prob = 1
for c in components:
    prob *= F(*c)
    prob %= MOD

answer = prob
for n in range(1,N+N+1):
    answer *= n
    answer %= MOD
print(answer)

Verdicts: Wrong Answer

Accepted