Problem:
Problem G: Strange Couple Alice and Bob are going to drive from their home to a theater for a date. They are very challenging - they have no maps with them even though they don’t know the route at all (since they have just moved to
their new home). Yes, they will be going just by their feeling. The town they drive can be considered as an undirected graph with a number of intersections (vertices)
and roads (edges). Each intersection may or may not have a sign. On intersections with signs, Alice and
Bob will enter the road for the shortest route. When there is more than one such roads, they will go into
one of them at random. On intersections without signs, they will just make a random choice. Each random selection is made with
equal probabilities. They can even choose to go back to the road they have just come along, on a random
selection. Calculate the expected distance Alice and Bob will drive before reaching the theater. Input The input consists of multiple datasets. Each dataset has the following format:  n s t  q 1  q 2 ... q n  a 11  a 12 ... a 1 n  a 21  a 22 ... a 2 n . . .  a n 1  a n 2 ... a nn   n is the number of intersections ( n ≤ 100). s and t are the intersections the home and the theater are
located respectively (1 ≤ s , t ≤ n , s ≠ t ); q i (for 1 ≤ i ≤ n ) is either 1 or 0, where 1 denotes there is a
sign at the i -th intersection and 0 denotes there is not; a ij (for 1 ≤ i , j ≤ n ) is a positive integer denoting
the distance of the road connecting the i -th and j -th intersections, or 0 indicating there is no road directly
connecting the intersections. The distance of each road does not exceed 10. Since the graph is undirectional, it holds a ij = a ji for any 1 ≤ i , j ≤ n . There can be roads connecting the
same intersection, that is, it does not always hold a ii = 0. Also, note that the graph is not always planar. The last dataset is followed by a line containing three zeros. This line is not a part of any dataset and should not be processed. Output For each dataset, print the expected distance accurate to 10 -8 , or " impossible " (without quotes) if there is no route to get
to the theater. The distance may be printed with any number of digits after the decimal point. Sample Input 5 1 5
1 0 1 0 0
0 2 1 0 0
2 0 0 1 0
1 0 0 1 0
0 1 1 0 1
0 0 0 1 0
0 0 0 Output for the Sample Input 8.50000000

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,copy,functools
import time,random

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
mod2 = 998244353
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(): return [list(map(int, l.split())) for l in sys.stdin.readlines()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)
def pe(s): return print(str(s), file=sys.stderr)
def JA(a, sep): return sep.join(map(str, a))
def JAA(a, s, t): return s.join(t.join(map(str, b)) for b in a)




def gauss_jordan(A, b):
    n = len(A)
    B = [A[i][:] for i in range(n)]
    for i in range(n):
        B[i].append(b[i])

    for i in range(n):
        pivot = i
        abp = abs(B[i][i])
        for j in range(i+1, n):
            if abp < abs(B[j][i]):
                abp = abs(B[j][i])
                pivot = j
        B[i],B[pivot] = B[pivot],B[i]
        if abp < eps:
            return

        for j in range(i+1, n+1):
            B[i][j] /= B[i][i]
        for j in range(n):
            if j == i:
                continue
            for k in range(i+1, n+1):
                B[j][k] -= B[j][i] * B[i][k]

    return list(map(lambda x: x[-1], B))

def main():
    def f():
        n,s,t = LI()
        if n < 1:
            return None

        s -= 1
        t -= 1
        p = LI()
        aa = [LI() for _ in range(n)]

        e = collections.defaultdict(list)
        for i in range(n):
            ai = aa[i]
            for j in range(n):
                if ai[j] == 0:
                    continue
                e[i].append((j, ai[j]))

        def search(s):
            d = collections.defaultdict(lambda: inf)
            d[s] = 0
            q = []
            heapq.heappush(q, (0, s))
            v = collections.defaultdict(bool)
            while len(q):
                k, u = heapq.heappop(q)
                if v[u]:
                    continue
                v[u] = True

                for uv, ud in e[u]:
                    if v[uv]:
                        continue
                    vd = k + ud
                    if d[uv] > vd:
                        d[uv] = vd
                        heapq.heappush(q, (vd, uv))

            return d

        d = search(t)
        if d[s] == inf:
            return "impossible"

        A = [[0]*n for _ in range(n)]
        b = [0] * n
        for i in range(n):
            if i == t or d[i] == inf:
                A[i][i] = 1
                b[i] = 0
                continue

            mv = 0
            mk = 0
            di = d[i]
            ai = aa[i]
            for j in range(n):
                if ai[j] < 1 or (p[i] == 1 and d[j] + ai[j] != di):
                    continue
                mv += 1
                mk += ai[j]
                A[i][j] = -1

            A[i][i] = mv
            b[i] = mk

        r = gauss_jordan(A, b)
        if r is None:
            return "impossible"

        return "{:0.9f}".format(r[s])

    rr = []
    while True:
        r = f()
        if r is None:
            break
        rr.append(r)

    return JA(rr, "\n")


print(main())

Solution (Python):
from heapq import heappush, heappop
import sys
readline = sys.stdin.readline
write = sys.stdout.write

def solve():
    N, s, t = map(int, readline().split())
    if N == 0:
        return False
    *Q, = map(int, readline().split())
    G = [[] for i in range(N)]
    for v in range(N):
        *A, = map(int, readline().split())
        for w, a in enumerate(A):
            if not a:
                continue
            G[v].append((w, a))
    INF = 10**9
    que = [(0, t-1)]
    dist = [INF]*N
    dist[t-1] = 0

    while que:
        cost, v = heappop(que)
        if dist[v] < cost:
            continue
        for w, d in G[v]:
            if cost + d < dist[w]:
                dist[w] = cost + d
                heappush(que, (cost + d, w))

    if dist[s-1] == INF:
        write("impossible\n")
        return

    G0 = [[] for i in range(N)]
    for v in range(N):
        if v == t-1:
            continue
        if Q[v]:
            cost = dist[v]
            for w, d in G[v]:
                if cost != dist[w] + d:
                    continue
                G0[v].append((w, d))
        else:
            G0[v] = G[v]

    MT = [[0]*(N+1) for i in range(N)]
    for v in range(N):
        Mv = MT[v]
        Mv[v] = 1
        r = 0
        g = len(G0[v])
        for w, d in G0[v]:
            Mv[w] -= 1 / g
            r += d / g
        Mv[N] = r

    for i in range(N):
        Mi = MT[i]
        v = Mi[i]
        for j in range(N+1):
            Mi[j] /= v
        for k in range(N):
            Mk = MT[k]
            e = Mk[i]
            if e == 0 or i == k:
                continue
            for j in range(N+1):
                Mk[j] -= e * Mi[j]
    write("%.16f\n" % MT[s-1][N])
    return True
while solve():
    ...