Jellyfish and EVA (Codeforces 1875F)
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 sys
input = sys.stdin.readline
import math
t = int(input())
for _ in range(t):
    [n, m] = [int(x) for x in input().split()]
    neighbors = [[] for i in range(n + 1)]
    for j in range(m):
        [a, b] = [int(x) for x in input().split()]
        neighbors[a].append(b)
    probabilities = [0 for i in range(n + 1)]
    probabilities[n] = 1
    coeffs = [[], [1], [0.5, 0]]
    for i in range(3, n):
        new_coeffs = [1 / i]
        for j in range(1, i):
            new_coeffs.append((j - 1) / i * coeffs[i - 2][j - 2] + (i - j - 1) / i * coeffs[i - 2][min(j - 1, i - 3)])
        coeffs.append(new_coeffs)
    for j in range(n - 1, 0, -1):
        vals = [probabilities[neighbors[j][i]] for i in range(len(neighbors[j]))]
        if len(vals) != 0:
            vals.sort()
            vals.reverse()
            val = 0
            for i in range(len(vals)):
                val += vals[i] * coeffs[len(vals)][i]
            probabilities[j] = val
    print(probabilities[1])

Solution (Python):
import sys
ipt = sys.stdin.readline

def iin():
    return int(ipt())

def lmin():
    return list(map(int, ipt().split()))
sz = 5001
pre = [0] * (sz * sz + 10)
pre[sz * 1 + 1] = 1
pre[sz * 2 + 1] = 0.5
for i in range(3, sz):
    pre[i * sz + 1] = 1 / i
    for j in range(2, i + 1):
        pre[i * sz + j] = pre[(i - 2) * sz + j - 2] * (j - 2) / i + pre[(i - 2) * sz + j - 1] * (i - j) / i

def solve():
    (N, M) = lmin()
    G = [[] for _ in range(N)]
    RG = [[] for _ in range(N)]
    deg = [0] * N
    for _ in range(M):
        (a, b) = lmin()
        (a, b) = (a - 1, b - 1)
        G[a].append(b)
        RG[b].append(a)
    st = [N - 1]
    while st:
        cur = st.pop()
        for nxt in RG[cur]:
            deg[nxt] += 1
            if deg[nxt] == 1:
                st.append(nxt)
    st = [N - 1]
    dp = [0] * N
    dp[N - 1] = 1
    st = [N - 1]
    while st:
        cur = st.pop()
        if cur != N - 1:
            l = len(G[cur])
            candidate = []
            for nxt in G[cur]:
                candidate.append(dp[nxt])
            candidate.sort(reverse=True)
            dp[cur] = sum((candidate[i] * pre[l * sz + i + 1] for i in range(l)))
        for nxt in RG[cur]:
            deg[nxt] -= 1
            if deg[nxt] == 0:
                st.append(nxt)
    print(dp[0])
for _ in range(iin()):
    solve()

Solution (Python):
from __future__ import print_function
from math import *
from collections import defaultdict, deque
import os
import random
import sys
from io import BytesIO, IOBase
from types import GeneratorType

def bootstrap(f, stack=[]):

    def wrappedfunc(*args, **kwargs):
        if stack:
            return f(*args, **kwargs)
        else:
            to = f(*args, **kwargs)
            while True:
                if type(to) is GeneratorType:
                    stack.append(to)
                    to = next(to)
                else:
                    stack.pop()
                    if not stack:
                        break
                    to = stack[-1].send(to)
            return to
    return wrappedfunc

def main():
    pass
BUFSIZE = 8192

def lcm(a, b):
    return a * b // gcd(a, b)

def ceilDiv(a, b):
    if a % b == 0:
        return a // b
    else:
        return a // b + 1

class FastIO(IOBase):
    newlines = 0

    def __init__(self, file):
        self._fd = file.fileno()
        self.buffer = BytesIO()
        self.writable = 'x' in file.mode or 'r' not in file.mode
        self.write = self.buffer.write if self.writable else None

    def read(self):
        while True:
            b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
            if not b:
                break
            ptr = self.buffer.tell()
            (self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
        self.newlines = 0
        return self.buffer.read()

    def readline(self):
        while self.newlines == 0:
            b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
            self.newlines = b.count(b'\n') + (not b)
            ptr = self.buffer.tell()
            (self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
        self.newlines -= 1
        return self.buffer.readline()

    def flush(self):
        if self.writable:
            os.write(self._fd, self.buffer.getvalue())
            (self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

    def __init__(self, file):
        self.buffer = FastIO(file)
        self.flush = self.buffer.flush
        self.writable = self.buffer.writable
        self.write = lambda s: self.buffer.write(s.encode('ascii'))
        self.read = lambda : self.buffer.read().decode('ascii')
        self.readline = lambda : self.buffer.readline().decode('ascii')

class myDict:

    def __init__(self, func):
        self.RANDOM = random.randint(0, 1 << 32)
        self.default = func
        self.dict = {}

    def items(self):
        return [(i ^ self.RANDOM, j) for (i, j) in self.dict.items()]

    def __getitem__(self, key):
        myKey = self.RANDOM ^ key
        if myKey not in self.dict:
            self.dict[myKey] = self.default()
        return self.dict[myKey]

    def __setitem__(self, key, item):
        myKey = self.RANDOM ^ key
        self.dict[myKey] = item

    def __iter__(self):
        return iter((i ^ self.RANDOM for i in self.dict))

    def __contains__(self, key):
        return key ^ self.RANDOM in self.dict

    def __repr__(self):
        return repr({i ^ self.RANDOM: j for (i, j) in self.dict.items()})

class mySet:

    def __init__(self):
        self.random = random.randint(0, 1 << 32)
        self.set = set()

    def add(self, key):
        self.set.add(key ^ self.random)

    def remove(self, key):
        self.set.remove(key ^ self.random)

    def __contains__(self, key):
        return key ^ self.random in self.set

    def __iter__(self):
        return iter((i ^ self.random for i in self.set))

    def __len__(self):
        return len(self.set)

    def __repr__(self):
        s = set()
        for item in self.set:
            s.add(item ^ self.random)
        return repr(s)
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')
mod = 10 ** 9 + 7

def inp():
    return int(input())

def inlt():
    return list(map(int, input().split()))

def insr():
    s = input()
    return list(s[:len(s)])

def invr():
    return map(int, input().split())
N = 5001
dpp = [[0] * (i + 2) for i in range(N)]
dpp[1][1] = 1
dpp[2][1] = 0.5
dpp[2][2] = 0
for i in range(3, N):
    for j in range(1, i + 1):
        if j == 1:
            dpp[i][j] = 1 / i
        else:
            dpp[i][j] = dpp[i - 2][j - 2] * (j - 2) / i + dpp[i - 2][j - 1] * (i - j) / i
for _ in range(inp()):
    (n, m) = invr()
    l = [[] for i in range(n)]
    for i in range(m):
        (a, b) = invr()
        l[a - 1].append(b - 1)
    dp = [0] * n
    dp[n - 1] = 1
    for i in range(n - 2, -1, -1):
        t = []
        for j in l[i]:
            t.append(dp[j])
        t.sort(reverse=1)
        c = 0
        for j in t:
            c += 1
            dp[i] += j * dpp[len(t)][c]
    print(dp[0])