Jellyfish and EVA (Codeforces 1874C)
Tags: Math, probabilities

Monsters have invaded the town again! Asuka invites her good friend, Jellyfish, to drive EVA with her.There are $$$n$$$ cities in the town. All the monsters are in city $$$n$$$. Jellyfish and Asuka are currently in city $$$1$$$ and need to move to city $$$n$$$ to defeat the monsters.There are $$$m$$$ roads. The $$$i$$$-th road allows one to travel from city $$$a_i$$$ to city $$$b_i$$$. All the roads are directed. That is, one cannot travel from city $$$b_i$$$ to $$$a_i$$$ using the $$$i$$$-th road. Interestingly, all roads satisfy $$$a_i<b_i$$$.Driving EVA requires two people to work together. However, Asuka and Jellyfish have not done any training together before.Suppose that EVA is currently in city $$$u$$$. Jellyfish and Asuka will both choose an undestroyed road that starts at city $$$u$$$. Suppose Jellyfish and Asuka choose roads that end at cities $$$v_1$$$ and $$$v_2$$$ respectively. If $$$v_1 = v_2$$$, EVA moves to city $$$v_1$$$ successfully. Otherwise, EVA stays in city $$$u$$$ and both roads that they have chosen will be destroyed.It is possible that EVA is currently in city $$$u$$$ ($$$u \neq n$$$) and there are no undestroyed roads that start at city $$$u$$$. In that case, the mission will be a failure. Otherwise, if they reach city $$$n$$$ in the end, the mission is considered a success.Every time they choose the roads, Jellyfish knows that Asuka will choose a road randomly. Now, Jellyfish wants to know, if she chooses the roads optimally, what is the maximum probability of the mission being successful.

Input: Each test contains multiple test cases. The first line contains the number of test cases $$$t$$$ ($$$1 \leq t \leq 2000$$$). The description of the test cases follows.The first line of each test case contains two integers, $$$n$$$ and $$$m$$$ ($$$2 \leq n \leq 5000$$$, $$$0 \leq m \leq \min(\frac{n(n-1)}{2}, 2 \cdot 10^5)$$$) — the number of the cities and the number of the roads.In the following $$$m$$$ lines of each test case, each line contains two integers $$$a$$$ and $$$b$$$ ($$$1 \leq a < b \leq n$$$) — representing a road from city $$$a$$$ to city $$$b$$$.It is guaranteed that for each test case, each road appears at most once.It is guaranteed that the sum of $$$n$$$ over all test cases will not exceed $$$5000$$$ and that the sum of $$$m$$$ over all test cases will not exceed $$$2 \cdot 10^5$$$.

Output: Output the maximum probability of the mission being successful if Jellyfish chooses the roads optimally.Your answer will be accepted if the absolute or relative error does not exceed $$$10^{-9}$$$. Formally, let your answer be $$$a$$$, and the jury's answer be $$$b$$$. Your answer is considered correct if $$$\frac{|a-b|}{\max(1,|b|)} \leq 10^{-9}$$$.

Note: In the first test case, Jellyfish will choose $$$v_1=3$$$, and Asuka will choose $$$v_2=2$$$ with a possibility of $$$0.5$$$ and $$$v_2=3$$$ with a possibility of $$$0.5$$$. The mission is considered a success with a possibility of $$$0.5$$$.

Sample Inputs:
3
3 2
1 2
1 3
7 8
1 2
1 3
1 4
1 5
2 6
3 6
4 6
6 7
10 20
1 2
1 3
1 4
1 5
1 6
2 6
2 7
2 8
2 9
3 4
3 7
3 8
3 10
4 6
4 8
4 10
6 10
7 8
7 9
7 10

Sample Outputs:
0.500000000000
0.625000000000
0.491666666667

Solution (Python):
import io
import os
import sys
input = io.BytesIO(os.read(0, os.fstat(0).st_size)).readline
cache = [None] * 5001
cache[0] = [0]
cache[1] = [1]
cache[2] = [0.5, 0]
last_n = 2

def simulate_ps(n):
    if n % 2 == 1:
        return [1 / n] * n
    probs = [0] * n
    p0 = 1 / n
    probs[0] = p0
    previous_ps = cache[n - 2]
    for i in range(1, n):
        p = (i - 1) * previous_ps[i - 2]
        if n - 1 - i > 0:
            p += (n - 1 - i) * previous_ps[i - 1]
        probs[i] = p0 * p
    probs.sort(reverse=True)
    return probs

def estimate_success(town_no):
    global towns
    global connections
    global last_n
    no_of_roads = len(connections[town_no])
    if no_of_roads % 2 == 1:
        ps = [1 / no_of_roads] * no_of_roads
    elif no_of_roads > last_n:
        for i in range(last_n, no_of_roads + 1, 2):
            cache[i] = simulate_ps(i)
        last_n = no_of_roads
        ps = cache[no_of_roads]
    else:
        ps = cache[no_of_roads]
    successes = [0] * no_of_roads
    for i in range(no_of_roads):
        successes[i] = towns[connections[town_no][i]]
    successes.sort(reverse=True)
    success = 0
    for i in range(no_of_roads):
        success += successes[i] * ps[i]
    return success
for _ in range(int(input())):
    (n, m) = map(int, input().split())
    connections = [[] for _ in range(n)]
    for _ in range(m):
        (a, b) = map(int, input().split())
        connections[a] += [b]
    towns = [None] * (n + 1)
    towns[n] = 1
    for i in range(n - 1, 0, -1):
        towns[i] = estimate_success(i)
    sys.stdout.write(f'{towns[1]}\n')

Solution (Python):
nuch = [[0.5, 0]]
a2 = [0.5, 0]
for i in range(2500):
    la = a2.copy()
    a2 = [0] * (i * 2 + 4)
    a2[0] = 1 / len(a2)
    n = len(a2)
    for j in range(1, n):
        if j == 1:
            a2[j] = la[j - 1] * (n - j - 1) / n
        elif j == n - 1:
            a2[j] = la[j - 2] * (j - 1) / n
        else:
            a2[j] = (la[j - 1] * (n - j - 1) + la[j - 2] * (j - 1)) / n
    nuch.append(a2.copy())
for _ in range(int(input())):
    (n, m) = map(int, input().split())
    roads = [[] for i in range(n)]
    for i in range(m):
        (a, b) = map(int, input().split())
        roads[a - 1].append(b - 1)
    ver = [1] * n
    for i in range(n - 2, -1, -1):
        if len(roads[i]) == 0:
            ver[i] = 0
        elif len(roads[i]) % 2 == 0:
            mas = []
            for x in roads[i]:
                mas.append(ver[x])
            mas.sort(reverse=True)
            s = 0
            for j in range(len(mas)):
                s += mas[j] * nuch[len(roads[i]) // 2 - 1][j]
            ver[i] = s
        else:
            s = 0
            for x in roads[i]:
                s += ver[x]
            ver[i] = s / len(roads[i])
    print(ver[0])

Solution (Python):
import sys

def input_general():
    return sys.stdin.readline().rstrip('\r\n')

def input_num():
    return int(sys.stdin.readline().rstrip('\r\n'))

def input_multi(x=int):
    return map(x, sys.stdin.readline().rstrip('\r\n').split())

def input_list(x=int):
    return list(input_multi(x))
coeff = [[], [1.0], [0.5, 0.0], [1 / 3, 1 / 3, 1 / 3]]

def main():
    (n, m) = input_multi()
    edge_dict = [[] for _ in range(n + 1)]
    for _ in range(m):
        (u, v) = input_multi()
        edge_dict[u].append(v)
    prob = [0.0] * n + [1.0]
    for u in range(n - 1, 0, -1):
        nbd = [prob[v] for v in edge_dict[u]]
        if not nbd:
            prob[u] = 0
            continue
        nbd.sort(reverse=True)
        d = len(nbd)
        if len(coeff) <= d:
            for nd in range(len(coeff), d + 1):
                if nd % 2 == 1:
                    coeff.append([1 / nd] * nd)
                else:
                    it = coeff[-2]
                    new = [1 / nd]
                    for (i, (p1, p2)) in enumerate(zip([0] + it, it + [0])):
                        new.append((p1 * i + p2 * (nd - 2 - i)) / nd)
                    coeff.append(new)
        prob[u] = sum((p1 * p2 for (p1, p2) in zip(nbd, coeff[d])))
    print(prob[1])
if __name__ == '__main__':
    cases = input_num()
    for _ in range(cases):
        main()