def write_graph(dotfile, graph, directed=False,
                node_label=None, arc_label=None, comment="",
                node_mark=set(), arc_mark=set()):
    """Writes a graph to a file in the DOT format

    :param dotfile: the filename.
    :param graph: directed graph in listlist or listdict format
    :param directed: true if graph is directed, false if undirected
    :param weight: in matrix format or same listdict graph or None
    :param node_label: vertex label table or None
    :param arc_label: arc label matrix or None
    :param comment: comment string for the dot file or None
    :param node_mark: set of nodes to be shown in gray
    :param arc_marc: set of arcs to be shown in red
    :complexity: `O(|V| + |E|)`
    """
    with open(dotfile, 'w') as f:
        if directed:
            f.write("digraph G{\n")
        else:
            f.write("graph G{\n")
        if comment:
            f.write('label="%s";\n' % comment)
        V = range(len(graph))
        #                              -- vertices
        for u in V:
            if node_mark and u in node_mark:
                f.write('%d [style=filled, color="lightgrey", ' % u)
            else:
                f.write('%d [' % u)
            if node_label:
                f.write('label="%u [%s]"];\n' % (u, node_label[u]))
            else:
                f.write('shape=circle, label="%u"];\n' % u)
        #                              -- edges
        if isinstance(arc_mark, list):
            arc_mark = set((u, arc_mark[u]) for u in V)
        for u in V:
            for v in graph[u]:
                if not directed and u > v:
                    continue   # don't show twice the edge
                if arc_label and arc_label[u][v] is None:
                    continue   # suppress arcs with no label
                if directed:
                    arc = "%d -> %d " % (u, v)
                else:
                    arc = "%d -- %d " % (u, v)
                if arc_mark and ((v, u) in arc_mark or
                                 (not directed and (u, v) in arc_mark)):
                    pen = 'color="red"'
                else:
                    pen = ""
                if arc_label:
                    tag = 'label="%s"' % arc_label[u][v]
                else:
                    tag = ""
                if tag and pen:
                    sep = ", "
                else:
                    sep = ""
                f.write(arc + "[" + tag + sep + pen + "];\n")
        f.write("}") def extract_path(prec, v):
    """extracts a path in form of vertex list from source to vertex v
       given a precedence table prec leading to the source

    :param prec: precedence table of a tree
    :param v: vertex on the tree
    :returns: path from root to v, in form of a list
    :complexity: linear
    """
    L = []
    while v is not None:
        L.append(v)
        v = prec[v]
        assert v not in L  # prevent infinite loops for a bad formed table prec
    return L[::-1] def make_flow_labels(graph, flow, capac):
    """Generate arc labels for a flow in a graph with capacities.

    :param graph: adjacency list or adjacency dictionary
    :param flow:  flow matrix or adjacency dictionary
    :param capac: capacity matrix or adjacency dictionary
    :returns: listdic graph representation, with the arc label strings
    """
    V = range(len(graph))
    arc_label = [{v: "" for v in graph[u]} for u in V]
    for u in V:
        for v in graph[u]:
            if flow[u][v] >= 0:
                arc_label[u][v] = "%s/%s" % (flow[u][v], capac[u][v])
            else:
                arc_label[u][v] = None   # do not show negative flow arcs
    return arc_label def tree_adj_to_prec(graph, root=0):
    """Transforms a tree given as adjacency list into predecessor table form.
    if graph is not a tree: will return a DFS spanning tree

    :param graph: directed graph in listlist or listdict format
    :returns: tree in predecessor table representation
    :complexity: linear
    """
    prec = [None] * len(graph)
    prec[root] = root            # mark to visit root only once
    to_visit = [root]
    while to_visit:              # DFS
        node = to_visit.pop()
        for neighbor in graph[node]:
            if prec[neighbor] is None:
                prec[neighbor] = node
                to_visit.append(neighbor)
    prec[root] = None            # put the standard mark for root
    return prec def tree_prec_to_adj(prec, root=0):
    """Transforms a tree given as predecessor table into adjacency list form

    :param prec: predecessor table representing a tree, prec[u] == v iff u is
    successor of v,
                 except for the root where prec[root] == root
    :param root: root vertex of the tree
    :returns: undirected graph in listlist representation
    :complexity: linear
    """
    n = len(prec)
    graph = [[prec[u]] for u in range(n)]   # add predecessors
    graph[root] = []
    for u in range(n):                      # add successors
        if u != root:
            graph[prec[u]].append(u)
    return graph def matrix_to_listlist(weight):
    """transforms a squared weight matrix in a adjacency table of type listlist
    encoding the directed graph corresponding to the entries of the matrix
    different from None

    :param weight: squared weight matrix, weight[u][v] != None iff arc (u, v)
    exists
    :complexity: linear
    :returns: the unweighted directed graph in the listlist representation,
                       listlist[u] contains all v for which arc (u,v) exists.
    """
    graph = [[] for _ in range(len(weight))]
    for u, _ in enumerate(graph):
        for v in range(len(graph)):
            if weight[u][v] is not None:
                graph[u].append(v)
    return graph def listlist_and_matrix_to_listdict(graph, weight=None):
    """Transforms the weighted adjacency list representation of a graph
    of type listlist + optional weight matrix
    into the listdict representation

    :param graph: in listlist representation
    :param weight: optional weight matrix
    :returns: graph in listdict representation
    :complexity: linear
    """
    if weight:
        return [{v: weight[u][v] for v in graph[u]} for u in range(len(graph))]
    else:
        return [{v: None for v in graph[u]} for u in range(len(graph))] def listdict_to_listlist_and_matrix(sparse):
    """Transforms the adjacency list representation of a graph
    of type listdict into the listlist + weight matrix representation

    :param sparse: graph in listdict representation
    :returns: couple with listlist representation, and weight matrix
    :complexity: linear
    """
    V = range(len(sparse))
    graph = [[] for _ in V]
    weight = [[None for v in V] for u in V]
    for u in V:
        for v in sparse[u]:
            graph[u].append(v)
            weight[u][v] = sparse[u][v]
    return graph, weight def dictdict_to_listdict(dictgraph):
    """Transforms a dict-dict graph representation into a
    adjacency dictionary representation (list-dict)

    :param dictgraph: dictionary mapping vertices to dictionary
           such that dictgraph[u][v] is weight of arc (u,v)
    :complexity: linear
    :returns: tuple with graph (listdict), name_to_node (dict),
    node_to_name (list)
    """
    n = len(dictgraph)                            # vertices
    node_to_name = list(dictgraph.keys())         # bijection indices <-> names
    node_to_name.sort()                           # to make it more readable
    name_to_node = {}
    for i in range(n):
        name_to_node[node_to_name[i]] = i
    sparse = [{} for _ in range(n)]               # build sparse graph
    for u in dictgraph:
        for v in dictgraph[u]:
            sparse[name_to_node[u]][name_to_node[v]] = dictgraph[u][v]
    return sparse, name_to_node, node_to_name def anagrams(S):                   # S is a set of strings
    """group a set of words into anagrams

    :param S: set of strings
    :returns: list of lists of strings

    :complexity:
        :math:`O(n k log k)` in average, for n words of length at most k.
        :math:`O(n^2 k log k)` in worst case due to the usage of a dictionary.
    """
    d = {}                         # maps s to list of words with signature s
    for word in S:                 # group words according to the signature
        s = ''.join(sorted(word))  # calculate the signature
        if s in d:
            d[s].append(word)      # append a word to an existing signature
        else:
            d[s] = [word]          # add a new signature and its first word
    # -- extract anagrams, ingoring anagram groups of size 1
    return [d[s] for s in d if len(d[s]) > 1] def arithm_expr_eval(cell, expr):
    """Evaluates a given expression

    :param expr: expression
    :param cell: dictionary variable name -> expression

    :returns: numerical value of expression

    :complexity: linear
    """
    if isinstance(expr, tuple):
        (left, operand, right) = expr
        lval = arithm_expr_eval(cell, left)
        rval = arithm_expr_eval(cell, right)
        if operand == '+':
            return lval + rval
        if operand == '-':
            return lval - rval
        if operand == '*':
            return lval * rval
        if operand == '/':
            return lval // rval
    elif isinstance(expr, int):
        return expr
    else:
        cell[expr] = arithm_expr_eval(cell, cell[expr])
        return cell[expr] def arithm_expr_parse(line_tokens):
    """Constructs an arithmetic expression tree

    :param line_tokens: list of token strings containing the expression
    :returns: expression tree

    :complexity: linear
    """
    vals = []
    ops = []
    for tok in line_tokens + [';']:
        if tok in PRIORITY:  # tok is an operator
            while (tok != '(' and ops and
                   PRIORITY[ops[-1]] >= PRIORITY[tok]):
                right = vals.pop()
                left = vals.pop()
                vals.append((left, ops.pop(), right))
            if tok == ')':
                ops.pop()    # this is the corresponding '('
            else:
                ops.append(tok)
        elif tok.isdigit():  # tok is an integer
            vals.append(int(tok))
        else:                # tok is an identifier
            vals.append(tok)
    return vals.pop() def arithm_expr_target(x, target):
    """ Create arithmetic expression approaching target value
    :param x: allowed constants
    :param target: target value
    :returns: string in form 'expression=value'
    :complexity: huge
    """
    n = len(x)
    expr = [{} for _ in range(1 << n)]
    # expr[S][val]
    # = string solely composed of values in set S that evaluates to val
    for i in range(n):
        expr[1 << i] = {x[i]: str(x[i])}   # store singletons
    all_ = (1 << n) - 1
    for S in range(3, all_ + 1):  # 3: first num that isn't a power of 2
        if expr[S] != {}:
            continue            # in that case S is a power of 2
        for L in range(1, S):   # decompose set S into non-empty sets L, R
            if L & S == L:
                R = S ^ L
                for vL in expr[L]:         # combine expressions from L
                    for vR in expr[R]:     # with expressions from R
                        eL = expr[L][vL]
                        eR = expr[R][vR]
                        expr[S][vL] = eL
                        if vL > vR:    # difference cannot become negative
                            expr[S][vL - vR] = "(%s-%s)" % (eL, eR)
                        if L < R:      # break symmetry
                            expr[S][vL + vR] = "(%s+%s)" % (eL, eR)
                            expr[S][vL * vR] = "(%s*%s)" % (eL, eR)
                        if vR != 0 and vL % vR == 0:  # only integer div
                            expr[S][vL // vR] = "(%s/%s)" % (eL, eR)
    # look for the closest expression from the target
    for dist in range(target + 1):
        for sign in [-1, +1]:
            val = target + sign * dist
            if val in expr[all_]:
                return "%s=%i" % (expr[all_][val], val)
    # never reaches here if x contains integers between 0 and target
    pass def bellman_ford(graph, weight, source=0):
    """ Single source shortest paths by Bellman-Ford

    :param graph: directed graph in listlist or listdict format
    :param weight: can be negative.
                   in matrix format or same listdict graph
    :returns: distance table, precedence table, bool
    :explanation: bool is True if a negative circuit is
                  reachable from the source, circuits
                  can have length 2.
    :complexity: `O(|V|*|E|)`
    """
    n = len(graph)
    dist = [float('inf')] * n
    prec = [None] * n
    dist[source] = 0
    for nb_iterations in range(n):
        changed = False
        for node in range(n):
            for neighbor in graph[node]:
                alt = dist[node] + weight[node][neighbor]
                if alt < dist[neighbor]:
                    dist[neighbor] = alt
                    prec[neighbor] = node
                    changed = True
        if not changed:                   # fixed point
            return dist, prec, False
    return dist, prec, True def bellman_ford2(graph, weight, source):
    """ Single source shortest paths by Bellman-Ford

    :param graph: directed graph in listlist or listdict format
    :param weight: can be negative.
                   in matrix format or same listdict graph
    :returns: distance table, precedence table, bool
    :explanation: bool is true if there is a negative cycle 
                  reachable from the source.
                  distance[v] is +inf if vertex v is not reachable 
                  from source and -inf if there are paths from the 
                  source to v of arbitrary small weight.
    :complexity: `O(|V|*|E|)`
    """
    n = len(graph)
    dist = [float('inf')] * n
    prec = [None] * n
    dist[source] = 0

    def relax():
        for nb_iterations in  range(n-1):
            for node in range(n):
                for neighbor in graph[node]:
                    alt = dist[node] + weight[node][neighbor]
                    if alt < dist[neighbor]:
                        dist[neighbor] = alt
                        prec[neighbor] = node

    relax()
    intermediate = dist[:]  # is fixpoint in absence of neg cycles
    relax()
    for node in range(n):
        if dist[node] < intermediate[node]:
            dist[node] = float('-inf')
    return dist, prec, min(dist) == float('-inf') def knapsack(p, v, cmax):
    """Knapsack problem: select maximum value set of items if total size not
    more than capacity

    :param p: table with size of items
    :param v: table with value of items
    :param cmax: capacity of bag
    :requires: number of items non-zero
    :returns: value optimal solution, list of item indexes in solution
    :complexity: O(n * cmax), for n = number of items
    """
    n = len(p)
    opt = [[0] * (cmax + 1) for _ in range(n + 1)]
    sel = [[False] * (cmax + 1) for _ in range(n + 1)]
    #                               --- basic case
    for cap in range(p[0], cmax + 1):
        opt[0][cap] = v[0]
        sel[0][cap] = True
    #                               --- induction case
    for i in range(1, n):
        for cap in range(cmax + 1):
            if cap >= p[i] and opt[i-1][cap - p[i]] + v[i] > opt[i-1][cap]:
                opt[i][cap] = opt[i-1][cap - p[i]] + v[i]
                sel[i][cap] = True
            else:
                opt[i][cap] = opt[i-1][cap]
                sel[i][cap] = False
    #                               --- reading solution
    cap = cmax
    solution = []
    for i in range(n-1, -1, -1):
        if sel[i][cap]:
            solution.append(i)
            cap -= p[i]
    return (opt[n - 1][cmax], solution) def knapsack2(p, v, cmax):
    """Knapsack problem: select maximum value set of items if total size not
    more than capacity.
    alternative implementation with same behavior.

    :param p: table with size of items
    :param v: table with value of items
    :param cmax: capacity of bag
    :requires: number of items non-zero
    :returns: value optimal solution, list of item indexes in solution
    :complexity: O(n * cmax), for n = number of items
    """
    n = len(p)
    # Plus grande valeur obtenable avec objets ≤ i et capacité c
    pgv = [[0] * (cmax + 1) for _ in range(n)]
    for c in range(cmax + 1):  # Initialisation
        pgv[0][c] = v[0] if c >= p[0] else 0
    pred = {}  # Prédécesseurs pour mémoriser les choix faits
    for i in range(1, n):
        for c in range(cmax + 1):
            pgv[i][c] = pgv[i - 1][c]  # Si on ne prend pas l'objet i
            pred[(i, c)] = (i - 1, c)
            # Est-ce que prendre l'objet i est préférable ?
            if c >= p[i] and pgv[i - 1][c - p[i]] + v[i] > pgv[i][c]:
                pgv[i][c] = pgv[i - 1][c - p[i]] + v[i]
                pred[(i, c)] = (i - 1, c - p[i])  # On marque le prédécesseur
    # On pourrait s'arrêter là, mais si on veut un sous-ensemble d'objets
    # optimal, il faut remonter les marquages
    cursor = (n - 1, cmax)
    chosen = []
    while cursor in pred:
        # Si la case prédécesseur a une capacité inférieure
        if pred[cursor][1] < cursor[1]:
            # C'est qu'on a ramassé l'objet sur le chemin
            chosen.append(cursor[0])
        cursor = pred[cursor]
    if cursor[1] > 0:  # A-t-on pris le premier objet ?
        # (La première ligne n'a pas de prédécesseur.)
        chosen.append(cursor[0])
    return pgv[n - 1][cmax], chosen def bfs(graph, start=0):
    """Shortest path in unweighted graph by BFS

       :param graph: directed graph in listlist or listdict format
       :param int start: source vertex
       :returns: distance table, precedence table
       :complexity: `O(|V|+|E|)`
       """
    to_visit = deque()
    dist = [float('inf')] * len(graph)
    prec = [None] * len(graph)
    dist[start] = 0
    to_visit.appendleft(start)
    while to_visit:              # an empty queue is considered False
        node = to_visit.pop()
        for neighbor in graph[node]:
            if dist[neighbor] == float('inf'):
                dist[neighbor] = dist[node] + 1
                prec[neighbor] = node
                to_visit.appendleft(neighbor)
    return dist, prec def cut_nodes_edges(graph):
    """Bi-connected components

    :param graph: undirected graph. in listlist format. Cannot be in listdict
    format.
    :returns: a tuple with the list of cut-nodes and the list of cut-edges
    :complexity: `O(|V|+|E|)`
    """
    n = len(graph)
    time = 0
    num = [None] * n
    low = [n] * n
    parent = [None] * n        # parent[v] = None if root else parent of v
    critical_children = [0] * n  # cc[u] = #{children v | low[v] >= num[u]}
    times_seen = [-1] * n
    for start in range(n):
        if times_seen[start] == -1:               # init DFS path
            times_seen[start] = 0
            to_visit = [start]
            while to_visit:
                node = to_visit[-1]
                if times_seen[node] == 0:         # start processing
                    num[node] = time
                    time += 1
                    low[node] = float('inf')
                children = graph[node]
                if times_seen[node] == len(children):  # end processing
                    to_visit.pop()
                    up = parent[node]            # propagate low to parent
                    if up is not None:
                        low[up] = min(low[up], low[node])
                        if low[node] >= num[up]:
                            critical_children[up] += 1
                else:
                    child = children[times_seen[node]]   # next arrow
                    times_seen[node] += 1
                    if times_seen[child] == -1:   # not visited yet
                        parent[child] = node      # link arrow
                        times_seen[child] = 0
                        to_visit.append(child)    # (below) back arrow
                    elif num[child] < num[node] and parent[node] != child:
                        low[node] = min(low[node], num[child])
    cut_edges = []
    cut_nodes = []                                # extract solution
    for node in range(n):
        if parent[node] is None:                  # characteristics
            if critical_children[node] >= 2:
                cut_nodes.append(node)
        else:                                     # internal nodes
            if critical_children[node] >= 1:
                cut_nodes.append(node)
            if low[node] >= num[node]:
                cut_edges.append((parent[node], node))
    return cut_nodes, cut_edges def cut_nodes_edges2(graph):
    """Bi-connected components, alternative recursive implementation

    :param graph: undirected graph. in listlist format. Cannot be in listdict
    format.
    :assumes: graph has about 5000 vertices at most, otherwise memory limit is
    reached
    :returns: a tuple with the list of cut-nodes and the list of cut-edges
    :complexity: `O(|V|+|E|)` in average, `O(|V|+|E|^2)` in worst case due to
    use of dictionary
    """
    N = len(graph)
    assert N <= 5000
    recursionlimit = getrecursionlimit()
    setrecursionlimit(max(recursionlimit, N + 42))
    edges = set((i, j) for i in range(N) for j in graph[i] if i <= j)
    nodes = set()
    NOT = -2  # not visited yet; -1 would be buggy `marked[v] != prof - 1`
    FIN = -3  # already visited
    marked = [NOT] * N  # if >= 0, it means depth within the DFS

    # pylint: disable=inconsistent-return-statements
    def DFS(n, prof=0):
        """
        Recursively search graph, update edge list and returns the first
        node the first edge within search to which we can come back.
        """
        if marked[n] == FIN:
            return  # only when there are several connected components
        if marked[n] != NOT:
            return marked[n]
        marked[n] = prof
        m = float('inf')
        count = 0  # useful only for prof == 0
        for v in graph[n]:
            if marked[v] != FIN and marked[v] != prof - 1:
                count += 1
                r = DFS(v, prof+1)
                if r <= prof:
                    edges.discard(tuple(sorted((n, v))))
                if prof and r >= prof:  # only if we are not at root
                    nodes.add(n)
                m = min(m, r)
        # root is an articulation point iff it has more than 2 childs
        if prof == 0 and count >= 2:
            nodes.add(n)
        marked[n] = FIN
        return m
    for r in range(N):
        DFS(r)  # we can count connected components by nb += DFS(r)
    setrecursionlimit(recursionlimit)
    return nodes, edges def continuous_binary_search(f, lo, hi, gap=1e-4):
    """Binary search for a function

    :param f: boolean monotone function with f(hi) = True
    :param int lo:
    :param int hi: with hi >= lo
    :param float gap:
    :returns: first value x in [lo,hi] such that f(x),
             x is computed up to some precision
    :complexity: `O(log((hi-lo)/gap))`
    """
    while hi - lo > gap:
        mid = (lo + hi) / 2.0
        if f(mid):
            hi = mid
        else:
            lo = mid
    return lo def discrete_binary_search(tab, lo, hi):
    """Binary search in a table

    :param tab: boolean monotone table with tab[hi] = True
    :param int lo:
    :param int hi: with hi >= lo
    :returns: first index i in [lo,hi] such that tab[i]
    :complexity: `O(log(hi-lo))`
    """
    while lo < hi:
        mid = lo + (hi - lo) // 2
        if tab[mid]:
            hi = mid
        else:
            lo = mid + 1
    return lo def optimized_binary_search(tab, logsize):
    """Binary search in a table using bit operations

    :param tab: boolean monotone table
       of size :math:`2^\\textrm{logsize}`
       with tab[hi] = True
    :param int logsize:
    :returns: first i such that tab[i]
    :complexity: O(logsize)
    """
    hi = (1 << logsize) - 1
    intervalsize = (1 << logsize) >> 1
    while intervalsize > 0:
        if tab[hi ^ intervalsize]:
            hi ^= intervalsize
        intervalsize >>= 1
    return hi def ternary_search(f, lo, hi, gap=1e-10):
    """Ternary maximum search for a bitonic function

    :param f: boolean bitonic function (increasing then decreasing,
    not necessarily strictly)
    :param int lo:
    :param int hi: with hi >= lo
    :param float gap:
    :returns: value x in [lo,hi] maximizing f(x),
             x is computed up to some precision
    :complexity: `O(log((hi-lo)/gap))`
    """
    while hi - lo > gap:
        step = (hi - lo) / 3.
        if f(lo + step) < f(lo + 2 * step):
            lo += step
        else:
            hi -= step
    return lo def max_bipartite_matching(bigraph):
    """Bipartie maximum matching

    :param bigraph: adjacency list, index = vertex in U,
                                    value = neighbor list in V
    :assumption: U = V = {0, 1, 2, ..., n - 1} for n = len(bigraph)
    :returns: matching list, match[v] == u iff (u, v) in matching
    :complexity: `O(|V|*|E|)`
    """
    n = len(bigraph)               # same domain for U and V
    match = [None] * n
    for u in range(n):
        augment(u, bigraph, [False] * n, match)
    return match def max_bipartite_matching2(bigraph):
    """Bipartie maximum matching

    :param bigraph: adjacency list, index = vertex in U,
                                    value = neighbor list in V
    :comment: U and V can have different cardinalities
    :returns: matching list, match[v] == u iff (u, v) in matching
    :complexity: `O(|V|*|E|)`
    """
    nU = len(bigraph)
    # the following line works only in Python version ≥ 2.5
    # nV = max(max(adjlist, default=-1) for adjlist in bigraph) + 1
    nV = 0
    for adjlist in bigraph:
        for v in adjlist:
            if v + 1 > nV:
                nV = v + 1
    match = [None] * nV
    for u in range(nU):
        augment(u, bigraph, [False] * nV, match)
    return match def bipartite_vertex_cover(bigraph):
    """Bipartite minimum vertex cover by Koenig's theorem

    :param bigraph: adjacency list, index = vertex in U,
                                    value = neighbor list in V
    :assumption: U = V = {0, 1, 2, ..., n - 1} for n = len(bigraph)
    :returns: boolean table for U, boolean table for V
    :comment: selected vertices form a minimum vertex cover,
              i.e. every edge is adjacent to at least one selected vertex
              and number of selected vertices is minimum
    :complexity: `O(|V|*|E|)`
    """
    V = range(len(bigraph))
    matchV = max_bipartite_matching(bigraph)
    matchU = [None for u in V]
    for v in V:                      # -- build the mapping from U to V
        if matchV[v] is not None:
            matchU[matchV[v]] = v
    visitU = [False for u in V]      # -- build max alternating forest
    visitV = [False for v in V]
    for u in V:
        if matchU[u] is None:        # -- starting with free vertices in U
            _alternate(u, bigraph, visitU, visitV, matchV)
    inverse = [not b for b in visitU]
    return (inverse, visitV) def closest_points(S):
    """Closest pair of points

    :param S: list of points
    :requires: size of S at least 2
    :modifies: changes the order in S
    :returns: pair of points p,q from S with minimum Euclidean distance
    :complexity: expected linear time
    """
    shuffle(S)
    assert len(S) >= 2
    p = S[0]                # start with distance between
    q = S[1]                # first two points
    d = dist(p, q)
    while d > 0:            # distance 0 cannot be improved
        r = improve(S, d)
        if r:               # distance improved
            d, p, q = r
        else:               # r is None: could not improve
            break
    return p, q def closest_values(L):
    """Closest values

    :param L: list of values
    :returns: two values from L with minimal distance
    :modifies: the order of L
    :complexity: O(n log n), for n=len(L)
    """
    assert len(L) >= 2
    L.sort()
    valmin, argmin = min((L[i] - L[i - 1], i) for i in range(1, len(L)))
    return L[argmin - 1], L[argmin] def left_turn(a, b, c):
    """function left-turn"""
    return ((a[0] - c[0]) * (b[1] - c[1]) -
            (a[1] - c[1]) * (b[0] - c[0]) > 0) def andrew(S):
    """Convex hull by Andrew

    :param S: list of points as coordinate pairs
    :requires: S has at least 2 points
    :returns: list of points of the convex hull
    :complexity: `O(n log n)`
    """
    S.sort()
    top = []
    bot = []
    for p in S:
        while len(top) >= 2 and not left_turn(p, top[-1], top[-2]):
            top.pop()
        top.append(p)
        while len(bot) >= 2 and not left_turn(bot[-2], bot[-1], p):
            bot.pop()
        bot.append(p)
    return bot[:-1] + top[:0:-1] def dancing_links(size_universe, sets):
    """Exact set cover by the dancing links algorithm

    :param size_universe: universe = {0, 1, ..., size_universe - 1}
    :param sets: list of sets
    :returns: list of set indices partitioning the universe, or None
    :complexity: huge
    """
    header = Cell(None, None, 0, None)  # building the cell structure
    col = []
    for j in range(size_universe):
        col.append(Cell(header, None, 0, None))
    for i, _ in enumerate(sets):
        row = None
        for j in sets[i]:
            col[j].S += 1               # one more entry in this column
            row = Cell(row, col[j], i, col[j])
    sol = []
    if solve(header, sol):
        return sol
    return None def find_cycle(graph):
    """find a cycle in an undirected graph

    :param graph: undirected graph in listlist or listdict format
    :returns: list of vertices in a cycle or None
    :complexity: `O(|V|+|E|)`
    """
    n = len(graph)
    prec = [None] * n  # ancestor marks for visited vertices
    for u in range(n):
        if prec[u] is None:  # unvisited vertex
            S = [u]  # start new DFS
            prec[u] = u  # mark root (not necessary for this algorithm)
            while S:
                u = S.pop()
                for v in graph[u]:  # for all neighbors
                    if v != prec[u]:  # except arcs to father in DFS tree
                        if prec[v] is not None:
                            cycle = [v, u]  # cycle found, (u,v) back edge
                            while u not in (prec[v], prec[u]):  # directed
                                u = prec[u]  # climb up the tree
                                cycle.append(u)
                            return cycle
                        else:
                            prec[v] = u  # v is new vertex in tree
                            S.append(v)
    return None def dfs_recursive(graph, node, seen):
    """DFS, detect connected component, recursive implementation

    :param graph: directed graph in listlist or listdict format
    :param int node: to start graph exploration
    :param boolean-table seen: will be set true for the connected component
          containing node.
    :complexity: `O(|V|+|E|)`
    """
    seen[node] = True
    for neighbor in graph[node]:
        if not seen[neighbor]:
            dfs_recursive(graph, neighbor, seen) def dfs_iterative(graph, start, seen):
    """DFS, detect connected component, iterative implementation

    :param graph: directed graph in listlist or listdict format
    :param int node: to start graph exploration
    :param boolean-table seen: will be set true for the connected component
          containing node.
    :complexity: `O(|V|+|E|)`
    """
    seen[start] = True
    to_visit = [start]
    while to_visit:
        node = to_visit.pop()
        for neighbor in graph[node]:
            if not seen[neighbor]:
                seen[neighbor] = True
                to_visit.append(neighbor) def dfs_grid(grid, i, j, mark='X', free='.'):
    """DFS on a grid, mark connected component, iterative version

    :param grid: matrix, 4-neighborhood
    :param i,j: cell in this matrix, start of DFS exploration
    :param free: symbol for walkable cells
    :param mark: symbol to overwrite visited vertices
    :complexity: linear
    """
    height = len(grid)
    width = len(grid[0])
    to_visit = [(i, j)]
    grid[i][j] = mark
    while to_visit:
        i1, j1 = to_visit.pop()
        for i2, j2 in [(i1 + 1, j1), (i1, j1 + 1),
                       (i1 - 1, j1), (i1, j1 - 1)]:
            if (0 <= i2 < height and 0 <= j2 < width and
                    grid[i2][j2] == free):
                grid[i2][j2] = mark  # mark path
                to_visit.append((i2, j2)) def dfs_grid_recursive(grid, i, j, mark='X', free='.'):
    """DFS on a grid, mark connected component, recursive version

    :param grid: matrix, 4-neighborhood
    :param i,j: cell in this matrix, start of DFS exploration
    :param free: symbol for walkable cells
    :param mark: symbol to overwrite visited vertices
    :complexity: linear
    """
    height = len(grid)
    width = len(grid[0])
    grid[i][j] = mark              # mark path
    for ni, nj in [(i + 1, j), (i, j + 1),
                   (i - 1, j), (i, j - 1)]:
        if 0 <= ni < height and 0 <= nj < width:
            if grid[ni][nj] == free:
                dfs_grid(grid, ni, nj) def dijkstra_update_heap(graph, weight, source=0, target=None):
    """single source shortest paths by Dijkstra
       with a heap implementing item updates

       :param graph: adjacency list or adjacency dictionary of a directed graph
       :param weight: matrix or adjacency dictionary
       :assumes: weights are non-negatif and weights are infinite for non edges
       :param source: source vertex
       :type source: int
       :param target: if given, stops once distance to target found
       :type target: int
       :returns: distance table, precedence table
       :complexity: `O(|V| + |E|log|V|)`
    """
    n = len(graph)
    assert all(weight[u][v] >= 0 for u in range(n) for v in graph[u])
    prec = [None] * n
    dist = [float('inf')] * n
    dist[source] = 0
    heap = OurHeap([(dist[node], node) for node in range(n)])
    while heap:
        dist_node, node = heap.pop()       # Closest node from source
        if node == target:
            break
        for neighbor in graph[node]:
            old = dist[neighbor]
            new = dist_node + weight[node][neighbor]
            if new < old:
                dist[neighbor] = new
                prec[neighbor] = node
                heap.update((old, neighbor), (new, neighbor))
    return dist, prec def dijkstra(graph, weight, source=0, target=None):
    """single source shortest paths by Dijkstra

       :param graph: directed graph in listlist or listdict format
       :param weight: in matrix format or same listdict graph
       :assumes: weights are non-negative
       :param source: source vertex
       :type source: int
       :param target: if given, stops once distance to target found
       :type target: int

       :returns: distance table, precedence table
       :complexity: `O(|V| + |E|log|V|)`
    """
    n = len(graph)
    assert all(weight[u][v] >= 0 for u in range(n) for v in graph[u])
    prec = [None] * n
    black = [False] * n
    dist = [float('inf')] * n
    dist[source] = 0
    heap = [(0, source)]
    while heap:
        dist_node, node = heappop(heap)       # Closest node from source
        if not black[node]:
            black[node] = True
            if node == target:
                break
            for neighbor in graph[node]:
                dist_neighbor = dist_node + weight[node][neighbor]
                if dist_neighbor < dist[neighbor]:
                    dist[neighbor] = dist_neighbor
                    prec[neighbor] = node
                    heappush(heap, (dist_neighbor, neighbor))
    return dist, prec def dilworth(graph):
    """Decompose a DAG into a minimum number of chains by Dilworth

    :param graph: directed graph in listlist or listdict format
    :assumes: graph is acyclic
    :returns: table giving for each vertex the number of its chains
    :complexity: same as matching
    """
    n = len(graph)
    match = max_bipartite_matching(graph)  # maximum matching
    part = [None] * n                      # partition into chains
    nb_chains = 0
    for v in range(n - 1, -1, -1):         # in inverse topological order
        if part[v] is None:                # start of chain
            u = v
            while u is not None:           # follow the chain
                part[u] = nb_chains        # mark
                u = match[u]
            nb_chains += 1
    return part def dinic(graph, capacity, source, target):
    """Maximum flow by Dinic

    :param graph: directed graph in listlist or listdict format
    :param capacity: in matrix format or same listdict graph
    :param int source: vertex
    :param int target: vertex
    :returns: skew symmetric flow matrix, flow value
    :complexity: :math:`O(|V|^2 |E|)`
    """
    assert source != target
    add_reverse_arcs(graph, capacity)
    Q = deque()
    total = 0
    n = len(graph)
    flow = [[0] * n for u in range(n)]   # flow initially empty
    while True:                   # repeat while we can increase
        Q.appendleft(source)
        level = [None] * n          # build levels, None = inaccessible
        level[source] = 0           # by BFS
        while Q:
            u = Q.pop()
            for v in graph[u]:
                if level[v] is None and capacity[u][v] > flow[u][v]:
                    level[v] = level[u] + 1
                    Q.appendleft(v)

        if level[target] is None:   # stop if sink is not reachable
            return flow, total
        up_bound = sum(capacity[source][v] for v in graph[source]) - total
        total += _dinic_step(graph, capacity, level, flow, source, target,
                             up_bound) def dist_grid(grid, source, target=None):
    """Distances in a grid by BFS

    :param grid: matrix with 4-neighborhood
    :param (int,int) source: pair of row, column ind_ices
    :param (int,int) target: exploration stops if target is reached
    :complexity: linear in grid size
    """
    rows = len(grid)
    cols = len(grid[0])
    i, j = source
    grid[i][j] = 's'
    Q = deque()
    Q.append(source)
    while Q:
        i1, j1 = Q.popleft()
        for di, dj, symbol in [(0, +1, '>'),
                               (0, -1, '<'),
                               (+1, 0, 'v'),
                               (-1, 0, '^')]:   # explore all directions
            i2 = i1 + di
            j2 = j1 + dj
            if not (0 <= i2 < rows and 0 <= j2 < cols):
                continue              # reached the bounds of the grid
            if grid[i2][j2] != ' ':   # inaccessible or already visited
                continue
            grid[i2][j2] = symbol     # mark visit
            if (i2, j2) == target:
                grid[i2][j2] = 't'    # goal is reached
                return
            Q.append((i2, j2)) def edmonds_karp(graph, capacity, source, target):
    """Maximum flow by Edmonds-Karp

    :param graph: directed graph in listlist or listdict format
    :param capacity: in matrix format or same listdict graph
    :param int source: vertex
    :param int target: vertex
    :returns: flow matrix, flow value
    :complexity: :math:`O(|V|*|E|^2)`
    """
    add_reverse_arcs(graph, capacity)
    V = range(len(graph))
    flow = [[0 for v in V] for u in V]
    while True:
        augm_path, delta = _augment(graph, capacity, flow, source, target)
        if delta == 0:
            break
        v = target                    # go back to source
        while v != source:
            u = augm_path[v]          # augment flow
            flow[u][v] += delta
            flow[v][u] -= delta
            v = u
    return (flow, sum(flow[source]))  # flow network, amount of flow def eulerian_tour_undirected(graph):
    """Eulerian tour on an undirected graph

       :param graph: directed graph in listlist format, cannot be listdict
       :assumes: graph is eulerian
       :returns: eulerian cycle as a vertex list
       :complexity: `O(|V|+|E|)`
    """
    P = []                            # resulting tour
    Q = [0]                           # vertices to be explored, start at 0
    R = []                            # path from start node
    next_ = [0] * len(graph)          # initialize next_ to 0 for each node
    seen = [set() for _ in graph]     # mark backward arcs
    while Q:
        start = Q.pop()               # explore a cycle from start node
        node = start                            # current node on cycle
        while next_[node] < len(graph[node]):   # visit all allowable arcs
            neighbor = graph[node][next_[node]]  # traverse an arc
            next_[node] += 1                     # mark arc traversed
            if neighbor not in seen[node]:      # not yet traversed
                seen[neighbor].add(node)        # mark backward arc
                R.append(neighbor)              # append to path from start
                node = neighbor                 # move on
        while R:
            Q.append(R.pop())         # add to Q the discovered cycle R
        P.append(start)               # resulting path P is extended
    return P def eulerian_tour_directed(graph):
    """Eulerian tour on a directed graph

       :param graph: directed graph in listlist format, cannot be listdict
       :assumes: graph is eulerian
       :returns: eulerian cycle as a vertex list
       :complexity: `O(|V|+|E|)`
    """
    P = []                            # resulting tour
    Q = [0]                           # vertices to be explored, start at 0
    R = []                            # path from start node
    next_ = [0] * len(graph)          # initialize next_ to 0 for each node
    while Q:
        start = Q.pop()               # explore a cycle from start node
        node = start                            # current node on cycle
        while next_[node] < len(graph[node]):   # visit all allowable arcs
            neighbor = graph[node][next_[node]]  # traverse an arc
            next_[node] += 1                     # mark arc traversed
            R.append(neighbor)                  # append to path from start
            node = neighbor                     # move on
        while R:
            Q.append(R.pop())         # add to Q the discovered cycle R
        P.append(start)               # resulting path P is extended
    return P def random_eulerien_graph(n):
    """Generates some random eulerian graph

       :param int n: number of vertices
       :returns: undirected graph in listlist representation
       :complexity: linear
    """
    graphe = [[] for _ in range(n)]
    for v in range(n - 1):
        noeuds = random.sample(range(v + 1, n), random.choice(
            range(0 if len(graphe[v]) % 2 == 0 else 1, (n - v), 2)))
        graphe[v].extend(noeuds)
        for w in graphe[v]:
            if w > v:
                graphe[w].append(v)
    return graphe def is_eulerian_tour(graph, tour):
    """Eulerian tour on an undirected graph

       :param graph: directed graph in listlist format, cannot be listdict
       :param tour: vertex list
       :returns: test if tour is eulerian
       :complexity: `O(|V|*|E|)` under the assumption that
       set membership is in constant time
    """
    m = len(tour)-1
    arcs = set((tour[i], tour[i+1]) for i in range(m))
    if len(arcs) != m:
        return False
    for (u, v) in arcs:
        if v not in graph[u]:
            return False
    return True def fast_exponentiation(a, b, q):
    """Compute (a pow b) % q, alternative shorter implementation

    :param int a b: non negative
    :param int q: positive
    :complexity: O(log b)
    """
    assert a >= 0 and b >= 0 and q >= 1
    result = 1
    while b:
        if b % 2 == 1:
            result = (result * a) % q
        a = (a * a) % q
        b >>= 1
    return result def fast_exponentiation2(a, b, q):
    """Compute (a pow b) % q

    :param int a b: non negative
    :param int q: positive
    :complexity: O(log b)
    """
    assert a >= 0 and b >= 0 and q >= 1
    p = 0               # only for documentation
    p2 = 1              # 2 ** p
    ap2 = a % q         # a ** (2 ** p)
    result = 1
    while b > 0:
        if p2 & b > 0:  # b's binary decomposition contains 2 ** p
            b -= p2
            result = (result * ap2) % q
        p += 1
        p2 *= 2
        ap2 = (ap2 * ap2) % q
    return result class Fenwick:
    """maintains a tree to allow quick updates and queries
    """
    def __init__(self, t):
        """stores a table t and allows updates and queries
        of prefix sums in logarithmic time.

        :param array t: with numerical values
        """
        self.s = [0] * (len(t) + 1)  # create internal storage
        for a, v in enumerate(t):
            self.add(a, v)           # initialize

    # pylint: disable=redefined-builtin
    def prefixSum(self, a):
        """
        :param int a: index in t, negative a will return 0
        :returns: t[0] + ... + t[a]
        """
        i = a + 1                  # internal index starts at 1
        total = 0
        while i > 0:               # loops over neighbors
            total += self.s[i]     # cumulative sum
            i -= (i & -i)          # left neighbor
        return total

    def intervalSum(self, a, b):
        """
        :param int a b: with 0 <= a <= b
        :returns: t[a] + ... + t[b]
        """
        return self.prefixSum(b) - self.prefixSum(a-1)

    def add(self, a, val):
        """
        :param int a: index in t
        :modifies: adds val to t[a]
        """
        i = a + 1                  # internal index starts at 1
        while i < len(self.s):     # loops over parents
            self.s[i] += val       # update node
            i += (i & -i)          # parent

    # variante:
    # pylint: disable=bad-whitespace
    def intervalAdd(self, a, b, val):
        """Variant, adds val to t[a], to t[a + 1] ... and to t[b]

        :param int a b: with 0 <= a <= b < len(t)
        """
        self.add(a,     +val)
        self.add(b + 1, -val)

    def get(self, a):
        """Variant, reads t[a]

        :param int i: negative a will return 0
        """
        return self.prefixSum(a) def floyd_warshall(weight):
    """All pairs shortest paths by Floyd-Warshall

    :param weight: edge weight matrix
    :modifies: weight matrix to contain distances in graph
    :returns: True if there are negative cycles
    :complexity: :math:`O(|V|^3)`
    """
    V = range(len(weight))
    for k in V:
        for u in V:
            for v in V:
                weight[u][v] = min(weight[u][v],
                                   weight[u][k] + weight[k][v])
    for v in V:
        if weight[v][v] < 0:      # negative cycle found
            return True
    return False def floyd_warshall2(weight):
    """All pairs shortest paths by Floyd-Warshall.
    An improved implementation by Pascal-Ortiz

    :param weight: edge weight matrix
    :modifies: weight matrix to contain distances in graph
    :returns: True if there are negative cycles
    :complexity: :math:`O(|V|^3)`
    """
    for k, Wk in enumerate(weight):
        for _, Wu in enumerate(weight):
            for v, Wuv in enumerate(Wu):
                alt = Wu[k] + Wk[v]
                if alt < Wuv:
                    Wu[v] = alt
    for v, Wv in enumerate(weight):
        if Wv[v] < 0:      # negative cycle found
            return True
    return False def ford_fulkerson(graph, capacity, s, t):
    """Maximum flow by Ford-Fulkerson

    :param graph: directed graph in listlist or listdict format
    :param capacity: in matrix format or same listdict graph
    :param int s: source vertex
    :param int t: target vertex

    :returns: flow matrix, flow value
    :complexity: `O(|V|*|E|*max_capacity)`
    """
    add_reverse_arcs(graph, capacity)
    n = len(graph)
    flow = [[0] * n for _ in range(n)]
    INF = float('inf')
    while _augment(graph, capacity, flow, INF, s, t, [False] * n) > 0:
        pass                         # work already done in _augment
    return (flow, sum(flow[s]))      # flow network, amount of flow def gale_shapley(men, women):
    """Stable matching by Gale-Shapley

    :param men: table of size n, men[i] is preference list of women for men i
    :param women: similar
    :returns: matching table, from women to men
    :complexity: :math:`O(n^2)`
    """
    n = len(men)
    assert n == len(women)
    current_suitor = [0] * n
    spouse = [None] * n
    rank = [[0] * n for j in range(n)]  # build rank
    for j in range(n):
        for r in range(n):
            rank[j][women[j][r]] = r
    singles = deque(range(n))  # all men are single and get in the queue
    while singles:
        i = singles.popleft()
        j = men[i][current_suitor[i]]
        current_suitor[i] += 1
        if spouse[j] is None:
            spouse[j] = i
        elif rank[j][spouse[j]] < rank[j][i]:
            singles.append(i)
        else:
            singles.put(spouse[j])  # sorry for spouse[j]
            spouse[j] = i
    return spouse def gauss_jordan(A, x, b):
    """Linear equation system Ax=b by Gauss-Jordan

    :param A: m by n matrix
    :param x: table of size n
    :param b: table of size m
    :modifies: x will contain solution if any
    :returns int:
          0 if no solution,
          1 if solution unique,
          2 otherwise
    :complexity: :math:`O(n^2m)`
    """
    n = len(x)
    m = len(b)
    assert len(A) == m and len(A[0]) == n
    S = []                        # put linear system in a single matrix S
    for i in range(m):
        S.append(A[i][:] + [b[i]])
    S.append(list(range(n)))      # indices in x
    k = diagonalize(S, n, m)
    if k < m:
        for i in range(k, m):
            if not is_zero(S[i][n]):
                return GJ_ZERO_SOLUTIONS
    for j in range(k):
        x[S[m][j]] = S[j][n]
    if k < n:
        for j in range(k, n):
            x[S[m][j]] = 0
        return GJ_SEVERAL_SOLUTIONS
    return GJ_SINGLE_SOLUTION GJ_ZERO_SOLUTIONS = 0 GJ_SINGLE_SOLUTION = 1 GJ_SEVERAL_SOLUTIONS = 2 class Graph:
    def __init__(self):
        self.neighbors = []
        self.name2node = {}
        self.node2name = []
        self.weight = []

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

    def __getitem__(self, v):
        return self.neighbors[v]

    def add_node(self, name):
        assert name not in self.name2node
        self.name2node[name] = len(self.name2node)
        self.node2name.append(name)
        self.neighbors.append([])
        self.weight.append({})
        return self.name2node[name]

    def add_edge(self, name_u, name_v, weight_uv=None):
        self.add_arc(name_u, name_v, weight_uv)
        self.add_arc(name_v, name_u, weight_uv)

    def add_arc(self, name_u, name_v, weight_uv=None):
        u = self.name2node[name_u]
        v = self.name2node[name_v]
        self.neighbors[u].append(v)
        self.weight[u][v] = weight_uv def dist01(graph, weight, source=0, target=None):
    """Shortest path in a 0,1 weighted graph

    :param graph: directed graph in listlist or listdict format
    :param weight: matrix or adjacency dictionary
    :param int source: vertex
    :param target: exploration stops once distance to target is found
    :returns: distance table, predecessor table
    :complexity: `O(|V|+|E|)`
    """
    n = len(graph)
    dist = [float('inf')] * n
    prec = [None] * n
    black = [False] * n
    dist[source] = 0
    gray = deque([source])
    while gray:
        node = gray.pop()
        if black[node]:
            continue
        black[node] = True
        if node == target:
            break
        for neighbor in graph[node]:
            ell = dist[node] + weight[node][neighbor]
            if black[neighbor] or dist[neighbor] <= ell:
                continue
            dist[neighbor] = ell
            prec[neighbor] = node
            if weight[node][neighbor] == 0:
                gray.append(neighbor)
            else:
                gray.appendleft(neighbor)
    return dist, prec def horn_sat(formula):
    """ Solving a HORN Sat formula

    :param formula: list of couple(posvar, negvars).
                    negvars is a list of the negative variables (can be empty)
                    posvar is the positive variable (can be None)
                    Variables can be any hashable objects: integers, strings...
    :returns: None if formula is not satisfiable, else a minimal set of vars
              that have to be set to true in order to satisfy the formula.
    :complexity: linear
    """
    # --- construct data structures
    CLAUSES = range(len(formula))
    score = [0 for c in CLAUSES]  # number of neg vars not yet in solution
    # the unique positive variable of a clause (if any)
    posvar_in_clause = [None for c in CLAUSES]
    # all clauses where a variable appears negatively
    clauses_with_negvar = defaultdict(set)
    for c in CLAUSES:
        posvar, negvars = formula[c]
        score[c] = len(set(negvars))  # do not count twice negative variables
        posvar_in_clause[c] = posvar
        for v in negvars:
            clauses_with_negvar[v].add(c)
    pool = [set() for s in range(max(score) + 1)]   # create the pool
    for c in CLAUSES:
        pool[score[c]].add(c)           # pool[s] = set of clauses with score s

    # --- solve Horn SAT formula
    solution = set()                    # contains all variables set to True
    while pool[0]:
        curr = pool[0].pop()               # arbitrary zero score clause
        v = posvar_in_clause[curr]
        if v is None:                      # formula is not satisfiable
            return None
        if v in solution or curr in clauses_with_negvar[v]:
            continue                       # clause is already satisfied
        solution.add(v)
        for c in clauses_with_negvar[v]:   # update score
            pool[score[c]].remove(c)
            score[c] -= 1
            pool[score[c]].add(c)          # change c to lower score in pool
    return solution def huffman(freq):
    """Huffman code

    :param freq: dictionary with frequencies for each item
    :returns: dictionary with binary code string for each item
    :complexity: O(n log n), where n = len(freq)
    """
    forest = []                              # build forest with singletons
    for item in freq:
        heappush(forest, (freq[item], item))
    while len(forest) > 1:
        (f_l, left) = heappop(forest)        # merge two trees
        (f_r, right) = heappop(forest)
        heappush(forest, (f_l + f_r, [left, right]))
    code = {}                                # build code from unique tree
    extract(code, forest[0][1], [])
    return code def interval_tree(intervals):
    """Construct an interval tree

    :param intervals: list of half-open intervals
                      encoded as value pairs *[left, right)*
    :assumes: intervals are lexicographically ordered
              ``>>> assert intervals == sorted(intervals)``
    :returns: the root of the interval tree
    :complexity: :math:`O(n)`
    """
    if intervals == []:
        return None
    center = intervals[len(intervals) // 2][0]
    L = []
    R = []
    C = []
    for I in intervals:
        if I[1] <= center:
            L.append(I)
        elif center < I[0]:
            R.append(I)
        else:
            C.append(I)
    by_low = sorted((I[0], I) for I in C)
    by_high = sorted((I[1], I) for I in C)
    IL = interval_tree(L)
    IR = interval_tree(R)
    return _Node(center, by_low, by_high, IL, IR) def intervals_containing(t, p):
    """Query the interval tree

    :param t: root of the interval tree
    :param p: value
    :returns: a list of intervals containing p
    :complexity: O(log n + m), where n is the number of intervals in t,
                and m the length of the returned list
    """
    INF = float('inf')
    if t is None:
        return []
    if p < t.center:
        retval = intervals_containing(t.left, p)
        j = bisect_right(t.by_low, (p, (INF, INF)))
        for i in range(j):
            retval.append(t.by_low[i][1])
    else:
        retval = intervals_containing(t.right, p)
        i = bisect_right(t.by_high, (p, (INF, INF)))
        for j in range(i, len(t.by_high)):
            retval.append(t.by_high[j][1])
    return retval def interval_cover(I):
    """Minimum interval cover

    :param I: list of closed intervals
    :returns: minimum list of points covering all intervals
    :complexity: O(n log n)
    """
    S = []
    # sort by right endpoints
    for start, end in sorted(I, key=lambda v: v[1]):
        if not S or S[-1] < start:
            S.append(end)
    return S def intervals_union(S):
    """Union of intervals

    :param S: list of pairs (low, high) defining intervals [low, high)
    :returns: ordered list of disjoint intervals with the same union as S
    :complexity: O(n log n)
    """
    E = [(low, -1) for (low, high) in S]
    E += [(high, +1) for (low, high) in S]
    nb_open = 0
    last = None
    retval = []
    for x, _dir in sorted(E):
        if _dir == -1:
            if nb_open == 0:
                last = x
            nb_open += 1
        else:
            nb_open -= 1
            if nb_open == 0:
                retval.append((last, x))
    return retval def maximum_border_length(w):
    """Maximum string borders by Knuth-Morris-Pratt

    :param w: string
    :returns: table f such that f[i] is the longest border length of w[:i + 1]
    :complexity: linear
    """
    n = len(w)
    f = [0] * n                # init f[0] = 0
    k = 0                      # current longest border length
    for i in range(1, n):      # compute f[i]
        while w[k] != w[i] and k > 0:
            k = f[k - 1]       # mismatch: try the next border
        if w[k] == w[i]:       # last characters match
            k += 1             # we can increment the border length
        f[i] = k               # we found the maximal border of w[:i + 1]
    return f def knuth_morris_pratt(s, t):
    """Find a substring by Knuth-Morris-Pratt

    :param s: the haystack string
    :param t: the needle string
    :returns: index i such that s[i: i + len(t)] == t, or -1
    :complexity: O(len(s) + len(t))
    """
    sep = '\x00'                   # special unused character
    assert sep not in t and sep not in s
    f = maximum_border_length(t + sep + s)
    n = len(t)
    for i, fi in enumerate(f):
        if fi == n:                # found a border of the length of t
            return i - 2 * n       # beginning of the border in s
    return -1 def powerstring_by_border(u):
    """Power string by Knuth-Morris-Pratt

    :param x: string
    :returns: largest k such that there is a string y with x = y^k
    :complexity: O(len(x))
    """
    f = maximum_border_length(u)
    n = len(u)
    if n % (n - f[-1]) == 0:       # does the alignment shift divide n ?
        return n // (n - f[-1])    # we found a power decomposition
    return 1 def powerstring_by_find(u):
    """Power string using the python find method

    :param x: string
    :returns: largest k such that there is a string y with x = y^k
    :complexity: O(len(x))
    """
    return len(x) // (x + x).find(x, 1) def kruskal(graph, weight):
    """Minimum spanning tree by Kruskal

    :param graph: undirected graph in listlist or listdict format
    :param weight: in matrix format or same listdict graph
    :returns: list of edges of the tree
    :complexity: ``O(|E|log|E|)``
    """
    u_f = UnionFind(len(graph))
    edges = []
    for u, _ in enumerate(graph):
        for v in graph[u]:
            edges.append((weight[u][v], u, v))
    edges.sort()
    min_span_tree = []
    for w_idx, u_idx, v_idx in edges:
        if u_f.union(u_idx, v_idx):
            min_span_tree.append((u_idx, v_idx))
    return min_span_tree def kuhn_munkres(G):      # maximum profit bipartite matching in O(n^4)
    """Maximum profit perfect matching

    for minimum cost perfect matching just inverse the weights

    :param G: squared weight matrix of a complete bipartite graph
    :complexity: :math:`O(n^4)`
    """
    assert len(G) == len(G[0])
    n = len(G)
    mu = [None] * n                 # Empty matching
    mv = [None] * n
    lu = [max(row) for row in G]    # Trivial labels
    lv = [0] * n
    for u0 in range(n):
        if mu[u0] is None:          # Free node
            while True:
                au = [False] * n    # Empty alternating tree
                av = [False] * n
                if improve_matching(G, u0, mu, mv, au, av, lu, lv):
                    break
                improve_labels(G, au, av, lu, lv)
    return (mu, sum(lu) + sum(lv)) def kuhn_munkres(G, TOLERANCE=1e-6):
    """Maximum profit bipartite matching by Kuhn-Munkres

    :param G: weight matrix where G[u][v] is the weight of edge (u,v),
    :param TOLERANCE: a value with absolute value below tolerance
                      is considered as being zero.
                      If G consists of integer or fractional values
                      then TOLERANCE can be chosen 0.
    :requires: graph (U,V,E) is complete bi-partite graph with len(U) <= len(V)
               float('-inf') or float('inf') entries in G
               are allowed but not None.
    :returns: matching table from U to V, value of matching
    :complexity: :math:`O(|U|^2 |V|)`
    """
    nU = len(G)
    U = range(nU)
    nV = len(G[0])
    V = range(nV)
    assert nU <= nV
    mu = [None] * nU                # empty matching
    mv = [None] * nV
    lu = [max(row) for row in G]    # trivial labels
    lv = [0] * nV
    for root in U:                  # build an alternate tree
        au = [False] * nU           # au, av mark nodes...
        au[root] = True             # ... covered by the tree
        Av = [None] * nV            # Av[v] successor of v in the tree
        # for every vertex u, slack[u] := (val, v) such that
        # val is the smallest slack on the constraints (*)
        # with fixed u and v being the corresponding vertex
        slack = [(lu[root] + lv[v] - G[root][v], root) for v in V]
        while True:
            (delta, u), v = min((slack[v], v) for v in V if Av[v] is None)
            assert au[u]
            if delta > TOLERANCE:   # tree is full
                for u0 in U:        # improve labels
                    if au[u0]:
                        lu[u0] -= delta
                for v0 in V:
                    if Av[v0] is not None:
                        lv[v0] += delta
                    else:
                        (val, arg) = slack[v0]
                        slack[v0] = (val - delta, arg)
            assert abs(lu[u] + lv[v] - G[u][v]) <= TOLERANCE  # equality
            Av[v] = u                # add (u, v) to A
            if mv[v] is None:
                break                # alternating path found
            u1 = mv[v]
            assert not au[u1]
            au[u1] = True            # add (u1, v) to A
            for v1 in V:
                if Av[v1] is None:   # update margins
                    alt = (lu[u1] + lv[v1] - G[u1][v1], u1)
                    if slack[v1] > alt:
                        slack[v1] = alt
        while v is not None:         # ... alternating path found
            u = Av[v]                # along path to root
            prec = mu[u]
            mv[v] = u                # augment matching
            mu[u] = v
            v = prec
    return (mu, sum(lu) + sum(lv)) def rabin_karp_matching(s, t):
    """Find a substring by Rabin-Karp

    :param s: the haystack string
    :param t: the needle string

    :returns: index i such that s[i: i + len(t)] == t, or -1
    :complexity: O(len(s) + len(t)) in expected time,
                and O(len(s) * len(t)) in worst case
    """
    hash_s = 0
    hash_t = 0
    len_s = len(s)
    len_t = len(t)
    last_pos = pow(DOMAIN, len_t - 1) % PRIME
    if len_s < len_t:              # substring too long
        return -1
    for i in range(len_t):         # preprocessing
        hash_s = (DOMAIN * hash_s + ord(s[i])) % PRIME
        hash_t = (DOMAIN * hash_t + ord(t[i])) % PRIME
    for i in range(len_s - len_t + 1):
        if hash_s == hash_t:                # hashes match
            # check character by character
            if matches(s, t, i, 0, len_t):
                return i
        if i < len_s - len_t:
            # shift window and calculate new hash on s
            hash_s = roll_hash(hash_s, ord(s[i]), ord(s[i + len_t]),
                               last_pos)
    return -1                               # no match def roman2int(s):
    """Decode roman number

    :param s: string representing a roman number between 1 and 9999
    :returns: the decoded roman number
    :complexity: linear (if that makes sense for constant bounded input size)
    """
    val = 0
    pos10 = 1000
    beg = 0
    for pos in range(3, -1, -1):
        for digit in range(9, -1, -1):
            r = roman[pos][digit]
            if s.startswith(r, beg):  # footnote 1
                beg += len(r)
                val += digit * pos10
                break
        pos10 //= 10
    return val def int2roman(val):
    """Code roman number

    :param val: integer between 1 and 9999
    :returns: the corresponding roman number
    :complexity: linear (if that makes sense for constant bounded input size)
    """
    s = ''
    pos10 = 1000
    for pos in range(3, -1, -1):
        digit = val // pos10
        s += roman[pos][digit]
        val %= pos10
        pos10 //= 10
    return s def laser_mirrors(rows, cols, mir):
    """Orienting mirrors to allow reachability by laser beam

    :param int rows:
    :param int cols: rows and cols are the dimension of the grid
    :param mir: list of mirror coordinates, except
                mir[0]= laser entrance,
                mir[-1]= laser exit.
    :complexity: :math:`O(2^n)`
    """
    # build structures
    n = len(mir)
    orien = [None] * (n + 2)
    orien[n] = 0      # arbitrary orientations
    orien[n + 1] = 0
    succ = [[None for direc in range(4)] for i in range(n + 2)]
    L = [(mir[i][0], mir[i][1], i) for i in range(n)]
    L.append((0, -1, n))                  # enter
    L.append((0, cols, n + 1))            # exit
    last_r, last_i = None, None
    for (r, c, i) in sorted(L):           # sweep by row
        if last_r == r:
            succ[i][LEFT] = last_i
            succ[last_i][RIGHT] = i
        last_r, last_i = r, i
    last_c = None
    for (r, c, i) in sorted(L, key=lambda rci: (rci[1], rci[0])):
        if last_c == c:                   # sweep by column
            succ[i][UP] = last_i
            succ[last_i][DOWN] = i
        last_c, last_i = c, i
    if solve(succ, orien, n, RIGHT):      # exploration
        return orien[:n]
    return None def left_right_inversions(tab):
    """ Compute left and right inversions of each element of a table.

    :param tab: list with comparable elements
    :returns: lists left and right. left[j] = the number of
    i<j such that tab[i] > tab[j].
              right[i] = the number of i<j such that tab[i] > tab[j].
    :complexity: `O(n log n)`
    """
    n = len(tab)
    left = [0] * n
    right = [0] * n
    tmp = [None] * n      # temporary table
    rank = list(range(n))
    _merge_sort(tab, tmp, rank, left, right, 0, n)
    return left, right def levenshtein(x, y):
    """Levenshtein edit distance

    :param x:
    :param y: strings
    :returns: distance
    :complexity: `O(|x|*|y|)`
    """
    n = len(x)
    m = len(y)
    # Create the table A
    #   Row 0 and column 0 are initialized as required
    #   The remaining entries will be overwritten during the computation
    A = [[i + j for j in range(m + 1)] for i in range(n + 1)]
    for i in range(n):
        for j in range(m):
            A[i + 1][j + 1] = min(A[i][j + 1] + 1,              # insert
                                  A[i + 1][j] + 1,              # delete
                                  A[i][j] + int(x[i] != y[j]))  # subst.
    return A[n][m] def longest_common_subsequence(x, y):
    """Longest common subsequence

    Dynamic programming

    :param x:
    :param y: x, y are lists or strings
    :returns: longest common subsequence in form of a string
    :complexity: `O(|x|*|y|)`
    """
    n = len(x)
    m = len(y)

    #                      -- compute optimal length
    A = [[0 for j in range(m + 1)] for i in range(n + 1)]
    for i in range(n):
        for j in range(m):
            if x[i] == y[j]:
                A[i + 1][j + 1] = A[i][j] + 1
            else:
                A[i + 1][j + 1] = max(A[i][j + 1],  A[i + 1][j])

    #                      -- extract solution in reverse order
    sol = []
    i, j = n, m
    while A[i][j] > 0:
        if A[i][j] == A[i - 1][j]:
            i -= 1
        elif A[i][j] == A[i][j - 1]:
            j -= 1
        else:
            i -= 1
            j -= 1
            sol.append(x[i])
    return ''.join(sol[::-1])  # reverse the list to obtain the solution def longest_increasing_subsequence(x):
    """Longest increasing subsequence

    :param x: sequence
    :returns: longest strictly increasing subsequence y
    :complexity: `O(|x|*log(|y|))`
    """
    n = len(x)
    p = [None] * n
    h = [None]
    b = [float('-inf')]  # - infinity
    for i in range(n):
        if x[i] > b[-1]:
            p[i] = h[-1]
            h.append(i)
            b.append(x[i])
        else:
            #   -- binary search: b[k - 1] < x[i] <= b[k]
            k = bisect_left(b, x[i])
            h[k] = i
            b[k] = x[i]
            p[i] = h[k - 1]
    # extract solution in reverse order
    q = h[-1]
    s = []
    while q is not None:
        s.append(x[q])
        q = p[q]
    return s[::-1]  # reverse the list to obtain the solution class LowestCommonAncestorShortcuts:
    """Lowest common ancestor data structure using shortcuts to ancestors
    """
    def __init__(self, prec):
        """builds the structure from a given tree

        :param prec: father for every node, with prec[0] = 0
        :assumes: prec[node] < node
        :complexity: O(n log n), with n = len(nodes)
        """
        n = len(prec)
        self.level = [None] * n        # build levels
        self.level[0] = 0
        for u in range(1, n):
            self.level[u] = 1 + self.level[prec[u]]
        depth = log2ceil(max(self.level[u] for u in range(n))) + 1
        self.anc = [[0] * n for _ in range(depth)]
        for u in range(n):
            self.anc[0][u] = prec[u]
        for k in range(1, depth):
            for u in range(n):
                self.anc[k][u] = self.anc[k - 1][self.anc[k - 1][u]]

    def query(self, u, v):
        """:returns: the lowest common ancestor of u and v
        :complexity: O(log n)
        """
        # -- assume w.l.o.g. that v is not higher than u in the tree
        if self.level[u] > self.level[v]:
            u, v = v, u
        # -- put v at the same level as u
        depth = len(self.anc)
        for k in range(depth-1, -1, -1):
            if self.level[u] <= self.level[v] - (1 << k):
                v = self.anc[k][v]
        assert self.level[u] == self.level[v]
        if u == v:
            return u
        # -- climb until the lowest common ancestor
        for k in range(depth-1, -1, -1):
            if self.anc[k][u] != self.anc[k][v]:
                u = self.anc[k][u]
                v = self.anc[k][v]
        assert self.anc[0][u] == self.anc[0][v]
        return self.anc[0][u] class LowestCommonAncestorRMQ:
    """Lowest common ancestor data structure using a reduction to
       range minimum query
    """
    def __init__(self, graph):
        """builds the structure from a given tree

        :param graph: adjacency matrix of a tree
        :complexity: O(n log n), with n = len(graph)
        """
        n = len(graph)
        dfs_trace = []
        self.last = [None] * n
        to_visit = [(0, 0, None)]            # node 0 is root
        succ = [0] * n
        while to_visit:
            level, node, father = to_visit[-1]
            self.last[node] = len(dfs_trace)
            dfs_trace.append((level, node))
            if succ[node] < len(graph[node]) and \
               graph[node][succ[node]] == father:
                succ[node] += 1
            if succ[node] == len(graph[node]):
                to_visit.pop()
            else:
                neighbor = graph[node][succ[node]]
                succ[node] += 1
                to_visit.append((level + 1, neighbor, node))
        self.rmq = RangeMinQuery(dfs_trace, (float('inf'), None))

    def query(self, u, v):
        """:returns: the lowest common ancestor of u and v
        :complexity: O(log n)
        """
        lu = self.last[u]
        lv = self.last[v]
        if lu > lv:
            lu, lv = lv, lu
        return self.rmq.range_min(lu, lv + 1)[1] def majority(L):
    # snip}
    """Majority

    :param L: list of elements
    :returns: element that appears most in L,
             tie breaking with smallest element
    :complexity: :math:`O(nk)` in average,
                 where n = len(L) and k = max(w for w in L)
                 :math:`O(n^2k)` in worst case due to the use of a dictionary
    """
    assert L    # majority is undefined on the empty set
    # snip{
    count = {}
    for word in L:
        if word in count:
            count[word] += 1
        else:
            count[word] = 1
    # Using min() like this gives the first word with
    #   maximal count "for free"
    val_1st_max, arg_1st_max = min((-count[word], word) for word in count)
    return arg_1st_max def manacher(s):
    """Longest palindrome in a string by Manacher

    :param s: string
    :requires: s is not empty
    :returns: i,j such that s[i:j] is the longest palindrome in s
    :complexity: O(len(s))
    """
    assert set.isdisjoint({'$', '^', '#'}, s)  # Forbidden letters
    if s == "":
        return (0, 1)
    t = "^#" + "#".join(s) + "#$"
    c = 1
    d = 1
    p = [0] * len(t)
    for i in range(2, len(t) - 1):
        #                        -- reflect index i with respect to c
        mirror = 2 * c - i         # = c - (i-c)
        p[i] = max(0, min(d - i, p[mirror]))
        #                        -- grow palindrome centered in i
        while t[i + 1 + p[i]] == t[i - 1 - p[i]]:
            p[i] += 1
        #                        -- adjust center if necessary
        if i + p[i] > d:
            c = i
            d = i + p[i]
    (k, i) = max((p[i], i) for i in range(1, len(t) - 1))
    return ((i - k) // 2, (i + k) // 2)  # extract solution def matrix_mult_opt_order(M):
    """Matrix chain multiplication optimal order

    :param M: list of matrices
    :returns: matrices opt, arg, such that opt[i][j] is the optimal number of
              operations to compute M[i] * ... * M[j] when done in the order
              (M[i] * ... * M[k]) * (M[k + 1] * ... * M[j]) for k = arg[i][j]
    :complexity: :math:`O(n^2)`
    """
    n = len(M)
    r = [len(Mi) for Mi in M]
    c = [len(Mi[0]) for Mi in M]
    opt = [[0 for j in range(n)] for i in range(n)]
    arg = [[None for j in range(n)] for i in range(n)]
    for j_i in range(1, n):   # loop on i, j of increasing j - i = j_i
        for i in range(n - j_i):
            j = i + j_i
            opt[i][j] = float('inf')
            for k in range(i, j):
                alt = opt[i][k] + opt[k + 1][j] + r[i] * c[k] * c[j]
                if opt[i][j] > alt:
                    opt[i][j] = alt
                    arg[i][j] = k
    return opt, arg def matrix_chain_mult(M):
    """Matrix chain multiplication

    :param M: list of matrices
    :returns: M[0] * ... * M[-1], computed in time optimal order
    :complexity: whatever is needed by the multiplications
    """
    opt, arg = matrix_mult_opt_order(M)
    return _apply_order(M, arg, 0, len(M)-1) def max_interval_intersec(S):
    """determine a value that is contained in a largest number
    of given intervals

    :param S: list of half open intervals
    :complexity: O(n log n), where n = len(S)
    """
    B = ([(left,  +1) for left, right in S] +
         [(right, -1) for left, right in S])
    B.sort()
    c = 0
    best = (c, None)
    for x, d in B:
        c += d
        if best[0] < c:
            best = (c, x)
    return best def merge(x, y):
    """Merge two ordered lists

    :param x:
    :param y: x, y are nondecreasing ordered lists
    :returns: union of x and y in order
    :complexity: linear
    """
    z = []
    i = 0
    j = 0
    while i < len(x) or j < len(y):
        if j == len(y) or i < len(x) and x[i] <= y[j]:  # priority on x
            z.append(x[i])
            i += 1
        else:
            z.append(y[j])                              # now switch to y
            j += 1
    return z def min_mean_cycle(graph, weight, start=0):
    """Minimum mean cycle by Karp

    :param graph: directed graph in listlist or listdict format
    :param weight: in matrix format or same listdict graph
    :param int start: vertex that should be contained in cycle
    :returns: cycle as vertex list, average arc weights
              or None if there is no cycle from start
    :complexity:  `O(|V|*|E|)`
    """
    INF = float('inf')
    n = len(graph)                  # compute distances
    dist = [[INF] * n]
    prec = [[None] * n]
    dist[0][start] = 0
    for ell in range(1, n + 1):
        dist.append([INF] * n)
        prec.append([None] * n)
        for node in range(n):
            for neighbor in graph[node]:
                alt = dist[ell - 1][node] + weight[node][neighbor]
                if alt < dist[ell][neighbor]:
                    dist[ell][neighbor] = alt
                    prec[ell][neighbor] = node
    #                               -- find the optimal value
    valmin = INF
    argmin = None
    for node in range(n):
        valmax = -INF
        argmax = None
        for k in range(n):
            alt = (dist[n][node] - dist[k][node]) / float(n - k)
            # do not divide by float(n-k) => cycle of minimal total weight
            if alt >= valmax:     # with >= we get simple cycles
                valmax = alt
                argmax = k
        if argmax is not None and valmax < valmin:
            valmin = valmax
            argmin = (node, argmax)
    #                               -- extract cycle
    if valmin == INF:             # -- there is no cycle
        return None
    C = []
    node, k = argmin
    for l in range(n, k, -1):
        C.append(node)
        node = prec[l][node]
    return C[::-1], valmin class OurHeap:
    """ min heap

    * heap: is the actual heap, heap[1] = index of the smallest element
    * rank: inverse of heap with rank[x]=i iff heap[i]=x
    * n: size of the heap

    :complexity: init O(n log n), len O(1),
                other operations O(log n) in expectation
                and O(n) in worst case, due to the usage of a dictionary
    """
    def __init__(self, items):
        self.heap = [None]  # index 0 will be ignored
        self.rank = {}
        for x in items:
            self.push(x)

    def __len__(self):
        return len(self.heap) - 1

    def push(self, x):
        """Insert new element x in the heap.
           Assumption: x is not already in the heap"""
        assert x not in self.rank
        i = len(self.heap)
        self.heap.append(x)    # add a new leaf
        self.rank[x] = i
        self.up(i)             # maintain heap order

    def pop(self):
        """Remove and return smallest element"""
        root = self.heap[1]
        del self.rank[root]
        x = self.heap.pop()    # remove last leaf
        if self:               # if heap is not empty
            self.heap[1] = x    # move the last leaf
            self.rank[x] = 1    # to the root
            self.down(1)        # maintain heap order
        return root
    # snip}

    # snip{ our_heap_up_down
    def up(self, i):
        """The value of heap[i] has decreased. Maintain heap invariant."""
        x = self.heap[i]
        while i > 1 and x < self.heap[i // 2]:
            self.heap[i] = self.heap[i // 2]
            self.rank[self.heap[i // 2]] = i
            i //= 2
        self.heap[i] = x       # insertion index found
        self.rank[x] = i

    def down(self, i):
        """the value of heap[i] has increased. Maintain heap invariant."""
        x = self.heap[i]
        n = len(self.heap)
        while True:
            left = 2 * i       # climb down the tree
            right = left + 1
            if (right < n and self.heap[right] < x and
                    self.heap[right] < self.heap[left]):
                self.heap[i] = self.heap[right]
                self.rank[self.heap[right]] = i   # move right child up
                i = right
            elif left < n and self.heap[left] < x:
                self.heap[i] = self.heap[left]
                self.rank[self.heap[left]] = i    # move left child up
                i = left
            else:
                self.heap[i] = x   # insertion index found
                self.rank[x] = i
                return

    def update(self, old, new):
        """Replace an element in the heap
        """
        i = self.rank[old]     # change value at index i
        del self.rank[old]
        self.heap[i] = new
        self.rank[new] = i
        if old < new:          # maintain heap order
            self.down(i)
        else:
            self.up(i) class OurQueue:
    """A FIFO queue

    Complexity:
        all operators in amortized constant time,
        except __str__ which is linear
    """
    def __init__(self):
        self.in_stack = []        # tail
        self.out_stack = []       # head

    def __len__(self):
        return len(self.in_stack) + len(self.out_stack)

    def push(self, obj):
        self.in_stack.append(obj)

    def pop(self):
        if not self.out_stack:    # head is empty
            # Note that the in_stack is assigned to the out_stack
            #  in reverse order. This is because the in_stack stores
            #  elements from oldest to newest whereas the out_stack
            #  needs to pop elements from newest to oldest
            self.out_stack = self.in_stack[::-1]
            self.in_stack = []
        return self.out_stack.pop()
# snip}

    def __str__(self):
        return str(self.out_stack[::-1] + self.in_stack) def permutation_rank(p):
    """Given a permutation of {0,..,n-1} find its rank according to
    lexicographical order

       :param p: list of length n containing all integers from 0 to n-1
       :returns: rank between 0 and n! -1
       :beware: computation with big numbers
       :complexity: `O(n^2)`
    """
    n = len(p)
    fact = 1                                # compute (n-1) factorial
    for i in range(2, n):
        fact *= i
    r = 0                                   # compute rank of p
    digits = list(range(n))                 # all yet unused digits
    for i in range(n-1):                    # for all digits except last one
        q = digits.index(p[i])
        r += fact * q
        del digits[q]                       # remove this digit p[i]
        fact //= (n - 1 - i)                # weight of next digit
    return r def rank_permutation(r, n):
    """Given r and n find the permutation of {0,..,n-1} with rank according to
    lexicographical order equal to r

       :param r n: integers with 0 ≤ r < n!
       :returns: permutation p as a list of n integers
       :beware: computation with big numbers
       :complexity: `O(n^2)`
    """
    fact = 1                               # compute (n-1) factorial
    for i in range(2, n):
        fact *= i
    digits = list(range(n))                # all yet unused digits
    p = []                                 # build permutation
    for i in range(n):
        q = r // fact                      # by decomposing r = q * fact + rest
        r %= fact
        p.append(digits[q])
        del digits[q]                      # remove digit at position q
        if i != n - 1:
            fact //= (n - 1 - i)           # weight of next digit
    return p class PartitionRefinement:
    """This data structure implements an order preserving
    partition with refinements.
    """

    def __init__(self, n):
        """Start with the partition consisting of the unique class {0,1,..,n-1}
        complexity: O(n) both in time and space
        """
        c = PartitionClass()   # initially there is a single class of size n
        self.classes = c       # reference to first class in class list
        self.items = [PartitionItem(i, c) for i in range(n)]  # value-ordered

    def refine(self, pivot):
        """Split every class C in the partition into C intersection pivot
        and C setminus pivot complexity: linear in size of pivot
        """
        has_split = []                        # remember which classes split
        for i in pivot:
            if 0 <= i < len(self.items):      # ignore if outside of domain
                x = self.items[i]
                c = x.theclass                # c = class of x
                if not c.split:               # possibly create new split class
                    c.split = PartitionClass(c)
                    if self.classes is c:
                        self.classes = c.split   # always point to 1st class
                    has_split.append(c)
                x.remove()                    # remove from its class
                x.theclass = c.split
                c.split.append(x)             # append to the split class
        for c in has_split:             # clean information about split classes
            c.split = None
            if not c.items:                   # delete class if it became empty
                c.remove()
                del c

    def tolist(self):
        """produce a list representation of the partition
        """
        return [[x.val for x in theclass.items] for theclass in self.classes]

    def order(self):
        """Produce a flatten list of the partition, ordered by classes
        """
        return [x.val for theclass in self.classes for x in theclass.items] def area(p):
    """Area of a polygone

    :param p: list of the points taken in any orientation,
              p[0] can differ from p[-1]
    :returns: area
    :complexity: linear
    """
    A = 0
    for i, _ in enumerate(p):
        A += p[i - 1][0] * p[i][1] - p[i][0] * p[i - 1][1]
    return A / 2. def is_simple(polygon):
    """Test if a rectilinear polygon is is_simple

    :param polygon: list of points as (x,y) pairs along the closed polygon
    :returns: True if the segements do not intersect
    :complexity: O(n log n) for n=len(polygon)
    """
    n = len(polygon)
    order = list(range(n))
    order.sort(key=lambda i: polygon[i])         # lexicographic order
    rank_to_y = list(set(p[1] for p in polygon))
    rank_to_y.sort()
    y_to_rank = {y: rank for rank, y in enumerate(rank_to_y)}
    S = RangeMinQuery([0] * len(rank_to_y))      # sweep structure
    last_y = None
    for i in order:
        x, y = polygon[i]
        rank = y_to_rank[y]
        #                             -- type of point
        right_x = max(polygon[i - 1][0], polygon[(i + 1) % n][0])
        left = x < right_x
        below_y = min(polygon[i - 1][1], polygon[(i + 1) % n][1])
        high = y > below_y
        if left:                      # y does not need to be in S yet
            if S[rank]:
                return False          # two horizontal segments intersect
            S[rank] = -1              # add y to S
        else:
            S[rank] = 0               # remove y from S
        if high:
            lo = y_to_rank[below_y]   # check S between [lo + 1, rank - 1]
            if (below_y != last_y or last_y == y or
                    rank - lo >= 2 and S.range_min(lo + 1, rank)):
                return False          # horiz. & vert. segments intersect
        last_y = y                    # remember for next iteration
    return True def predictive_text(dic):
    """Predictive text for mobile phones

    :param dic: associates weights to words from [a-z]*
    :returns: a dictionary associating to words from [2-9]*
             a corresponding word from the dictionary with highest weight
    :complexity: linear in total word length
    """
    # total_weight[p] = total weight of words having prefix p
    total_weight = {}
    for word, weight in dic:
        prefix = ""
        for x in word:
            prefix += x
            if prefix in total_weight:
                total_weight[prefix] += weight
            else:
                total_weight[prefix] = weight
    #   prop[s] = prefix to display for s
    prop = {}
    for prefix in total_weight:
        code = code_word(prefix)
        if (code not in prop
                or total_weight[prop[code]] < total_weight[prefix]):
            prop[code] = prefix
    return prop def propose(prop, seq):
    """wrapper to access a dictionary even for non-present keys"""
    if seq in prop:
        return prop[seq]
    return None def eratosthene(n):
    """Prime numbers by sieve of Eratosthene

    :param n: positive integer
    :assumes: n > 2
    :returns: list of prime numbers <n
    :complexity: O(n loglog n)
    """
    P = [True] * n
    answ = [2]
    for i in range(3, n, 2):
        if P[i]:
            answ.append(i)
            for j in range(2 * i, n, i):
                P[j] = False
    return answ def gries_misra(n):
    """Prime numbers by the sieve of Gries-Misra
    Computes both the list of all prime numbers less than n,
    and a table mapping every integer 2 ≤ x < n to its smallest prime factor

    :param n: positive integer
    :returns: list of prime numbers, and list of prime factors
    :complexity: O(n)
    """
    primes = []
    factor = [0] * n
    for x in range(2, n):
        if not factor[x]:      # no factor found
            factor[x] = x      # meaning x is prime
            primes.append(x)
        for p in primes:       # loop over primes found so far
            if p > factor[x] or p * x >= n:
                break
            factor[p * x] = p  # p is the smallest factor of p * x
    return primes, factor def rabin_karp_factor(s, t, k):
    """Find a common factor by Rabin-Karp

    :param string s: haystack
    :param string t: needle
    :param int k: factor length
    :returns: (i, j) such that s[i:i + k] == t[j:j + k] or None.
              In case of tie, lexicographical minimum (i, j) is returned
    :complexity: O(len(s) + len(t)) in expected time,
                and O(len(s) + len(t) * k) in worst case
    """
    last_pos = pow(DOMAIN, k - 1) % PRIME
    pos = {}
    assert k > 0
    if len(s) < k or len(t) < k:
        return None
    hash_t = 0

    # First calculate hash values of factors of t
    for j in range(k):
        hash_t = (DOMAIN * hash_t + ord(t[j])) % PRIME
    for j in range(len(t) - k + 1):
        # store the start position with the hash value
        if hash_t in pos:
            pos[hash_t].append(j)
        else:
            pos[hash_t] = [j]
        if j < len(t) - k:
            hash_t = roll_hash(hash_t, ord(t[j]), ord(t[j + k]), last_pos)

    hash_s = 0
    # Now check for matching factors in s
    for i in range(k):         # preprocessing
        hash_s = (DOMAIN * hash_s + ord(s[i])) % PRIME
    for i in range(len(s) - k + 1):
        if hash_s in pos:      # is this signature in s?
            for j in pos[hash_s]:
                if matches(s, t, i, j, k):
                    return (i, j)
        if i < len(s) - k:
            hash_s = roll_hash(hash_s, ord(s[i]), ord(s[i + k]), last_pos)
    return None class RangeMinQuery:
    """Range minimum query

    maintains a table t, can read/write items t[i],
    and query range_min(i,k) = min{ t[i], t[i + 1], ..., t[k - 1]}
    :complexity: all operations in O(log n), for n = len(t)
    """
    def __init__(self, t, INF=float('inf')):
        self.INF = INF
        self.N = 1
        while self.N < len(t):                     # find size N
            self.N *= 2
        self.s = [self.INF] * (2 * self.N)
        for i in range(len(t)):                    # store values of t
            self.s[self.N + i] = t[i]              # in the leaf nodes
        for p in range(self.N - 1, 0, -1):         # fill inner nodes
            self.s[p] = min(self.s[2 * p], self.s[2 * p + 1])

    def __getitem__(self, i):
        return self.s[self.N + i]

    def __setitem__(self, i, v):
        """ sets t[i] to v.
            :complexity: O(log len(t))
        """
        p = self.N + i
        self.s[p] = v
        p //= 2                                    # climb up the tree
        while p > 0:                               # update node
            self.s[p] = min(self.s[2 * p], self.s[2 * p + 1])
            p //= 2

    def range_min(self, i, k):
        """:returns:  min{ t[i], t[i + 1], ..., t[k - 1]}
        :complexity: O(log len(t))
        """
        return self._range_min(1, 0, self.N, i, k)

    def _range_min(self, p, start, span, i, k):
        """returns the minimum in t in the indexes [i, k) intersected
           with [start, start + span).
           p is the node associated to the later interval.
        """
        if start + span <= i or k <= start:        # disjoint intervals
            return self.INF
        if i <= start and start + span <= k:       # contained intervals
            return self.s[p]
        left = self._range_min(2 * p, start, span // 2,
                               i, k)
        right = self._range_min(2 * p + 1, start + span // 2, span // 2,
                                i, k)
        return min(left, right) class LazySegmentTree:
    """maintains a tree to allow quick updates and queries on a table.

    This is more general than a Fenwick tree or a tree for MinRangeQuery. Here
    queries and updates act on index ranges. Updates can be set a range to a
    value or add a value to a range. Queries can be max, min and sum over an
    index range. All operations run in time O(log n) for a the table size n.
    The given ranges are in the form [i,j] where i is included and j excluded.
    In the recursive calls, node is the index of a node in the tree, and left,
    right its range. Values can be any numerical values allowing max, min, and
    sum, such as integers, floating point numbers or fractions (from the class
    Fraction). Updates over an empty range is valid and does nothing. Queries
    over an empty range is valid and returns the neutral value -inf, +inf or
    0.

    If the node is cleared, then maxval, minval, sumval represent for each
    node the query responses over the corresponding index ranges.  If the node
    is not clean, it means that lazyset and/or lazyadd contain suspendet
    update instructions for that node. Clearing a node means propagating these
    values to the descents in the subtrees, and updating maxval,minval and
    sumval for that node.
    """
    def __init__(self, tab):
        """stores an integer table tab.
        will be padded to get a table with a size of a power of 2.

        :param array tab: of positive length
        """
        self.N = 1
        while self.N < len(tab):
            self.N *= 2
        self.maxval = [float('-inf')] * 2 * self.N  # init with neutral values
        self.minval = [float('+inf')] * 2 * self.N
        self.sumval = [0] * 2 * self.N
        self.lazyset = [None] * 2 * self.N
        self.lazyadd = [0] * 2 * self.N
        for i, tabi in enumerate(tab):            # initialize with given table
            j = self.N + i
            self.maxval[j] = self.minval[j] = self.sumval[j] = tabi
        for node in range(self.N - 1, 0, -1):
            self._maintain(node)                    # maintain invariant

    def _maintain(self, node):
        """maintains the invariant for the given node
        :promize: the lazy values are None/0 for this node
        """
        # requires node and its direct descends to be clean
        ll = 2 * node
        r = 2 * node + 1
        assert self.lazyset[node] is None
        assert self.lazyadd[node] == 0
        assert self.lazyset[ll] is None
        assert self.lazyadd[ll] == 0
        assert self.lazyset[r] is None
        assert self.lazyadd[r] == 0
        self.maxval[node] = max(self.maxval[ll], self.maxval[r])
        self.minval[node] = min(self.minval[ll], self.minval[r])
        self.sumval[node] = self.sumval[ll] + self.sumval[r]

    def _clear(self, node, left, right):
        """propagates the lazy updates for this node to the subtrees.
        as a result the maxval, minval, sumval values for the node
        are up to date.
        """
        if self.lazyset[node] is not None:  # first do the pending set
            val = self.lazyset[node]
            self.minval[node] = val
            self.maxval[node] = val
            self.sumval[node] = val * (right - left)
            self.lazyset[node] = None
            if left < right - 1:            # not a leaf
                self.lazyset[2 * node] = val  # propagate to direct childs
                self.lazyadd[2 * node] = 0
                self.lazyset[2 * node + 1] = val
                self.lazyadd[2 * node + 1] = 0
        if self.lazyadd[node] != 0:        # then do the pending add
            val = self.lazyadd[node]
            self.minval[node] += val
            self.maxval[node] += val
            self.sumval[node] += val * (right - left)
            self.lazyadd[node] = 0
            if left < right - 1:            # not at a leaf
                self.lazyadd[2 * node] += val  # propagate to direct childs
                self.lazyadd[2 * node + 1] += val

    def add(self, i, j, val):
        self._add(i, j, val, 1, 0, self.N)

    def set(self, i, j, val):
        self._set(i, j, val, 1, 0, self.N)

    def max(self, i, j):
        return self._max(i, j, 1, 0, self.N)

    def min(self, i, j):
        return self._min(i, j, 1, 0, self.N)

    def sum(self, i, j):
        return self._sum(i, j, 1, 0, self.N)

    def _add(self, i, j, val, node, left, right):
        self._clear(node, left, right)
        if j <= left or right <= i:
            return   # disjoint intervals, nothing to do
        if i <= left and right <= j:
            self.lazyadd[node] += val
            self._clear(node, left, right)
        else:
            mid = (right + left) // 2
            self._add(i, j, val, 2 * node, left, mid)
            self._add(i, j, val, 2 * node + 1, mid, right)
            self._maintain(node)

    def _set(self, i, j, val, node, left, right):
        self._clear(node, left, right)
        if j <= left or right <= i:
            return   # disjoint intervals, nothing to do
        if i <= left and right <= j:
            self.lazyset[node] = val
            self.lazyadd[node] = 0
            self._clear(node, left, right)
        else:
            mid = (right + left) // 2
            self._set(i, j, val, 2 * node, left, mid)
            self._set(i, j, val, 2 * node + 1, mid, right)
            self._maintain(node)

    def _max(self, i, j, node, left, right):
        if j <= left or right <= i:
            return float('-inf')   # neutral value for max
        self._clear(node, left, right)
        if i <= left and right <= j:
            return self.maxval[node]
        else:
            mid = (right + left) // 2
            a = self._max(i, j, 2 * node, left, mid)
            b = self._max(i, j, 2 * node + 1, mid, right)
            return max(a, b)

    def _min(self, i, j, node, left, right):
        if j <= left or right <= i:
            return float('+inf')   # neutral value for min
        self._clear(node, left, right)
        if i <= left and right <= j:
            return self.minval[node]
        else:
            mid = (right + left) // 2
            a = self._min(i, j, 2 * node, left, mid)
            b = self._min(i, j, 2 * node + 1, mid, right)
            return min(a, b)

    def _sum(self, i, j, node, left, right):
        if j <= left or right <= i:
            return 0               # neutral value for sum
        self._clear(node, left, right)
        if i <= left and right <= j:
            return self.sumval[node]
        else:
            mid = (right + left) // 2
            a = self._sum(i, j, 2 * node, left, mid)
            b = self._sum(i, j, 2 * node + 1, mid, right)
            return a + b

    def _dump(self):
        f = open("tmp.dot", "w")
        print("digraph G{", file=f)
        print('0 [label="lazyset/lazyadd/maxval/minval/sumval"]', file=f)
        for node in range(1, 2 * self.N):
            s = '%i [label="%s/%i/%s/%s/%s"]' % \
                (node, self.lazyset[node], self.lazyadd[node],
                    self.maxval[node], self.minval[node], self.sumval[node])
            print(s.replace('inf', '∞'), file=f)
        for node in range(1, self.N):
            print("%i -> %i" % (node, 2 * node), file=f)
            print("%i -> %i" % (node, 2 * node + 1), file=f)
        print("}", file=f)
        f.close() def rectangles_from_grid(P, black=1):
    """Largest area rectangle in a binary matrix

    :param P: matrix
    :param black: search for rectangles filled with value black
    :returns: area, left, top, right, bottom of optimal rectangle
             consisting of all (i,j) with
             left <= j < right and top <= i <= bottom
    :complexity: linear
    """
    rows = len(P)
    cols = len(P[0])
    t = [0] * cols
    best = None
    for i in range(rows):
        for j in range(cols):
            if P[i][j] == black:
                t[j] += 1
            else:
                t[j] = 0
        (area, left, height, right) = rectangles_from_histogram(t)
        alt = (area, left, i, right, i-height)
        if best is None or alt > best:
            best = alt
    return best def rectangles_from_histogram(H):
    """Largest Rectangular Area in a Histogram

    :param H: histogram table
    :returns: area, left, height, right, rect. is [0, height] * [left, right)
    :complexity: linear
    """
    best = (float('-inf'), 0, 0, 0)
    S = []
    H2 = H + [float('-inf')]  # extra element to empty the queue
    for right, _ in enumerate(H2):
        x = H2[right]
        left = right
        while len(S) > 0 and S[-1][1] >= x:
            left, height = S.pop()
            # first element is area of candidate
            rect = (height * (right - left), left, height, right)
            if rect > best:
                best = rect
        S.append((left, x))
    return best def rectangles_from_points(S):
    """How many rectangles can be formed from a set of points

    :param S: list of points, as coordinate pairs
    :returns: the number of rectangles
    :complexity: :math:`O(n^2)`
    """
    answ = 0
    pairs = {}
    for j, _ in enumerate(S):
        for i in range(j):      # loop over point pairs (p,q)
            px, py = S[i]
            qx, qy = S[j]
            center = (px + qx, py + qy)
            dist = (px - qx) ** 2 + (py - qy) ** 2
            signature = (center, dist)
            if signature in pairs:
                answ += len(pairs[signature])
                pairs[signature].append((i, j))
            else:
                pairs[signature] = [(i, j)]
    return answ def min_scalar_prod(x, y):
    """Permute vector to minimize scalar product

    :param x:
    :param y: x, y are vectors of same size
    :returns: min sum x[i] * y[sigma[i]] over all permutations sigma
    :complexity: O(n log n)
    """
    x1 = sorted(x)  # make copies to preserve the input arguments
    y1 = sorted(y)
    return sum(x1[i] * y1[-i - 1] for i in range(len(x1))) def shortest_cycle(graph):
    """ Finding a shortest cycle in an undirected unweighted graph

    :param graph: undirected graph in listlist or listdict format
    :returns: None or a list C describing a shortest cycle
    :complexity: `O(|V|*|E|)`
    """
    best_cycle = float('inf')
    best_u = None
    best_v = None
    best_tree = None
    V = list(range(len(graph)))
    for root in V:
        tree, cycle_len, u, v = bfs(graph, root, best_cycle // 2)
        if cycle_len < best_cycle:
            best_cycle = cycle_len
            best_u = u
            best_v = v
            best_tree = tree
    if best_cycle == float('inf'):
        return None                   # no cycle found
    Pu = path(best_tree, best_u)      # combine path to make a cycle
    Pv = path(best_tree, best_v)
    cycle = Pu[::-1] + Pv   # last vertex equals first vertex
    return cycle[1:]        # remove duplicate vertex def powergraph(graph, k):
    """Compute the k-th
       `powergraph <https://en.wikipedia.org/wiki/Graph_power>`_
       which has an edge u,v for every vertex pair
       of distance at most k in the given graph.

    :param graph: undirected graph in listlist or listdict format
    :param k: non-negative integer.
    :complexity: O(V^3)
    """
    V = range(len(graph))
    # create weight matrix for paths of length 1
    M = [[float('inf') for v in V] for u in V]
    for u in V:
        for v in graph[u]:
            M[u][v] = M[v][u] = 1
        M[u][u] = 0
    floyd_warshall(M)
    return [[v for v in V if M[u][v] <= k] for u in V] class SortedSet(AbstractSkipList):

    Node = namedtuple('Node', 'key next count')

    def __init__(self, iterable=()):
        self.head = self.Node(None, [], [])
        self._len = 0
        for i in iterable:
            self.add(i)

    def pop(self):
        """Pops the first element"""
        try:
            x = next(iter(self))
            self.remove(x)
            return x
        except StopIteration:
            raise KeyError('pop from an empty set')
            # raise KeyError('pop from an empty set') from None

    def add(self, key):
        height = self.randomHeight()
        self.insert(self.Node(key, [None] * height, [0] * height))

    def __repr__(self):
        return '{%s}' % (', '.join(map(str, iter(self)))) class SortedDict(AbstractSkipList):

    Node = namedtuple('Node', 'key val next count')

    def __init__(self, iterable=()):
        self.head = self.Node(None, None, [], [])
        self._len = 0
        if hasattr(iterable, 'keys') and hasattr(iterable, '__getitem__'):
            for i in iterable.keys():
                self[i] = iterable[i]
        else:
            for i, j in iterable:
                self[i] = j

    def keys(self):
        return list(self)

    def __getitem__(self, key):
        x = self.find(key)
        if x is None:
            raise KeyError
        else:
            return x.val

    def __setitem__(self, key, value):
        try:
            self.remove(key)
        except KeyError:
            pass
        height = self.randomHeight()
        self.insert(self.Node(key, value, [None] * height, [0] * height))

    def __delitem__(self, key):
        self.remove(key)

    def __repr__(self):
        return '{%s}' % (', '.join('%s: %s' % (k, self[k]) for k in self)) def tarjan(graph):
    """Strongly connected components by Tarjan, iterative implementation

    :param graph: directed graph in listlist format, cannot be listdict
    :returns: list of lists for each component
    :complexity: linear
    """
    n = len(graph)
    dfs_num = [None] * n
    dfs_min = [n] * n
    waiting = []
    waits = [False] * n  # invariant: waits[v] iff v in waiting
    sccp = []          # list of detected components
    dfs_time = 0
    times_seen = [-1] * n
    for start in range(n):
        if times_seen[start] == -1:                    # initiate path
            times_seen[start] = 0
            to_visit = [start]
            while to_visit:
                node = to_visit[-1]                    # top of stack
                if times_seen[node] == 0:              # start process
                    dfs_num[node] = dfs_time
                    dfs_min[node] = dfs_time
                    dfs_time += 1
                    waiting.append(node)
                    waits[node] = True
                children = graph[node]
                if times_seen[node] == len(children):  # end of process
                    to_visit.pop()                     # remove from stack
                    dfs_min[node] = dfs_num[node]      # compute dfs_min
                    for child in children:
                        if waits[child] and dfs_min[child] < dfs_min[node]:
                            dfs_min[node] = dfs_min[child]
                    if dfs_min[node] == dfs_num[node]:  # representative
                        component = []                 # make component
                        while True:                    # add nodes
                            u = waiting.pop()
                            waits[u] = False
                            component.append(u)
                            if u == node:              # until repr.
                                break
                        sccp.append(component)
                else:
                    child = children[times_seen[node]]
                    times_seen[node] += 1
                    if times_seen[child] == -1:        # not visited yet
                        times_seen[child] = 0
                        to_visit.append(child)
    return sccp def kosaraju(graph):
    """Strongly connected components by Kosaraju

    :param graph: directed graph in listlist format, cannot be listdict
    :returns: list of lists for each component
    :complexity: linear
    """
    n = len(graph)
    order = []
    sccp = []
    kosaraju_dfs(graph, range(n), order, [])
    kosaraju_dfs(reverse(graph), order[::-1], [], sccp)
    return sccp[::-1]  # follow inverse topological order def tarjan_recursif(graph):
    """Strongly connected components by Tarjan, recursive implementation

    :param graph: directed graph in listlist format, cannot be listdict
    :returns: list of lists for each component
    :complexity: linear
    """
    global sccp, waiting, dfs_time, dfs_num
    sccp = []
    waiting = []
    waits = [False] * len(graph)
    dfs_time = 0
    dfs_num = [None] * len(graph)

    def dfs(node):
        global sccp, waiting, dfs_time, dfs_num
        waiting.append(node)           # new node is waiting
        waits[node] = True
        dfs_num[node] = dfs_time       # mark visit
        dfs_time += 1
        dfs_min = dfs_num[node]        # compute dfs_min
        for neighbor in graph[node]:
            if dfs_num[neighbor] is None:
                dfs_min = min(dfs_min, dfs(neighbor))
            elif waits[neighbor] and dfs_min > dfs_num[neighbor]:
                dfs_min = dfs_num[neighbor]
        if dfs_min == dfs_num[node]:   # representative of a component
            sccp.append([])            # make a component
            while True:                # add waiting nodes
                u = waiting.pop()
                waits[u] = False
                sccp[-1].append(u)
                if u == node:          # until representative
                    break
        return dfs_min

    for node in range(len(graph)):
        if dfs_num[node] is None:
            dfs(node)
    return sccp def subset_sum(x_table, r_target):
    """Subsetsum by splitting

    :param x_table: table of values
    :param r_target: target value
    :returns bool: if there is a subsequence of x_table with total sum r_target
    :complexity: :math:`O(n^{\\lceil n/2 \\rceil})`
    """
    k = len(x_table) // 2              # divide input
    y_value = list(part_sum(x_table[:k]))
    z_value = [r_target - v for v in part_sum(x_table[k:])]
    y_value.sort()       # test of intersection between y_value and z_value
    z_value.sort()
    i = 0
    j = 0
    while i < len(y_value) and j < len(z_value):
        if y_value[i] == z_value[j]:
            return True
        if y_value[i] < z_value[j]:  # increment index of smallest element
            i += 1
        else:
            j += 1
    return False def subset_sum2(x_table, r_target):
    """Subsetsum by splitting

    :param x_table: table of values
    :param r_target: target value
    :returns bool: if there is a subsequence of x_table with total sum r_target
    :complexity: :math:`O(n^{\\lceil n/2 \\rceil})`
    """
    k = len(x_table) // 2              # divide input
    y_set = part_sum2(x_table[:k])
    z_set = set(r_target - value for value in part_sum2(x_table[k:]))
    return len(y_set & z_set) > 0        # test intersection def subset_sum(x, R):
    """Subsetsum

    :param x: table of non negative values
    :param R: target value
    :returns bool: True if a subset of x sums to R
    :complexity: O(n*R)
    """
    b = [False] * (R + 1)
    b[0] = True
    for xi in x:
        for s in range(R, xi - 1, -1):
            b[s] |= b[s - xi]
    return b[R] def coin_change(x, R):
    """Coin change

    :param x: table of non negative values
    :param R: target value
    :returns bool: True if there is a non negative linear combination
    of x that has value R
    :complexity: O(n*R)
    """
    b = [False] * (R + 1)
    b[0] = True
    for xi in x:
        for s in range(xi, R + 1):
            b[s] |= b[s - xi]
    return b[R] def sudoku(G):
    """Solving Sudoku

    :param G: integer matrix with 0 at empty cells
    :returns bool: True if grid could be solved
    :modifies: G will contain the solution
    :complexity: huge, but linear for usual published 9x9 grids
    """
    global N, N2, N4
    if len(G) == 16:              # for a 16 x 16 sudoku grid
        N, N2, N4 = 4, 16, 256
    e = N * N4
    universe = e + 1
    S = [[rc(a), rv(a), cv(a), bv(a)] for a in range(N4 * N2)]
    A = [e]
    for r in range(N2):
        for c in range(N2):
            if G[r][c] != 0:
                a = assignment(r, c, G[r][c] - 1)
                A += S[a]
    sol = dancing_links(universe, S + [A])
    if sol:
        for a in sol:
            if a < len(S):
                G[row(a)][col(a)] = val(a) + 1
        return True
    return False def three_partition(x):
    """partition a set of integers in 3 parts of same total value

    :param x: table of non negative values
    :returns: triplet of the integers encoding the sets, or None otherwise
    :complexity: :math:`O(2^{2n})`
    """
    f = [0] * (1 << len(x))
    for i, _ in enumerate(x):
        for S in range(1 << i):
            f[S | (1 << i)] = f[S] + x[i]
    for A in range(1 << len(x)):
        for B in range(1 << len(x)):
            if A & B == 0 and f[A] == f[B] and 3 * f[A] == f[-1]:
                return (A, B, ((1 << len(x)) - 1) ^ A ^ B)
    return None def topological_order_dfs(graph):
    """Topological sorting by depth first search

    :param graph: directed graph in listlist format, cannot be listdict
    :returns: list of vertices in order
    :complexity: `O(|V|+|E|)`
    """
    n = len(graph)
    order = []
    times_seen = [-1] * n
    for start in range(n):
        if times_seen[start] == -1:
            times_seen[start] = 0
            to_visit = [start]
            while to_visit:
                node = to_visit[-1]
                children = graph[node]
                if times_seen[node] == len(children):
                    to_visit.pop()
                    order.append(node)
                else:
                    child = children[times_seen[node]]
                    times_seen[node] += 1
                    if times_seen[child] == -1:
                        times_seen[child] = 0
                        to_visit.append(child)
    return order[::-1] def topological_order(graph):
    """Topological sorting by maintaining indegree

    :param graph: directed graph in listlist format, cannot be listdict
    :returns: list of vertices in order
    :complexity: `O(|V|+|E|)`
    """
    V = range(len(graph))
    indeg = [0 for _ in V]
    for node in V:            # compute indegree
        for neighbor in graph[node]:
            indeg[neighbor] += 1
    Q = [node for node in V if indeg[node] == 0]
    order = []
    while Q:
        node = Q.pop()        # node without incoming arrows
        order.append(node)
        for neighbor in graph[node]:
            indeg[neighbor] -= 1
            if indeg[neighbor] == 0:
                Q.append(neighbor)
    return order def Trie(S):  # Build the trie for the words in the dictionary S
    """
    :param S: set of words
    :returns: trie containing all words from S
    :complexity: linear in total word sizes from S
    """
    T = None
    for w in S:
        T = add(T, w)
    return T def spell_check(T, w):  # Spell check a word against the trie
    """Spellchecker

    :param T: trie encoding the dictionary
    :param w: given word
    :returns: a closest word from the dictionary
    :complexity: linear if distance was constant
    """
    assert T is not None
    dist = 0
    while True:
        u = search(T, dist, w)
        if u is not None:  # Match at distance dist
            return u
        dist += 1          # No match - try increasing the distance def two_sat(formula):
    """Solving a 2-SAT boolean formula

    :param formula: list of clauses, a clause is pair of literals
                    over X1,...,Xn for some n.
                    a literal is an integer, for example -1 = not X1, 3 = X3
    :returns: table with boolean assignment satisfying the formula or None
    :complexity: linear
    """
    # num_variables is the number of variables
    num_variables = max(abs(clause[p])
                        for p in (0, 1) for clause in formula)
    graph = [[] for node in range(2 * num_variables)]
    for x_idx, y_idx in formula:                           # x_idx or y_idx
        graph[_vertex(-x_idx)].append(_vertex(y_idx))     # -x_idx => y_idx
        graph[_vertex(-y_idx)].append(_vertex(x_idx))     # -y_idx => x_idx
    sccp = tarjan(graph)
    comp_id = [None] * (2 * num_variables)  # each node's component ID
    assignment = [None] * (2 * num_variables)
    for component in sccp:
        rep = min(component)             # representative of the component
        for vtx in component:
            comp_id[vtx] = rep
            if assignment[vtx] is None:
                assignment[vtx] = True
                assignment[vtx ^ 1] = False    # complementary literal
    for i in range(num_variables):
        if comp_id[2 * i] == comp_id[2 * i + 1]:
            return None                        # insatisfiable formula
    return assignment[::2] def union_rectangles_naive(R):
    """Area of union of rectangles

    :param R: list of rectangles defined by (x1, y1, x2, y2)
       where (x1, y1) is top left corner and (x2, y2) bottom right corner
    :returns: area
    :complexity: :math:`O(n^3)`
    """
    X = set()        # set of all x coordinates in the input
    Y = set()        # same for y
    for x1, y1, x2, y2 in R:
        assert x1 <= x2 and y1 <= y2
        X.add(x1)
        X.add(x2)
        Y.add(y1)
        Y.add(y2)
    j_to_x = list(sorted(X))
    i_to_y = list(sorted(Y))
    # X and Y partition space into a grid
    area = 0
    for j in range(len(j_to_x) - 1):      # loop over columns in grid
        x1 = j_to_x[j]
        x2 = j_to_x[j + 1]
        for i in range(len(i_to_y) - 1):  # loop over rows
            y1 = i_to_y[i]                # (x1,...,y2) is the grid cell
            y2 = i_to_y[i + 1]
            if rectangles_contains_point(R, x1, y1):
                area += (y2 - y1) * (x2 - x1)  # cell is covered
    return area def union_rectangles(R):
    """Area of union of rectangles.

    Sweep from top to bottom.
    Maintain in a set the horizontal projection of rectangles,
    for which the top border has been processed but not yet the bottom.

    :param R: list of rectangles defined by (x1, y1, x2, y2)
       where (x1, y1) is top left corner and (x2, y2) bottom right corner
    :returns: area
    :complexity: :math:`O(n^2 \\log n)`
    """
    events = []
    for x1, y1, x2, y2 in R:                     # initialize events
        assert x1 <= x2 and y1 <= y2
        events.append((y1, OPENING, x1, x2))
        events.append((y2, CLOSING, x1, x2))
    current_intervals = Counter()
    area = 0
    previous_y = 0  # arbitrary initial value,
    #                 ok, because union_intervals is 0 at first event
    for y, offset, x1, x2 in sorted(events):         # sweep top down
        area += (y - previous_y) * union_intervals(current_intervals)
        previous_y = y
        current_intervals[x1, x2] += offset
    return area def union_rectangles_fast(R):
    """Area of union of rectangles

    :param R: list of rectangles defined by (x1, y1, x2, y2)
       where (x1, y1) is top left corner and (x2, y2) bottom right corner
    :returns: area
    :complexity: :math:`O(n^2)`
    """
    X = set()                 # set of all x coordinates in the input
    events = []               # events for the sweep line
    for x1, y1, x2, y2 in R:
        assert x1 <= x2 and y1 <= y2
        X.add(x1)
        X.add(x2)
        events.append((y1, OPENING, x1, x2))
        events.append((y2, CLOSING, x1, x2))
    # array of x coordinates in left to right order
    i_to_x = list(sorted(X))
    # inverse dictionary maps x coordinate to its rank
    x_to_i = {xi: i for i, xi in enumerate(i_to_x)}
    # nb_current_rectangles[i] = number of rectangles intersected
    # by the sweepline in interval [i_to_x[i], i_to_x[i + 1]]
    nb_current_rectangles = [0] * (len(i_to_x) - 1)
    area = 0
    length_union_intervals = 0
    previous_y = 0  # arbitrary initial value,
    #                 because length is 0 at first iteration
    for y, offset, x1, x2 in sorted(events):
        area += (y - previous_y) * length_union_intervals
        i1 = x_to_i[x1]
        i2 = x_to_i[x2]         # update nb_current_rectangles
        for j in range(i1, i2):
            length_interval = i_to_x[j + 1] - i_to_x[j]
            if nb_current_rectangles[j] == 0:
                length_union_intervals += length_interval
            nb_current_rectangles[j] += offset
            if nb_current_rectangles[j] == 0:
                length_union_intervals -= length_interval
        previous_y = y
    return area def union_rectangles_fastest(R):
    """Area of union of rectangles

    :param R: list of rectangles defined by (x1, y1, x2, y2)
       where (x1, y1) is top left corner and (x2, y2) bottom right corner
    :returns: area
    :complexity: :math:`O(n \\log n)`
    """
    if R == []:               # segment tree would fail on an empty list
        return 0
    X = set()                 # set of all x coordinates in the input
    events = []               # events for the sweep line
    for Rj in R:
        (x1, y1, x2, y2) = Rj
        assert x1 <= x2 and y1 <= y2
        X.add(x1)
        X.add(x2)
        events.append((y1, OPENING, x1, x2))
        events.append((y2, CLOSING, x1, x2))
    i_to_x = list(sorted(X))
    # inverse dictionary
    x_to_i = {i_to_x[i]: i for i in range(len(i_to_x))}
    L = [i_to_x[i + 1] - i_to_x[i] for i in range(len(i_to_x) - 1)]
    C = CoverQuery(L)
    area = 0
    previous_y = 0  # arbitrary initial value,
    #                 because C.cover() is 0 at first iteration
    for y, offset, x1, x2 in sorted(events):
        area += (y - previous_y) * C.cover()
        i1 = x_to_i[x1]
        i2 = x_to_i[x2]
        C.change(i1, i2, offset)
        previous_y = y
    return area def windows_k_distinct(x, k):
    """Find all largest windows containing exactly k distinct elements

    :param x: list or string
    :param k: positive integer
    :yields: largest intervals [i, j) with len(set(x[i:j])) == k
    :complexity: `O(|x|)`
    """
    dist, i, j = 0, 0, 0                # dist = |{x[i], ..., x[j-1]}|
    occ = {xi: 0 for xi in x}           # number of occurrences in x[i:j]
    while j < len(x):
        while dist == k:                # move start of interval
            occ[x[i]] -= 1              # update counters
            if occ[x[i]] == 0:
                dist -= 1
            i += 1
        while j < len(x) and (dist < k or occ[x[j]]):
            if occ[x[j]] == 0:          # update counters
                dist += 1
            occ[x[j]] += 1
            j += 1                      # move end of interval
        if dist == k:
            yield (i, j)                # one interval found
import unittest
import random
from collections import deque
from tryalgo.graph import write_graph, extract_path, make_flow_labels
from tryalgo.graph import tree_adj_to_prec, tree_prec_to_adj
from tryalgo.graph import matrix_to_listlist, listlist_and_matrix_to_listdict
from tryalgo.graph import listdict_to_listlist_and_matrix, dictdict_to_listdict
from tryalgo.anagrams import anagrams
from tryalgo.arithm_expr_eval import arithm_expr_eval, arithm_expr_parse
from tryalgo.arithm_expr_target import arithm_expr_target
from tryalgo.bellman_ford import bellman_ford, bellman_ford2
from tryalgo.knapsack import knapsack, knapsack2
from tryalgo.bfs import bfs
from tryalgo.biconnected_components import cut_nodes_edges, cut_nodes_edges2
from tryalgo.binary_search import continuous_binary_search, discrete_binary_search, optimized_binary_search, ternary_search
from tryalgo.bipartite_matching import max_bipartite_matching, max_bipartite_matching2
from tryalgo.bipartite_vertex_cover import bipartite_vertex_cover
from tryalgo.closest_points import closest_points
from tryalgo.closest_values import closest_values
from tryalgo.convex_hull import left_turn, andrew
from tryalgo.dancing_links import dancing_links
from tryalgo.dfs import find_cycle, dfs_recursive, dfs_iterative, dfs_grid, dfs_grid_recursive
from tryalgo.dijkstra import dijkstra_update_heap, dijkstra
from tryalgo.dilworth import dilworth
from tryalgo.dinic import dinic
from tryalgo.dist_grid import dist_grid
from tryalgo.edmonds_karp import edmonds_karp
from tryalgo.eulerian_tour import eulerian_tour_undirected, eulerian_tour_directed, random_eulerien_graph, is_eulerian_tour
from tryalgo.fast_exponentiation import fast_exponentiation, fast_exponentiation2
from tryalgo.fenwick import Fenwick
from tryalgo.floyd_warshall import floyd_warshall, floyd_warshall2
from tryalgo.ford_fulkerson import ford_fulkerson
from tryalgo.gale_shapley import gale_shapley
from tryalgo.gauss_jordan import gauss_jordan, GJ_ZERO_SOLUTIONS, GJ_SINGLE_SOLUTION, GJ_SEVERAL_SOLUTIONS
from tryalgo.graph import Graph
from tryalgo.graph01 import dist01
from tryalgo.horn_sat import horn_sat
from tryalgo.huffman import huffman
from tryalgo.interval_tree import interval_tree, intervals_containing
from tryalgo.interval_cover import interval_cover
from tryalgo.intervals_union import intervals_union
from tryalgo.knuth_morris_pratt import maximum_border_length, knuth_morris_pratt, powerstring_by_border, powerstring_by_find
from tryalgo.kruskal import kruskal
from tryalgo.kuhn_munkres_n4 import kuhn_munkres as kuhn_munkres_n4
from tryalgo.kuhn_munkres import kuhn_munkres as kuhn_munkres_n3
from tryalgo.rabin_karp import rabin_karp_matching
from tryalgo.roman_numbers import roman2int, int2roman
from tryalgo.laser_mirrors import laser_mirrors
from tryalgo.left_right_inversions import left_right_inversions
from tryalgo.levenshtein import levenshtein
from tryalgo.longest_common_subsequence import longest_common_subsequence
from tryalgo.longest_increasing_subsequence import longest_increasing_subsequence
from tryalgo.lowest_common_ancestor import LowestCommonAncestorShortcuts, LowestCommonAncestorRMQ
from tryalgo.majority import majority
from tryalgo.manacher import manacher
from tryalgo.matrix_chain_mult import matrix_mult_opt_order, matrix_chain_mult
from tryalgo.max_interval_intersec import max_interval_intersec
from tryalgo.merge_ordered_lists import merge
from tryalgo.min_mean_cycle import min_mean_cycle
from tryalgo.our_heap import OurHeap
from tryalgo.our_queue import OurQueue
from tryalgo.permutation_rank import permutation_rank, rank_permutation
from tryalgo.partition_refinement import PartitionRefinement
from tryalgo.polygon import area, is_simple
from tryalgo.predictive_text import predictive_text, propose
from tryalgo.primes import eratosthene, gries_misra
from tryalgo.rabin_karp import rabin_karp_factor
from tryalgo.range_minimum_query import RangeMinQuery, LazySegmentTree
from tryalgo.rectangles_from_grid import rectangles_from_grid
from tryalgo.rectangles_from_histogram import rectangles_from_histogram
from tryalgo.rectangles_from_points import rectangles_from_points
from tryalgo.scalar import min_scalar_prod
from tryalgo.shortest_cycle import shortest_cycle, powergraph
from tryalgo.skip_list import SortedSet, SortedDict
from tryalgo.strongly_connected_components import tarjan, kosaraju, tarjan_recursif
from tryalgo.subsetsum_divide import subset_sum as subset_sum1
from tryalgo.subsetsum_divide import subset_sum2 as subset_sum2
from tryalgo.subsetsum import subset_sum as subset_sum3, coin_change
from tryalgo.sudoku import sudoku
from tryalgo.three_partition import three_partition
from tryalgo.topological_order import topological_order_dfs, topological_order
from tryalgo.trie import Trie, spell_check
from tryalgo.two_sat import two_sat
from tryalgo.union_rectangles import union_rectangles_naive, union_rectangles, union_rectangles_fast, union_rectangles_fastest
from tryalgo.windows_k_distinct import windows_k_distinct
            seen = [False] * n
            dfs(graph, start, seen)
            retval = [node for node in range(n) if seen[node]]
            seen = [False] * n
            dfs(listlist_and_matrix_to_listdict(graph), start, seen)
            self.assertEqual(retval, [node for node in range(n) if seen[node]])
            return retval

        n = 100000
        G = [[v for v in range(u + 1, min(u + 10, n))] for u in range(n)]
        self.assertEqual(reachable(G, 3, dfs_iterative), list(range(3, n)))

        # tests sur la profondeur possible de la récursion
        # from random import *
        # from sys    import *

        # G2 = [[] for u in range(n)]
        # for _ in range(15 * n):
        #     u = random.randint(0, n - 1)
        #     v = random.randint(0, n - 1)
        #     G2[u].append(v)
        #     G2[v].append(u)

        # setrecursionlimit(2 * n + 10)
        # graph = G2
        # seen = [None] * n
        # dfs(0)
        # print( "aha")
        # exit(0)

        for f in [dfs_recursive, dfs_iterative]:
            self.assertEqual(reachable([[]], 0, f), [0])
            self.assertEqual(reachable([[1], []], 0, f), [0, 1])
            self.assertEqual(reachable([[1], []], 1, f), [1])
            self.assertEqual(reachable([[1], [2], []], 1, f), [1, 2])
            self.assertEqual(reachable(
                [[1, 5], [2, 3, 5], [3], [4, 5], [5], []], 2, f), [2, 3, 4, 5])

    def test_dfs_grid(self):

        outTextGrid = """\
##########
XXXXX#...#
####X###.#
#..#X#...#
#..#X#.###
###XX#.#X#
#X#X####X#
#XXXXXXXX#
########X#\
"""
        inTextGrid = """\
##########
.....#...#
####.###.#
#..#.#...#
#..#.#.###
###..#.#.#
#.#.####.#
#........#
########.#\
"""
        out = [list(line.strip()) for line in outTextGrid.split('\n')]
        for f in [dfs_grid, dfs_grid_recursive]:
            grid = [list(line.strip()) for line in inTextGrid.split('\n')]
            f(grid, 1, 0)
            self.assertEqual(str(grid), str(out))

    def test_find_cycle(self):
        L = [([], None),
             ([[]], None),
             ([[], []], None),
             ([[1], [0]], None),
             ([[], [2], [1]], None),
             ([[1, 2], [0, 2], [0, 1]], {0, 1, 2}),
             ([[1, 2], [0], [0], [2, 4], [3]], [2, 3]),
             ([[1, 2], [0], [0]], None),
             ([[1, 2], [0], [0], [4, 5], [3, 5], [3, 4]], [4, 5, 3])]
        for graph, result in L:
            for g in [graph, listlist_and_matrix_to_listdict(graph)]:
                answer = find_cycle(g)
                if isinstance(result, set):
                    self.assertEqual(set(answer), result)
                else:
                    self.assertEqual(answer, result)

    def test_dijkstra(self):
        _ = None
        L_dir = [("reachable_cycle",
                  [[1], [2], [5], [6], [3], [4, 8], [7], [4], [2]],
                  #  0  1  2  3  4  5  6  7  8
                  [[_, 1, _, 0, _, _, _, _, _],  # 0
                   [1, 1, 4, _, _, _, _, _, _],  # 1
                   [_, 4, _, _, _, 3, _, _, 2],  # 2
                   [0, _, _, _, 1, _, 3, _, _],  # 3
                   [_, _, _, 1, _, 5, _, 1, _],  # 4
                   [_, _, 3, _, 5, _, _, _, 2],  # 5
                   [_, _, _, 3, _, _, _, 6, _],  # 6
                   [_, _, _, _, 1, _, 6, _, 2],  # 7
                   [_, _, 2, _, _, 2, _, 2, _],  # 8
                   ], [0, 1, 2, 5, 8]),
                 ("single_vertex",
                  [[]],
                  [[_]], [0]),
                 ("single_arc",
                  [[1], []],
                  [[_, 3],
                   [_, _]], [0, 1]),
                 ("two_nodes",
                  [[], []],
                  [[_, _],
                   [_, _]], None),
                 ]
        for f in [dijkstra, dijkstra_update_heap]:
            for title, graph, weight, shortest_path in L_dir:
                sparse = listlist_and_matrix_to_listdict(graph, weight)
                for g, w in [(graph, weight), (sparse, sparse)]:
                    source = 0
                    target = len(g) - 1
                    dist, prec = f(g, w, source, target)

                    path = extract_path(prec, target)