Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        # Dist = [INF] * (N+1)
        # Dist[start] = 0
        # q = [0<<17 | start]
        # mask = (1<<17)-1
        # while q:
        #     dist_v = heappop(q)
        #     v = dist_v & mask
        #     dist = dist_v >> 17
        #     if Dist[v] != dist:
        #         continue
        #     for u, d in E[v]:
        #         new_dist = dist + d
        #         if Dist[u] > new_dist:
        #             Dist[u] = new_dist
        #             heappush(q, new_dist<<17 | u)
        # return Dist

        import numpy as np
        from scipy.sparse import csr_matrix, csgraph
        V, U, D = np.array(edges).T
        graph = csr_matrix((D, (V, U)), (N+1, N+1))
        return list(csgraph.dijkstra(graph, indices=start, directed=False))

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res += calc_n_paths_S(u)
        res %= mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res += calc_n_paths_T(u)
        res %= mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res += calc_n_paths_S(u)
        res %= mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res += calc_n_paths_T(u)
        res %= mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = sum([calc_n_paths_S(u) for u in DAG_rev[v]]) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = sum([calc_n_paths_T(u) for u in DAG[v]]) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = sum(calc_n_paths_S(u) for u in DAG_rev[v]) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = sum(calc_n_paths_T(u) for u in DAG[v]) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for _ in range(M):
        u, v, d = map(int, input().split())
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
class RadixHeap:
    # from https://ikatakos.com/pot/programming_algorithm/data_structure/radix_heap
    v = [[] for _ in range(70)]

    def __init__(self):
        self.last = 0

    def push(self, x):
        self.v[(x ^ self.last).bit_length()].append(x)

    def pop(self):
        if not self.v[0]:
            cv = next(x for x in self.v if x)
            self.last = min(cv)
            for x in cv:
                self.push(x)
            cv.clear()
        return self.v[0].pop()


def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for _ in range(M):
        u, v, d = map(int, input().split())
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        #q = [0<<20 | start]
        q = RadixHeap()
        q.push(1<<68)  # 番兵
        q.push(0<<20 | start)
        mask = (1<<20)-1
        while True:
            #dist_v = heappop(q)
            dist_v = q.pop()
            if dist_v == 1<<68:
                break
            v = dist_v & mask
            dist = dist_v >> 20
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    #heappush(q, new_dist<<20 | u)
                    q.push(new_dist<<20 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)

    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res

    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)

    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for _ in range(M):
        u, v, d = map(int, input().split())
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))
    
    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<20 | start]
        mask = (1<<20)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 20
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<20 | u)
        return Dist
    
    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]
    
    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)
    
    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res
    
    ans = calc_n_paths_S(T)  # 全経路数
    calc_n_paths_T(S)
    
    ans = ans * ans % mod
    
    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod
    
    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod
    
    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
# heapq の中身を 1 つに
import sys
sys.setrecursionlimit(500000)
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [0<<20 | start]
    mask = (1<<20)-1
    while q:
        dist_v = heappop(q)
        v = dist_v & mask
        dist = dist_v >> 20
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, new_dist<<20 | u)
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v))
        DAG[u].append(v)
        DAG_rev[v].append(u)
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u))
        DAG[v].append(u)
        DAG_rev[u].append(v)

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
calc_n_paths_T(S)

ans = ans * ans % mod

for v, u in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist*2 == dist_st:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
sys.setrecursionlimit(500000)
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v))
        DAG[u].append(v)
        DAG_rev[v].append(u)
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u))
        DAG[v].append(u)
        DAG_rev[u].append(v)

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
calc_n_paths_T(S)

ans = ans * ans % mod

for v, u in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist*2 == dist_st:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
sys.setrecursionlimit(500000)
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v))
        DAG[u].append(v)
        DAG_rev[v].append(u)
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u))
        DAG[v].append(u)
        DAG_rev[u].append(v)

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

for v, u in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist*2 == dist_st:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
def main():
    import sys
    sys.setrecursionlimit(500000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    for _ in range(M):
        u, v, d = map(int, input().split())
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))
    
    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [(0, start)]
        while q:
            dist, v = heappop(q)
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, (new_dist, u))
        return Dist
    
    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]
    
    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v, d))
            DAG[u].append((v, d))
            DAG_rev[v].append((u, d))
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u, d))
            DAG[v].append((u, d))
            DAG_rev[u].append((v, d))
    
    n_paths_S = [-1] * (N+1)
    n_paths_T = [-1] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    def calc_n_paths_S(v):
        if n_paths_S[v] != -1:
            return n_paths_S[v]
        res = 0
        for u, _ in DAG_rev[v]:
            res = (res + calc_n_paths_S(u)) % mod
        n_paths_S[v] = res
        return res
    def calc_n_paths_T(v):
        if n_paths_T[v] != -1:
            return n_paths_T[v]
        res = 0
        for u, _ in DAG[v]:
            res = (res + calc_n_paths_T(u)) % mod
        n_paths_T[v] = res
        return res
    
    ans = calc_n_paths_S(T)  # 全経路数
    assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)
    
    ans = ans * ans % mod
    
    for v, u, _ in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod
    
    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod
    
    print(ans)

main()


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
sys.setrecursionlimit(500000)
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v, d))
        DAG[u].append((v, d))
        DAG_rev[v].append((u, d))
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u, d))
        DAG[v].append((u, d))
        DAG_rev[u].append((v, d))

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u, _ in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u, _ in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

for v, u, _ in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist*2 == dist_st:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Runtime Error):
import sys
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v, d))
        DAG[u].append((v, d))
        DAG_rev[v].append((u, d))
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u, d))
        DAG[v].append((u, d))
        DAG_rev[u].append((v, d))

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u, _ in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u, _ in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
calc_n_paths_T(S)
#assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

for v, u, _ in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist*2 == dist_st:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Runtime Error):
import sys
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v, d))
        DAG[u].append((v, d))
        DAG_rev[v].append((u, d))
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u, d))
        DAG[v].append((u, d))
        DAG_rev[u].append((v, d))

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u, _ in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u, _ in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

for v, u, _ in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist*2 == dist_st:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Runtime Error):
import sys
input = sys.stdin.readline
from heapq import heappop, heappush
mod = 10**9+7
INF = 1<<60
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [INF] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v, d))
        DAG[u].append((v, d))
        DAG_rev[v].append((u, d))
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u, d))
        DAG[v].append((u, d))
        DAG_rev[u].append((v, d))

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u, _ in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u, _ in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

half = dist_st / 2.0
assert half * 2 == dist_st
for v, u, _ in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v] < half and half < Dist_S[u]:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist == half:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Runtime Error):
from heapq import heappop, heappush
mod = 10**9+7
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [float("inf")] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v, d))
        DAG[u].append((v, d))
        DAG_rev[v].append((u, d))
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u, d))
        DAG[v].append((u, d))
        DAG_rev[u].append((v, d))

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u, _ in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u, _ in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

half = dist_st / 2.0
assert half * 2 == dist_st
for v, u, _ in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v] < half and half < Dist_S[u]:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

# 頂点ですれ違う場合
for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    if dist == half:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Runtime Error):
from heapq import heappop, heappush
mod = 10**9+7
N, M = map(int, input().split())
S, T = map(int, input().split())
E = [[] for _ in range(N+1)]
edges = []
for _ in range(M):
    u, v, d = map(int, input().split())
    E[u].append((v, d))
    E[v].append((u, d))
    edges.append((u, v, d))

def dijksrtra(start):
    Dist = [float("inf")] * (N+1)
    Dist[start] = 0
    q = [(0, start)]
    while q:
        dist, v = heappop(q)
        if Dist[v] != dist:
            continue
        for u, d in E[v]:
            new_dist = dist + d
            if Dist[u] > new_dist:
                Dist[u] = new_dist
                heappush(q, (new_dist, u))
    return Dist

Dist_S = dijksrtra(S)
Dist_T = dijksrtra(T)
dist_st = Dist_S[T]

DAG_edges = []  # S -> T
DAG = [[] for _ in range(N+1)]
DAG_rev = [[] for _ in range(N+1)]
for u, v, d in edges:
    if Dist_S[u] + Dist_T[v] + d == dist_st:
        DAG_edges.append((u, v, d))
        DAG[u].append((v, d))
        DAG_rev[v].append((u, d))
    elif Dist_T[u] + Dist_S[v] + d == dist_st:
        DAG_edges.append((v, u, d))
        DAG[v].append((u, d))
        DAG_rev[u].append((v, d))

n_paths_S = [-1] * (N+1)
n_paths_T = [-1] * (N+1)
n_paths_S[S] = 1
n_paths_T[T] = 1
def calc_n_paths_S(v):
    if n_paths_S[v] != -1:
        return n_paths_S[v]
    res = 0
    for u, _ in DAG_rev[v]:
        res = (res + calc_n_paths_S(u)) % mod
    n_paths_S[v] = res
    return res
def calc_n_paths_T(v):
    if n_paths_T[v] != -1:
        return n_paths_T[v]
    res = 0
    for u, _ in DAG[v]:
        res = (res + calc_n_paths_T(u)) % mod
    n_paths_T[v] = res
    return res

ans = calc_n_paths_S(T)  # 全経路数
assert ans == calc_n_paths_T(S), (n_paths_S, n_paths_T)

ans = ans * ans % mod

half = dist_st / 2.0
assert half * 2 == dist_st
for v, u, _ in DAG_edges:  # 辺ですれ違う場合
    if Dist_S[v] < half and half < Dist_S[u]:
        ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
    # 頂点ですれ違う場合
    if dist == half:
        ans = (ans - (ns*nt)**2) % mod

print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
import heapq
N, M = map(int, input().split())
S, T = map(int, input().split())
E_ = [list(map(int, input().split())) for _ in range(M)]
E = [[] for _ in range(N+1)]
for U, V, D in E_:
    E[U].append([V, D])
    E[V].append([U, D])

E_S = [[] for _ in range(N+1)]
dist_S = [float("inf") for _ in range(N+1)]
dist_S[S] = 0
cnt_S = [0] * (N+1)
cnt_S[S] = 1
q = [(0, S, 0, 0)]
dist_st = float("inf")
closed = [False] * (N+1)
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_S[v] < prov_cost:
        continue

    cnt_S[v] += cnt_S[from_]
    E_S[v].append([from_, d])

    if closed[v]:
        continue
    if dist_st < prov_cost:
        break
    if v == T:
        dist_st = prov_cost
    closed[v] = True
    for u, d in E[v]:
        cost_new = prov_cost + d
        if dist_S[u] >= cost_new:
            dist_S[u] = cost_new
            heapq.heappush(q, (cost_new, u, v, d))
cnt_st = cnt_S[T]

d2 = (dist_st+1)//2
dist_T = [float("inf") for _ in range(N+1)]
cnt_T = [0] * (N+1)
cnt_T[T] = 1
q = [(0, T, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_T[v] < prov_cost:
        continue

    cnt_T[v] += cnt_T[from_]

    if dist_T[v] == prov_cost:
        continue
    dist_T[v] = prov_cost
    if d2 < prov_cost:
        break
    for u, d in E[v]:
        if dist_T[u] > prov_cost + d:
            heapq.heappush(q, (prov_cost+d, u, v, d))

ans = cnt_st * cnt_st

if dist_st % 2 == 0:
    d2 = dist_st//2
    for i, (ds, dt) in enumerate(zip(dist_S, dist_T)):
        if d2==ds==dt:
            c = cnt_S[i] * cnt_T[i]
            ans -= c*c
d2 = (dist_st+1) // 2
for u, v, d in E_:
    dsu = dist_S[u]
    dtv = dist_T[v]
    if dsu < d2 and dtv < d2 and dsu + dtv + d == dist_st:
        c = cnt_S[u] * cnt_T[v]
        ans -= c*c
    dsv = dist_S[v]
    dtu = dist_T[u]
    if dsv < d2 and dtu < d2 and dsv +dtu + d == dist_st:
        c = cnt_S[v] * cnt_T[u]
        ans -= c*c
print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
import heapq
N, M = map(int, input().split())
S, T = map(int, input().split())
E_ = [list(map(int, input().split())) for _ in range(M)]
E = [[] for _ in range(N+1)]
for U, V, D in E_:
    E[U].append([V, D])
    E[V].append([U, D])

E_S = [[] for _ in range(N+1)]
dist_S = [float("inf") for _ in range(N+1)]
cnt_S = [0] * (N+1)
cnt_S[S] = 1
q = [(0, S, 0, 0)]
dist_st = float("inf")
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_S[v] < prov_cost:
        continue

    cnt_S[v] += cnt_S[from_]
    E_S[v].append([from_, d])

    if dist_S[v] == prov_cost:
        continue
    if dist_st < prov_cost:
        break
    dist_S[v] = prov_cost
    if v == T:
        dist_st = prov_cost
    for u, d in E[v]:
        if dist_S[u] > prov_cost + d:
            heapq.heappush(q, (prov_cost+d, u, v, d))

cnt_st = cnt_S[T]

d2 = (dist_st+1)//2
dist_T = [float("inf") for _ in range(N+1)]
cnt_T = [0] * (N+1)
cnt_T[T] = 1
q = [(0, T, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_T[v] < prov_cost:
        continue

    cnt_T[v] += cnt_T[from_]

    if dist_T[v] == prov_cost:
        continue
    dist_T[v] = prov_cost
    if d2 < prov_cost:
        break
    for u, d in E[v]:
        if dist_T[u] > prov_cost + d:
            heapq.heappush(q, (prov_cost+d, u, v, d))

ans = cnt_st * cnt_st

if dist_st % 2 == 0:
    d2 = dist_st//2
    for i, (ds, dt) in enumerate(zip(dist_S, dist_T)):
        if d2==ds==dt:
            c = cnt_S[i] * cnt_T[i]
            ans -= c*c
d2 = (dist_st+1) // 2
for u, v, d in E_:
    dsu = dist_S[u]
    dtv = dist_T[v]
    if dsu < d2 and dtv < d2 and dsu + dtv + d == dist_st:
        c = cnt_S[u] * cnt_T[v]
        ans -= c*c
    dsv = dist_S[v]
    dtu = dist_T[u]
    if dsv < d2 and dtu < d2 and dsv +dtu + d == dist_st:
        c = cnt_S[v] * cnt_T[u]
        ans -= c*c
print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
import heapq
N, M = map(int, input().split())
S, T = map(int, input().split())
E_ = [list(map(int, input().split())) for _ in range(M)]
E = [[] for _ in range(N+1)]
for U, V, D in E_:
    E[U].append([V, D])
    E[V].append([U, D])

E_S = [[] for _ in range(N+1)]
dist_S = [float("inf") for _ in range(N+1)]
cnt_S = [0] * (N+1)
cnt_S[S] = 1
q = [(0, S, 0, 0)]
dist_st = float("inf")
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_S[v] < prov_cost:
        continue

    cnt_S[v] += cnt_S[from_]
    E_S[v].append([from_, d])

    if dist_S[v] == prov_cost:
        continue
    if dist_st < prov_cost:
        break
    dist_S[v] = prov_cost
    if v == T:
        dist_st = prov_cost
    for u, d in E[v]:
        if dist_S[u] == float("inf"):
            heapq.heappush(q, (prov_cost+d, u, v, d))

cnt_st = cnt_S[T]

d2 = (dist_st+1)//2
dist_T = [float("inf") for _ in range(N+1)]
cnt_T = [0] * (N+1)
cnt_T[T] = 1
q = [(0, T, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_T[v] < prov_cost:
        continue

    cnt_T[v] += cnt_T[from_]

    if dist_T[v] == prov_cost:
        continue
    dist_T[v] = prov_cost
    if d2 < prov_cost:
        break
    for u, d in E[v]:
        if dist_T[u] == float("inf"):
            heapq.heappush(q, (prov_cost+d, u, v, d))

ans = cnt_st * cnt_st

if dist_st % 2 == 0:
    d2 = dist_st//2
    for i, (ds, dt) in enumerate(zip(dist_S, dist_T)):
        if d2==ds==dt:
            c = cnt_S[i] * cnt_T[i]
            ans -= c*c
d2 = (dist_st+1) // 2
for u, v, d in E_:
    dsu = dist_S[u]
    dtv = dist_T[v]
    if dsu < d2 and dtv < d2 and dsu + dtv + d == dist_st:
        c = cnt_S[u] * cnt_T[v]
        ans -= c*c
    dsv = dist_S[v]
    dtu = dist_T[u]
    if dsv < d2 and dtu < d2 and dsv +dtu + d == dist_st:
        c = cnt_S[v] * cnt_T[u]
        ans -= c*c
print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import heapq
import sys
input = sys.stdin.readline
N, M = map(int, input().split())
S, T = map(int, input().split())
E_ = [list(map(int, input().split())) for _ in range(M)]
E = [[] for _ in range(N+1)]
for U, V, D in E_:
    E[U].append([V, D])
    E[V].append([U, D])

E_S = [[] for _ in range(N+1)]
dist_S = [float("inf") for _ in range(N+1)]
cnt_S = [0] * (N+1)
cnt_S[S] = 1
q = [(0, S, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_S[v] < prov_cost:
        continue

    cnt_S[v] += cnt_S[from_]
    E_S[v].append([from_, d])

    if dist_S[v] == prov_cost:
        continue
    dist_S[v] = prov_cost

    for u, d in E[v]:
        if dist_S[u] == float("inf"):
            heapq.heappush(q, (prov_cost+d, u, v, d))

dist_st = dist_S[T]
cnt_st = cnt_S[T]

dist_T = [float("inf") for _ in range(N+1)]
cnt_T = [0] * (N+1)
cnt_T[T] = 1
q = [(0, T, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_T[v] < prov_cost:
        continue

    cnt_T[v] += cnt_T[from_]

    if dist_T[v] == prov_cost:
        continue
    dist_T[v] = prov_cost

    for u, d in E[v]:
        if dist_T[u] == float("inf"):
            heapq.heappush(q, (prov_cost+d, u, v, d))

ans = cnt_st * cnt_st

if dist_st % 2 == 0:
    d2 = dist_st//2
    for i, (ds, dt) in enumerate(zip(dist_S, dist_T)):
        if d2==ds==dt:
            c = cnt_S[i] * cnt_T[i]
            ans -= c*c
d2 = (dist_st+1) // 2
for u, v, d in E_:
    if dist_S[u] < d2 and dist_T[v] < d2 and dist_S[u] + dist_T[v] + d == dist_st:
        c = cnt_S[u] * cnt_T[v]
        ans -= c*c
    if dist_S[v] < d2 and dist_T[u] < d2 and dist_S[v] + dist_T[u] + d == dist_st:
        c = cnt_S[v] * cnt_T[u]
        ans -= c*c
print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Buggy python submission (verdict: Time Limit Exceeded):
import heapq
N, M = map(int, input().split())
S, T = map(int, input().split())
E_ = [list(map(int, input().split())) for _ in range(M)]
E = [[] for _ in range(N+1)]
for U, V, D in E_:
    E[U].append([V, D])
    E[V].append([U, D])

E_S = [[] for _ in range(N+1)]
dist_S = [float("inf") for _ in range(N+1)]
cnt_S = [0] * (N+1)
cnt_S[S] = 1
q = [(0, S, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_S[v] < prov_cost:
        continue

    cnt_S[v] += cnt_S[from_]
    E_S[v].append([from_, d])

    if dist_S[v] == prov_cost:
        continue
    dist_S[v] = prov_cost

    for u, d in E[v]:
        if dist_S[u] == float("inf"):
            heapq.heappush(q, (prov_cost+d, u, v, d))

dist_st = dist_S[T]
cnt_st = cnt_S[T]

dist_T = [float("inf") for _ in range(N+1)]
cnt_T = [0] * (N+1)
cnt_T[T] = 1
q = [(0, T, 0, 0)]
while len(q)!=0:
    prov_cost, v, from_, d = heapq.heappop(q)
    if dist_T[v] < prov_cost:
        continue

    cnt_T[v] += cnt_T[from_]

    if dist_T[v] == prov_cost:
        continue
    dist_T[v] = prov_cost

    for u, d in E[v]:
        if dist_T[u] == float("inf"):
            heapq.heappush(q, (prov_cost+d, u, v, d))

ans = cnt_st * cnt_st

if dist_st % 2 == 0:
    d2 = dist_st//2
    for i, (ds, dt) in enumerate(zip(dist_S, dist_T)):
        if d2==ds==dt:
            c = cnt_S[i] * cnt_T[i]
            ans -= c*c
d2 = (dist_st+1) // 2
for u, v, d in E_:
    if dist_S[u] < d2 and dist_T[v] < d2 and dist_S[u] + dist_T[v] + d == dist_st:
        c = cnt_S[u] * cnt_T[v]
        ans -= c*c
    if dist_S[v] < d2 and dist_T[u] < d2 and dist_S[v] + dist_T[u] + d == dist_st:
        c = cnt_S[v] * cnt_T[u]
        ans -= c*c
print(ans)


Fixed submission:
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    # for _ in range(M):
    #     u, v, d = map(int, input().split())
    #     E[u].append((v, d))
    #     E[v].append((u, d))
    #     edges.append((u, v, d))
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  # S -> T
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    # トポロジカルソート
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    # n_paths_S = [-1] * (N+1)
    # n_paths_T = [-1] * (N+1)
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    # def calc_n_paths_S(v):
    #     if n_paths_S[v] != -1:
    #         return n_paths_S[v]
    #     res = 0
    #     for u in DAG_rev[v]:
    #         res = (res + calc_n_paths_S(u)) % mod
    #     n_paths_S[v] = res
    #     return res
    # def calc_n_paths_T(v):
    #     if n_paths_T[v] != -1:
    #         return n_paths_T[v]
    #     res = 0
    #     for u in DAG[v]:
    #         res = (res + calc_n_paths_T(u)) % mod
    #     n_paths_T[v] = res
    #     return res

    # ans = calc_n_paths_S(T)  # 全経路数
    # calc_n_paths_T(S)

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  # 辺ですれ違う場合
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    # 頂点ですれ違う場合
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()