We'll count the number of ways to choose directions so that the resulting graph would NOT be strong-connected. Main idea: graph is not strong-connected if and only if it has at least one strong-connected subcomponent such that it is a leaf in DAG of components.

Let's call dp[mask] number of ways to choose directions of edges with both ends in mask so that it is not strong-connected.

The answer is 2^m — dp[2^n — 1].

How to count: Let's look at the first vertex in mask. Let's choose it's strong-subcomponent (now it's O(3^n)). Call this subcomponent M1; M2 := mask ^ M1. We suggest that M2 != 0. M1 is either leaf or not. 1) M1 is a leaf. Then we now the direction of all edges between M1 and M2. We need to add 2^{number of edges strictly inside M2} * (2^{number of edges strictly inside M1} — dp[M1]). 2) M1 is not a leaf but M2 is a component and it is a leaf. Add (2^{number of edges strictly inside M2} — dp[M2]) * (2^{number of edges strictly inside M1} — dp[M1]). 3) M1 is not a leaf and so M2 has subcomponent that is a leaf (or the same: M2 is not strong-connected). Add (2^{number of edges between M1 and M2} — 2) * dp[M2] * (2^{number of edges strictly inside M1} — dp[M1]).

Here's my solution.

Let dp[mask] be the number of ways to orient the edges to make each connected component in mask strongly connected. Then answer is dp[2^n-1] (if there is more than one connected component, the answer is zero).

Let's do this by complimentary counting. First, dp[mask] = 2^(edges in mask). We subtract any bad arrangement. For each bad arrangement, we can identify it with the subset of nodes that are in some source strongly connected component (i.e. strongly connected component with indegree 0).

Ideally, we would like to do something like dp[mask] -= dp[submask] * 2^(edges in submask^mask), and then we are done.

What's the problem with this? Well, some bad arrangements are subtracted multiple times. In particular, let's say submask has C connected components. Then, it's actually subtracted 2^C-1 different times right now.

We can fix this by adding a (-1)^(C) term in our subtraction. You can check that -(C choose 1) + (C choose 2) — (C choose 3) + ... (C choose C) = -1, so each bad arrangment is only subtracted exactly once.

The implementation of this is pretty short. After setting everything up, the main dp can be computed in just a couple of lines. Here's an example implementation:

    public int mod = 1000000007;
    public int count(int n, int[] a, int[] b) {
        int m = a.length;
        int[] pow2 = new int[m+1];
        pow2[0] = 1;
        for (int i = 1; i <= m; i++) pow2[i] = (pow2[i-1]<<1) % mod;
        int[] edges = new int[1<<n];
        int[] ncomp = new int[1<<n];
        for (int mask = 0; mask < 1 << n; mask++) {
            ncomp[mask] = Integer.bitCount(mask);
            int[] p = DisjointSets.createSets(n);
            for (int i = 0; i < m; i++) {
                if (((mask>>a[i])&1) == 1 && ((mask>>b[i])&1) == 1) {
                    edges[mask]++;
                    if (DisjointSets.unite(p, a[i], b[i])) {
                        ncomp[mask]--;
                    }
                }
            }
        }
        if (ncomp[(1<<n)-1] > 1) return 0;

        int[] scc = new int[1<<n];
        scc[0] = 1;
        for (int mask = 1; mask < 1 << n; mask++) {
            for (int submask = mask; submask != 0; submask = (submask - 1) & mask) {
                long sign = ncomp[submask^mask] % 2 == ncomp[mask] % 2 ? mod-1 : 1;
                scc[mask] = ((int)((scc[mask] + sign * pow2[edges[submask]] % mod * (long) scc[submask^mask]) % mod));
            }
        }
        return scc[(1<<n)-1];
    }

Found it here.

I solved it in the practice room in the following way.

1. First use divide and conquer to determine an x-coordinate that is within the rectangle: choose the middle x-coordinate, solve recursively for rectangles lying left of this x-coordinate and for rectangles lying right of this x-coordinate. Now all we need to implement is a method for finding rectangles containing a fixed x-coordinate.
2. Fix the lower y-coordinate.
3. Imagine we also fixed the higher y-coordinate. Now for every color it is easy to see we have at most two options: take the closest point from the left of the fixed x-coordinate or take the closest point from the right of the fixed x-coordinate. The answer will then be (hy-ly)*(max_left_dist+max_right_dist), where max_left_dist is the maximum distance from a point chosen from the left and max_right_dist is defined similarly.
4. We can visualize this by plotting for each color (left_dist,right_dist)=(a,b) in the plane. Now we have to choose a point (aopt,bopt) so that for each color either a<=aopt or b<=bopt (or both) and we want to minimize aopt+bopt.
5. From this we can immediately see that we can safely remove points that are completely 'covered' by other points. Therefore for the remaining points as a increases, b decreases.
6. We can also see that the only options for the optimal (aopt,bopt) are all combinations of two remaining adjacent points (if (a1,b1) and (a2,b2) are adjacent, an option for a solution is (a1,b2)).
7. Now instead of fixing the higher y-coordinate, instead we will use a sweepline, maintaining the set of remaining points and the options for the optimal (aopt,bopt) as we go. This sweepline will go from top to bottom, since that way the left_dist and right_dist for each color will only get worse, so we don't have to worry about 'old' values for a color still being in the datastructure.
8. Maintaining these datastructures is a little tricky, you can look at my code to see how I did it (basically use a set with some binary searching for maintaining the points, a priority queue for maintaining the options, and some bookkeeping for determining when an option becomes invalid (the corresponding points are not next to eachother anymore)).
9. Putting it all together, we see that it works, and the complexity seems to be O(n^2*log(n)). (It takes 557ms on the worst test case in the system tests).