zariski_vankampen.py

# sage.doctest: needs sage.geometry.polyhedron sage.graphs sage.groups sage.rings.number_field
r"""
Zariski-Van Kampen method implementation

This file contains functions to compute the fundamental group of
the complement of a curve in the complex affine or projective plane,
using Zariski-Van Kampen approach. It depends on the package ``sirocco``.

The current implementation allows to compute a presentation of the
fundamental group of curves over the rationals or number fields with
a fixed embedding on `\QQbar`.

Instead of computing a representation of the braid monodromy, we
choose several base points and a system of paths joining them that
generate all the necessary loops around the points of the discriminant.
The group is generated by the free groups over these points, and
braids over these paths give relations between these generators.
This big group presentation is simplified at the end.

AUTHORS:

- Miguel Marco (2015-09-30): Initial version

EXAMPLES::

    sage: # needs sirocco
    sage: from sage.schemes.curves.zariski_vankampen import fundamental_group, braid_monodromy
    sage: R.<x, y> = QQ[]
    sage: f = y^3 + x^3 - 1
    sage: braid_monodromy(f)
    ([s1*s0, s1*s0, s1*s0], {0: 0, 1: 0, 2: 0}, {}, 3)
    sage: fundamental_group(f)
    Finitely presented group < x0 |  >
"""
# ****************************************************************************
#       Copyright (C) 2015 Miguel Marco <mmarco@unizar.es>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#                  https://www.gnu.org/licenses/
# ****************************************************************************
import itertools

from copy import copy
from itertools import combinations

from sage.combinat.permutation import Permutation
from sage.functions.generalized import sign
from sage.geometry.voronoi_diagram import VoronoiDiagram
from sage.graphs.graph import Graph
from sage.groups.braid import BraidGroup
from sage.groups.finitely_presented import wrap_FpGroup
from sage.groups.free_group import FreeGroup
from sage.groups.perm_gps.permgroup_named import SymmetricGroup
from sage.matrix.constructor import matrix
from sage.misc.cachefunc import cached_function
from sage.misc.flatten import flatten
from sage.misc.lazy_import import lazy_import
from sage.misc.misc_c import prod
from sage.parallel.decorate import parallel
from sage.rings.complex_interval_field import ComplexIntervalField
from sage.rings.complex_mpfr import ComplexField
from sage.rings.integer_ring import ZZ
from sage.rings.number_field.number_field import NumberField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.qqbar import QQbar
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
from sage.schemes.curves.constructor import Curve

lazy_import('sage.libs.braiding', ['leftnormalform', 'rightnormalform'])

roots_interval_cache = {}


def braid_from_piecewise(strands):
    r"""
    Compute the braid corresponding to the piecewise linear curves strands.

    INPUT:

    - ``strands`` -- a list of lists of tuples ``(t, c1, c2)``, where ``t``
      is a number between 0 and 1, and ``c1`` and ``c2`` are rationals
      or algebraic reals.

    OUTPUT:

    The braid formed by the piecewise linear strands.

    EXAMPLES::

        sage: # needs sirocco
        sage: from sage.schemes.curves.zariski_vankampen import braid_from_piecewise
        sage: paths = [[(0, 0, 1), (0.2, -1, -0.5), (0.8, -1, 0), (1, 0, -1)],
        ....:          [(0, -1, 0), (0.5, 0, -1), (1, 1, 0)],
        ....:          [(0, 1, 0), (0.5, 1, 1), (1, 0, 1)]]
        sage: braid_from_piecewise(paths)
        s0*s1
    """
    L = strands
    i = min(val[1][0] for val in L)
    totalpoints = [[[a[0][1], a[0][2]]] for a in L]
    indices = [1 for a in range(len(L))]
    while i < 1:
        for j, val in enumerate(L):
            if val[indices[j]][0] > i:
                xauxr = val[indices[j] - 1][1]
                xauxi = val[indices[j] - 1][2]
                yauxr = val[indices[j]][1]
                yauxi = val[indices[j]][2]
                aaux = val[indices[j] - 1][0]
                baux = val[indices[j]][0]
                interpolar = xauxr + (yauxr - xauxr)*(i - aaux) / (baux - aaux)
                interpolai = xauxi + (yauxi - xauxi)*(i - aaux) / (baux - aaux)
                totalpoints[j].append([interpolar, interpolai])
            else:
                totalpoints[j].append([val[indices[j]][1],
                                       val[indices[j]][2]])
                indices[j] = indices[j] + 1
        i = min(val[indices[k]][0] for k, val in enumerate(L))

    for j, val in enumerate(L):
        totalpoints[j].append([val[-1][1], val[-1][2]])
    braid = []
    G = SymmetricGroup(len(totalpoints))

    def sgn(x, y):
        if x < y:
            return 1
        if x > y:
            return -1
        return 0

    for i in range(len(totalpoints[0]) - 1):
        l1 = [totalpoints[j][i] for j in range(len(L))]
        l2 = [totalpoints[j][i + 1] for j in range(len(L))]
        M = [[l1[s], l2[s]] for s in range(len(l1))]
        M.sort()
        l1 = [a[0] for a in M]
        l2 = [a[1] for a in M]
        cruces = []
        for j, l2j in enumerate(l2):
            l1j = l1[j]
            for k in range(j):
                if l2j < l2[k]:
                    t = (l1j[0] - l1[k][0]) / ((l2[k][0] - l2j[0]) + (l1j[0] - l1[k][0]))
                    s = sgn(l1[k][1] * (1 - t) + t * l2[k][1],
                            l1j[1] * (1 - t) + t * l2j[1])
                    cruces.append([t, k, j, s])
        if cruces:
            cruces.sort()
            P = G(Permutation([]))
            while cruces:
                # we select the crosses in the same t
                crucesl = [c for c in cruces if c[0] == cruces[0][0]]
                crossesl = [(P(c[2] + 1) - P(c[1] + 1), c[1], c[2], c[3])
                            for c in crucesl]
                cruces = cruces[len(crucesl):]
                while crossesl:
                    crossesl.sort()
                    c = crossesl.pop(0)
                    braid.append(c[3] * min(map(P, [c[1] + 1, c[2] + 1])))
                    P = G(Permutation([(c[1] + 1, c[2] + 1)])) * P
                    crossesl = [(P(cr[2] + 1) - P(cr[1] + 1),
                                 cr[1], cr[2], cr[3]) for cr in crossesl]

    B = BraidGroup(len(L))
    return B(braid)


def discrim(pols) -> tuple:
    r"""
    Return the points in the discriminant of the product of the polynomials
    of a list or tuple ``pols``.

    The result is the set of values of the first variable for which
    two roots in the second variable coincide.

    INPUT:

    - ``pols`` -- a list or tuple of polynomials in two variables with
      coefficients in a number field with a fixed embedding in `\QQbar`.

    OUTPUT:

    A tuple with the roots of the discriminant in `\QQbar`.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import discrim
        sage: R.<x, y> = QQ[]
        sage: flist = (y^3 + x^3 - 1, 2 * x + y)
        sage: sorted((discrim(flist)))
        [-0.522757958574711?,
         -0.500000000000000? - 0.866025403784439?*I,
         -0.500000000000000? + 0.866025403784439?*I,
         0.2613789792873551? - 0.4527216721561923?*I,
         0.2613789792873551? + 0.4527216721561923?*I,
         1]
    """
    x, y = pols[0].parent().gens()
    field = pols[0].base_ring()
    pol_ring = PolynomialRing(field, (x,))

    @parallel
    def discrim_pairs(f, g):
        if g is None:
            return pol_ring(f.discriminant(y))
        return pol_ring(f.resultant(g, y))

    pairs = [(f, None) for f in pols] + [tuple(t) for t
                                         in combinations(pols, 2)]
    fdiscrim = discrim_pairs(pairs)
    rts = ()
    poly = 1
    for u in fdiscrim:
        h0 = u[1].radical()
        h1 = h0 // h0.gcd(poly)
        rts += tuple(h1.roots(QQbar, multiplicities=False))
        poly = poly * h1
    return rts


@cached_function
def corrected_voronoi_diagram(points):
    r"""
    Compute a Voronoi diagram of a set of points with rational coordinates.
    The given points are granted to lie one in each bounded region.

    INPUT:

    - ``points`` -- a tuple of complex numbers

    OUTPUT:

    A Voronoi diagram constructed from rational approximations of the points,
    with the guarantee that each bounded region contains exactly one of the
    input points.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import corrected_voronoi_diagram
        sage: points = (2, I, 0.000001, 0, 0.000001*I)
        sage: V = corrected_voronoi_diagram(points)
        sage: V
        The Voronoi diagram of 9 points of dimension 2 in the Rational Field
        sage: V.regions()
        {P(-7, 0): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices and 2 rays,
        P(0, -7): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices and 2 rays,
        P(0, 0): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices,
        P(0, 1): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 5 vertices,
        P(0, 1/1000000): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices,
        P(0, 7): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 3 vertices and 2 rays,
        P(1/1000000, 0): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 5 vertices,
        P(2, 0): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 5 vertices,
        P(7, 0): A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 2 vertices and 2 rays}
    """
    prec = 53
    point_coordinates = [(p.real(), p.imag()) for p in points]
    while True:
        RF = RealField(prec)
        apprpoints = {(QQ(RF(p[0])), QQ(RF(p[1]))): p
                      for p in point_coordinates}
        added_points = 3 * max(map(abs, flatten(apprpoints))) + 1
        configuration = list(apprpoints.keys()) + [(added_points, 0),
                                                   (-added_points, 0),
                                                   (0, added_points),
                                                   (0, -added_points)]
        V = VoronoiDiagram(configuration)
        valid = True
        for r in V.regions().items():
            if (not r[1].rays() and
               not r[1].interior_contains(apprpoints[r[0].affine()])):
                prec += 53
                valid = False
                break
        if valid:
            break
    return V


def orient_circuit(circuit, convex=False, precision=53, verbose=False):
    r"""
    Reverse a circuit if it goes clockwise; otherwise leave it unchanged.

    INPUT:

    - ``circuit`` --  a circuit in the graph of a Voronoi Diagram, given
      by a list of edges

    - ``convex`` -- boolean (default: ``False``); if set to ``True`` a simpler
      computation is made

    -  ``precision`` -- bits of precision (default: 53)

    - ``verbose`` -- boolean (default: ``False``); for testing purposes

    OUTPUT:

    The same circuit if it goes counterclockwise, and its reversed otherwise,
    given as the ordered list of vertices with identic extremities.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import orient_circuit
        sage: points = [(-4, 0), (4, 0), (0, 4), (0, -4), (0, 0)]
        sage: V = VoronoiDiagram(points)
        sage: E = Graph()
        sage: for reg  in V.regions().values():
        ....:     if reg.rays() or reg.lines():
        ....:         E  = E.union(reg.vertex_graph())
        sage: E.vertices(sort=True)
        [A vertex at (-2, -2),
         A vertex at (-2, 2),
         A vertex at (2, -2),
         A vertex at (2, 2)]
        sage: cir = E.eulerian_circuit()
        sage: cir
        [(A vertex at (-2, -2), A vertex at (2, -2), None),
         (A vertex at (2, -2), A vertex at (2, 2), None),
         (A vertex at (2, 2), A vertex at (-2, 2), None),
         (A vertex at (-2, 2), A vertex at (-2, -2), None)]
        sage: cir_oriented = orient_circuit(cir); cir_oriented
        (A vertex at (-2, -2), A vertex at (2, -2), A vertex at (2, 2),
         A vertex at (-2, 2), A vertex at (-2, -2))
        sage: cirinv = list(reversed([(c[1],c[0],c[2]) for c in cir]))
        sage: cirinv
        [(A vertex at (-2, -2), A vertex at (-2, 2), None),
         (A vertex at (-2, 2), A vertex at (2, 2), None),
         (A vertex at (2, 2), A vertex at (2, -2), None),
         (A vertex at (2, -2), A vertex at (-2, -2), None)]
        sage: orient_circuit(cirinv) == cir_oriented
        True
        sage: cir_oriented == orient_circuit(cir, convex=True)
        True
        sage: P0=[(1,1/2),(0,1),(1,1)]; P1=[(0,3/2),(-1,0)]
        sage: Q=Polyhedron(P0).vertices()
        sage: Q = [Q[2], Q[0], Q[1]] + [_ for _ in reversed(Polyhedron(P1).vertices())]
        sage: Q
        [A vertex at (1, 1/2), A vertex at (0, 1), A vertex at (1, 1),
         A vertex at (0, 3/2), A vertex at (-1, 0)]
        sage: E = Graph()
        sage: for v, w in zip(Q, Q[1:] + [Q[0]]):
        ....:   E.add_edge((v, w))
        sage: cir = orient_circuit(E.eulerian_circuit(), precision=1, verbose=True)
        2
        sage: cir
        (A vertex at (1, 1/2), A vertex at (0, 1), A vertex at (1, 1),
         A vertex at (0, 3/2), A vertex at (-1, 0), A vertex at (1, 1/2))
    """
    vectors = [v[1].vector() - v[0].vector() for v in circuit]
    circuit_vertex = (circuit[0][0],) + tuple(e[1] for e in circuit)
    circuit_vertex = tuple(circuit_vertex)
    if convex:
        pr = matrix([vectors[0], vectors[1]]).determinant()
        if pr > 0:
            # return circuit
            return circuit_vertex
        elif pr < 0:
            return tuple(reversed(circuit_vertex))
    prec = precision
    while True:
        CIF = ComplexIntervalField(prec)
        totalangle = sum((CIF(*vectors[i]) / CIF(*vectors[i - 1])).argument()
                         for i in range(len(vectors)))
        if totalangle < 0:
            return tuple(reversed(circuit_vertex))
        if totalangle > 0:
            return circuit_vertex
        prec *= 2
        if verbose:
            print(prec)


def voronoi_cells(V, vertical_lines=frozenset()):
    r"""
    Compute the graph, the boundary graph, a base point, a positive orientation
    of the boundary graph, and the dual graph of a corrected Voronoi diagram.

    INPUT:

    - ``V`` -- a corrected Voronoi diagram

    - ``vertical_lines`` -- frozenset (default: ``frozenset()``); indices of the
      vertical lines

    OUTPUT:

    - ``G`` -- the graph of the 1-skeleton of ``V``
    - ``E`` -- the subgraph of the boundary
    - ``p`` -- a vertex in ``E``
    - ``EC`` -- a list of vertices (representing a counterclockwise orientation
      of ``E``) with identical first and last elements)
    - ``DG`` -- the dual graph of ``V``, where the vertices are labelled
      by the compact regions of ``V`` and the edges by their dual edges.
    - ``vertical_regions`` -- dictionary for the regions associated
      with vertical lines

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import corrected_voronoi_diagram, voronoi_cells
        sage: points = (2, I, 0.000001, 0, 0.000001*I)
        sage: V = corrected_voronoi_diagram(points)
        sage: G, E, p, EC, DG, VR = voronoi_cells(V, vertical_lines=frozenset((1,)))
        sage: Gv = G.vertices(sort=True)
        sage: Ge = G.edges(sort=True)
        sage: len(Gv), len(Ge)
        (12, 16)
        sage: Ev = E.vertices(sort=True); Ev
        [A vertex at (-4, 4),
        A vertex at (-49000001/14000000, 1000001/2000000),
        A vertex at (-7/2, -7/2),
        A vertex at (-7/2, 1/2000000),
        A vertex at (1/2000000, -7/2),
        A vertex at (2000001/2000000, -24500001/7000000),
        A vertex at (11/4, 4),
        A vertex at (9/2, -9/2),
        A vertex at (9/2, 9/2)]
        sage: Ev.index(p)
        7
        sage: EC
        (A vertex at (9/2, -9/2),
        A vertex at (9/2, 9/2),
        A vertex at (11/4, 4),
        A vertex at (-4, 4),
        A vertex at (-49000001/14000000, 1000001/2000000),
        A vertex at (-7/2, 1/2000000),
        A vertex at (-7/2, -7/2),
        A vertex at (1/2000000, -7/2),
        A vertex at (2000001/2000000, -24500001/7000000),
        A vertex at (9/2, -9/2))
        sage: len(DG.vertices(sort=True)), len(DG.edges(sort=True))
        (5, 7)
        sage: edg = DG.edges(sort=True)[0]; edg
        ((0,
          (A vertex at (9/2, -9/2),
           A vertex at (9/2, 9/2),
           A vertex at (11/4, 4),
           A vertex at (2000001/2000000, 500001/1000000),
           A vertex at (2000001/2000000, -24500001/7000000),
           A vertex at (9/2, -9/2))),
         (1,
          (A vertex at (-49000001/14000000, 1000001/2000000),
           A vertex at (1000001/2000000, 1000001/2000000),
           A vertex at (2000001/2000000, 500001/1000000),
           A vertex at (11/4, 4),
           A vertex at (-4, 4),
           A vertex at (-49000001/14000000, 1000001/2000000))),
         (A vertex at (2000001/2000000, 500001/1000000), A vertex at (11/4, 4), None))
        sage: edg[-1] in Ge
        True
        sage: VR
        {1: (A vertex at (-49000001/14000000, 1000001/2000000),
             A vertex at (1000001/2000000, 1000001/2000000),
             A vertex at (2000001/2000000, 500001/1000000),
             A vertex at (11/4, 4),
             A vertex at (-4, 4),
             A vertex at (-49000001/14000000, 1000001/2000000))}
    """
    regions = V.regions()
    points = [p for p in V.regions().keys() if V.regions()[p].is_compact()]
    compact_regions = [regions[p] for p in points]
    vertical_regions = {}
    non_compact_regions = [reg for reg in V.regions().values()
                           if not reg.is_compact()]
    G = Graph([u.vertices() for v in compact_regions for u in v.faces(1)],
              format='list_of_edges')
    E = Graph([u.vertices() for v in non_compact_regions for u in v.faces(1)
               if u.is_compact()], format='list_of_edges')
    p = next(E.vertex_iterator())
    EC = orient_circuit(E.eulerian_circuit())
    DG = Graph()
    for i, reg in enumerate(compact_regions):
        Greg0 = orient_circuit(reg.graph().eulerian_circuit(), convex=True)
        if i in vertical_lines:
            vertical_regions[i] = Greg0
        DG.add_vertex((i, Greg0))
    for e in G.edges(sort=True):
        a, b = e[:2]
        regs = [v for v in DG.vertices(sort=True) if a in v[1] and b in v[1]]
        if len(regs) == 2:
            DG.add_edge(regs[0], regs[1], e)
    return (G, E, p, EC, DG, vertical_regions)


def followstrand(f, factors, x0, x1, y0a, prec=53) -> list:
    r"""
    Return a piecewise linear approximation of the homotopy continuation
    of the root ``y0a`` from ``x0`` to ``x1``.

    INPUT:

    - ``f`` -- an irreducible polynomial in two variables
    - ``factors`` -- a list of irreducible polynomials in two variables
    - ``x0`` -- a complex value, where the homotopy starts
    - ``x1`` -- a complex value, where the homotopy ends
    - ``y0a`` -- an approximate solution of the polynomial `F(y) = f(x_0, y)`
    - ``prec`` -- the precision to use

    OUTPUT:

    A list of values `(t, y_{tr}, y_{ti})` such that:

    - ``t`` is a real number between zero and one
    - `f(t \cdot x_1 + (1-t) \cdot x_0, y_{tr} + I \cdot y_{ti})`
      is zero (or a good enough approximation)
    - the piecewise linear path determined by the points has a tubular
      neighborhood  where the actual homotopy continuation path lies, and
      no other root of ``f``, nor any root of the polynomials in ``factors``,
      intersects it.

    EXAMPLES::

        sage: # needs sirocco
        sage: from sage.schemes.curves.zariski_vankampen import followstrand
        sage: R.<x, y> = QQ[]
        sage: f = x^2 + y^3
        sage: x0 = CC(1, 0)
        sage: x1 = CC(1, 0.5)
        sage: followstrand(f, [], x0, x1, -1.0)             # abs tol 1e-15
        [(0.0, -1.0, 0.0),
         (0.7500000000000001, -1.015090921153253, -0.24752813818386948),
         (1.0, -1.026166099551513, -0.32768940253604323)]
        sage: fup = f.subs({y: y - 1/10})
        sage: fdown = f.subs({y: y + 1/10})
        sage: followstrand(f, [fup, fdown], x0, x1, -1.0)   # abs tol 1e-15
        [(0.0, -1.0, 0.0),
         (0.5303300858899107, -1.0076747107983448, -0.17588022709184917),
         (0.7651655429449553, -1.015686131039112, -0.25243563967299404),
         (1.0, -1.026166099551513, -0.3276894025360433)]
    """
    if f.degree() == 1:
        CF = ComplexField(prec)
        g = f.change_ring(CF)
        (x, y) = g.parent().gens()
        y0 = CF[y](g.subs({x: x0})).roots()[0][0]
        y1 = CF[y](g.subs({x: x1})).roots()[0][0]
        res = [(0.0, y0.real(), y0.imag()), (1.0, y1.real(), y1.imag())]
        return res
    CIF = ComplexIntervalField(prec)
    CC = ComplexField(prec)
    G = f.change_ring(QQbar).change_ring(CIF)
    x, y = G.parent().gens()
    g = G.subs({x: (1 - x) * CIF(x0) + x * CIF(x1)})
    coefs = []
    deg = g.total_degree()
    for d in range(deg + 1):
        for i in range(d + 1):
            c = CIF(g.coefficient({x: d - i, y: i}))
            cr = c.real()
            ci = c.imag()
            coefs += list(cr.endpoints())
            coefs += list(ci.endpoints())
    yr = CC(y0a).real()
    yi = CC(y0a).imag()
    coefsfactors = []
    degsfactors = []
    for fc in factors:
        degfc = fc.degree()
        degsfactors.append(degfc)
        G = fc.change_ring(QQbar).change_ring(CIF)
        g = G.subs({x: (1 - x) * CIF(x0) + x * CIF(x1)})
        for d in range(degfc + 1):
            for i in range(d + 1):
                c = CIF(g.coefficient({x: d - i, y: i}))
                cr = c.real()
                ci = c.imag()
                coefsfactors += list(cr.endpoints())
                coefsfactors += list(ci.endpoints())
    from sage.libs.sirocco import (contpath, contpath_mp, contpath_comps, contpath_mp_comps)
    try:
        if prec == 53:
            if factors:
                points = contpath_comps(deg, coefs, yr, yi, degsfactors, coefsfactors)
            else:
                points = contpath(deg, coefs, yr, yi)
        else:
            if factors:
                points = contpath_mp_comps(deg, coefs, yr, yi, prec, degsfactors, coefsfactors)
            else:
                points = contpath_mp(deg, coefs, yr, yi, prec)
        return points
    except Exception:
        return followstrand(f, factors, x0, x1, y0a, 2 * prec)


def newton(f, x0, i0):
    r"""
    Return the interval Newton operator.

    INPUT:

    - ``f`` -- a univariate polynomial
    - ``x0`` -- a number
    - ``I0`` -- an interval

    OUTPUT:

    The interval `x_0-\frac{f(x_0)}{f'(I_0)}`

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import newton
        sage: R.<x> = QQbar[]
        sage: f = x^3 + x
        sage: x0 = 1/10
        sage: I0 = RIF((-1/5,1/5))
        sage: n = newton(f, x0, I0)
        sage: n
        0.0?
        sage: n.real().endpoints()
        (-0.0147727272727274, 0.00982142857142862)
        sage: n.imag().endpoints()
        (0.000000000000000, -0.000000000000000)
    """
    return x0 - f(x0) / f.derivative()(i0)


def fieldI(field):
    r"""
    Return the (either double or trivial) extension of a number field which contains ``I``.

    INPUT:

    - ``field`` -- a number field with an embedding in `\QQbar`.

    OUTPUT:

    The extension ``F`` of ``field`` containing  ``I`` with  an embedding in `\QQbar`.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import fieldI
        sage: p = QQ[x](x^5 + 2 * x + 1)
        sage: a0 = p.roots(QQbar, multiplicities=False)[0]
        sage: F0.<a> = NumberField(p, embedding=a0)
        sage: fieldI(F0)
        Number Field in prim with defining polynomial
        x^10 + 5*x^8 + 14*x^6 - 2*x^5 - 10*x^4 + 20*x^3 - 11*x^2 - 14*x + 10
        with prim = 0.4863890359345430? + 1.000000000000000?*I
        sage: F0 = CyclotomicField(5)
        sage: fieldI(F0)
        Number Field in prim with defining polynomial
        x^8 - 2*x^7 + 7*x^6 - 10*x^5 + 16*x^4 - 10*x^3 - 2*x^2 + 4*x + 1
        with prim = -0.3090169943749474? + 0.04894348370484643?*I
        sage: fieldI(QuadraticField(3))
        Number Field in prim with defining polynomial x^4 - 4*x^2 + 16
        with prim = -1.732050807568878? + 1.000000000000000?*I
        sage: fieldI(QuadraticField(-3))
        Number Field in prim with defining polynomial x^4 + 8*x^2 + 4
        with prim = 0.?e-18 - 0.732050807568878?*I

    If ``I`` is already in the field, the result is the field itself::

        sage: from sage.schemes.curves.zariski_vankampen import fieldI
        sage: p = QQ[x](x^4 + 1)
        sage: a0 = p.roots(QQbar, multiplicities=False)[0]
        sage: F0.<a> = NumberField(p, embedding=a0)
        sage: F1 = fieldI(F0)
        sage: F0 == F1
        True
        sage: QuadraticField(-1) == fieldI(QuadraticField(-1))
        True
    """
    I0 = QQbar.gen()
    if I0 in field:
        return field
    field_a = field[I0]
    field_b = field_a.absolute_field('b0')
    b0 = field_b.gen()
    q = b0.minpoly()
    qembd = field_b.embeddings(QQbar)
    for h1 in qembd:
        b1 = h1(b0)
        b2 = h1(field_b(field_a.gen(0)))
        b3 = QQbar(field.gen(0))
        F1 = NumberField(q, 'prim', embedding=b1)
        if b3 in F1 and b2.imag() > 0:
            return F1


@parallel
def roots_interval(f, x0):
    """
    Find disjoint intervals that isolate the roots of a polynomial for a fixed
    value of the first variable.

    INPUT:

    - ``f`` -- a bivariate squarefree polynomial
    - ``x0`` -- a Gauss rational number corresponding to the first coordinate

    The intervals are taken as big as possible to be able to detect when two
    approximate roots of `f(x_0, y)` correspond to the same exact root, where
    `f` is the product of the polynomials in `flist`.

    The result is given as a dictionary, where the keys are
    approximations to the roots with rational real and imaginary
    parts, and the values are intervals containing them.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import roots_interval, fieldI
        sage: R.<x, y> = QQ[]
        sage: K = fieldI(QQ)
        sage: f = y^3 - x^2
        sage: f = f.change_ring(K)
        sage: ri = roots_interval(f, 1)
        sage: ri
        {-138907099/160396102*I - 1/2: -1.? - 1.?*I,
         138907099/160396102*I - 1/2: -1.? + 1.?*I,
         1: 1.? + 0.?*I}
        sage: [r.endpoints() for r in ri.values()]
        [(0.566987298107781 - 0.433012701892219*I,
          1.43301270189222 + 0.433012701892219*I,
          0.566987298107781 + 0.433012701892219*I,
          1.43301270189222 - 0.433012701892219*I),
         (-0.933012701892219 - 1.29903810567666*I,
          -0.0669872981077806 - 0.433012701892219*I,
          -0.933012701892219 - 0.433012701892219*I,
          -0.0669872981077806 - 1.29903810567666*I),
         (-0.933012701892219 + 0.433012701892219*I,
          -0.0669872981077806 + 1.29903810567666*I,
          -0.933012701892219 + 1.29903810567666*I,
          -0.0669872981077806 + 0.433012701892219*I)]
    """
    F1 = f.base_ring()
    x, y = f.parent().gens()
    fx = F1[y](f.subs({x: F1(x0)}))
    roots = fx.roots(QQbar, multiplicities=False)
    result = {}
    for i, r in enumerate(roots):
        prec = 53
        IF = ComplexIntervalField(prec)
        CF = ComplexField(prec)
        divisor = 4
        diam = min((CF(r) - CF(r0)).abs()
                   for r0 in roots[:i] + roots[i + 1:]) / divisor
        envelop = IF(diam) * IF((-1, 1), (-1, 1))
        while not newton(fx, r, r + envelop) in r + envelop:
            prec += 53
            IF = ComplexIntervalField(prec)
            CF = ComplexField(prec)
            divisor *= 2
            diam = min((CF(r) - CF(r0)).abs()
                       for r0 in roots[:i] + roots[i + 1:]) / divisor
            envelop = IF(diam) * IF((-1, 1), (-1, 1))
        qapr = QQ(CF(r).real()) + QQbar.gen() * QQ(CF(r).imag())
        if qapr not in r + envelop:
            raise ValueError("could not approximate roots with exact values")
        result[qapr] = r + envelop
    return result


def roots_interval_cached(f, x0):
    r"""
    Cached version of :func:`roots_interval`.

    TESTS::

        sage: from sage.schemes.curves.zariski_vankampen import roots_interval, roots_interval_cached, roots_interval_cache, fieldI
        sage: R.<x, y> = QQ[]
        sage: K = fieldI(QQ)
        sage: f = y^3 - x^2
        sage: f = f.change_ring(K)
        sage: (f, 1) in roots_interval_cache
        False
        sage: ri = roots_interval_cached(f, 1)
        sage: ri
        {-138907099/160396102*I - 1/2: -1.? - 1.?*I,
         138907099/160396102*I - 1/2: -1.? + 1.?*I,
         1: 1.? + 0.?*I}
        sage: (f, 1) in roots_interval_cache
        True
    """
    global roots_interval_cache
    try:
        return roots_interval_cache[(f, x0)]
    except KeyError:
        result = roots_interval(f, x0)
        roots_interval_cache[(f, x0)] = result
        return result


def populate_roots_interval_cache(inputs):
    r"""
    Call :func:`roots_interval` to the inputs that have not been
    computed previously, and cache them.

    INPUT:

    - ``inputs`` -- a list of tuples ``(f, x0)``

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import populate_roots_interval_cache, roots_interval_cache, fieldI
        sage: R.<x,y> = QQ[]
        sage: K=fieldI(QQ)
        sage: f = y^5 - x^2
        sage: f = f.change_ring(K)
        sage: (f, 3) in roots_interval_cache
        False
        sage: populate_roots_interval_cache([(f, 3)])
        sage: (f, 3) in roots_interval_cache
        True
        sage: roots_interval_cache[(f, 3)]
        {-1.255469441943070? - 0.9121519421827974?*I: -2.? - 1.?*I,
         -1.255469441943070? + 0.9121519421827974?*I: -2.? + 1.?*I,
         0.4795466549853897? - 1.475892845355996?*I: 1.? - 2.?*I,
         0.4795466549853897? + 1.475892845355996?*I: 1.? + 2.?*I,
         14421467174121563/9293107134194871: 2.? + 0.?*I}

    """
    global roots_interval_cache
    tocompute = [inp for inp in inputs if inp not in roots_interval_cache]
    problem_par = True
    while problem_par:  # hack to deal with random fails in parallelization
        try:
            result = roots_interval(tocompute)
            for r in result:
                roots_interval_cache[r[0][0]] = r[1]
            problem_par = False
        except TypeError:
            pass


@parallel
def braid_in_segment(glist, x0, x1, precision={}):
    """
    Return the braid formed by the `y` roots of ``f`` when `x` moves
    from ``x0`` to ``x1``.

    INPUT:

    - ``glist`` -- a tuple of polynomials in two variables
    - ``x0`` -- a Gauss rational
    - ``x1`` -- a Gauss rational
    - ``precision`` -- a dictionary (default: `{}`) which assigns a number
      precision bits to each element of ``glist``

    OUTPUT:

    A braid.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import braid_in_segment, fieldI
        sage: R.<x, y> = QQ[]
        sage: K = fieldI(QQ)
        sage: f = x^2 + y^3
        sage: f = f.change_ring(K)
        sage: x0 = 1
        sage: x1 = 1 + I / 2
        sage: braid_in_segment(tuple(_[0] for _ in f.factor()), x0, x1)     # needs sirocco
        s1

    TESTS:

    Check that :issue:`26503` is fixed::

        sage: # needs sage.rings.real_mpfr sage.symbolic
        sage: wp = QQ['t']([1, 1, 1]).roots(QQbar)[0][0]
        sage: Kw.<wp> = NumberField(wp.minpoly(), embedding=wp)
        sage: R.<x, y> = Kw[]
        sage: z = -wp - 1
        sage: f = y * (y + z) * x * (x - 1) * (x - y) * (x + z * y - 1) * (x + z * y + wp)
        sage: from sage.schemes.curves.zariski_vankampen import fieldI, braid_in_segment
        sage: Kw1 = fieldI(Kw)
        sage: g = f.subs({x: x + 2 * y})
        sage: g = g.change_ring(Kw1)
        sage: p1 = QQbar(sqrt(-1/3))
        sage: p1a = CC(p1)
        sage: p1b = QQ(p1a.real()) + I*QQ(p1a.imag())
        sage: p2 = QQbar(1/2 + sqrt(-1/3)/2)
        sage: p2a = CC(p2)
        sage: p2b = QQ(p2a.real()) + I*QQ(p2a.imag())
        sage: glist = tuple([_[0] for _ in g.factor()])
        sage: B = braid_in_segment(glist, p1b, p2b); B              # needs sirocco
        s5*s3^-1
    """
    precision1 = precision.copy()
    g = prod(glist)
    F1 = g.base_ring()
    x, y = g.parent().gens()
    intervals = {}
    if not precision1:
        precision1 = {f: 53 for f in glist}
    y0s = []
    for f in glist:
        if f.variables() == (y,):
            f0 = F1[y](f)
        else:
            f0 = F1[y](f.subs({x: F1(x0)}))
        y0sf = f0.roots(QQbar, multiplicities=False)
        y0s += list(y0sf)
        while True:
            CIFp = ComplexIntervalField(precision1[f])
            intervals[f] = [r.interval(CIFp) for r in y0sf]
            if not any(a.overlaps(b) for a, b in
                       itertools.combinations(intervals[f], 2)):
                break
            precision1[f] *= 2
    strands = []
    for f in glist:
        for i in intervals[f]:
            aux = followstrand(f, [p for p in glist if p != f],
                               x0, x1, i.center(), precision1[f])
            strands.append(aux)
    complexstrands = [[(QQ(a[0]), QQ(a[1]), QQ(a[2])) for a in b]
                      for b in strands]
    centralbraid = braid_from_piecewise(complexstrands)
    initialstrands = []
    finalstrands = []
    initialintervals = roots_interval_cached(g, x0)
    finalintervals = roots_interval_cached(g, x1)
    I1 = QQbar.gen()
    for cs in complexstrands:
        ip = cs[0][1] + I1 * cs[0][2]
        fp = cs[-1][1] + I1 * cs[-1][2]
        matched = 0
        for center, interval in initialintervals.items():
            if ip in interval:
                initialstrands.append([(0, center.real(), center.imag()),
                                       (1, cs[0][1], cs[0][2])])
                matched += 1
        if matched != 1:
            precision1 = {f: precision1[f] * 2 for f in glist}
            return braid_in_segment(glist, x0, x1, precision=precision1)

        matched = 0
        for center, interval in finalintervals.items():
            if fp in interval:
                finalstrands.append([(0, cs[-1][1], cs[-1][2]),
                                     (1, center.real(), center.imag())])
                matched += 1
        if matched != 1:
            precision1 = {f: precision1[f] * 2 for f in glist}
            return braid_in_segment(glist, x0, x1, precision=precision1)
    initialbraid = braid_from_piecewise(initialstrands)
    finalbraid = braid_from_piecewise(finalstrands)

    return initialbraid * centralbraid * finalbraid


def geometric_basis(G, E, EC0, p, dual_graph, vertical_regions={}) -> list:
    r"""
    Return a geometric basis, based on a vertex.

    INPUT:

    - ``G`` -- a graph with the bounded edges of a Voronoi Diagram

    - ``E`` -- a subgraph of ``G`` which is a cycle containing the bounded
      edges touching an unbounded region of a Voronoi Diagram

    - ``EC0`` -- A counterclockwise orientation of the vertices of ``E``

    - ``p`` -- a vertex of ``E``

    - ``dual_graph`` -- a dual graph for a plane embedding of ``G`` such that
      ``E`` is the boundary of the non-bounded component of the complement.
      The edges are labelled as the dual edges and the vertices are labelled
      by a tuple whose first element is the an integer for the position and the
      second one is the cyclic ordered list of vertices in the region

    - ``vertical_regions`` -- dictionary (default: `{}`); its keys are
      the vertices of ``dual_graph`` to fix regions associated with
      vertical lines

    OUTPUT: A geometric basis and a dictionary.

    The geometric basis is formed by a list of sequences of paths. Each path is a
    ist of vertices, that form a closed path in ``G``, based at ``p``, that goes
    to a region, surrounds it, and comes back by the same path it came. The
    concatenation of all these paths is equivalent to ``E``.

    The dictionary associates to each vertical line the index of the generator
    of the geometric basis associated to it.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import geometric_basis, corrected_voronoi_diagram, voronoi_cells
        sage: points = (0, -1, I, 1, -I)
        sage: V = corrected_voronoi_diagram(points)
        sage: G, E, p, EC, DG, VR = voronoi_cells(V, vertical_lines=frozenset((0 .. 4)))
        sage: gb, vd = geometric_basis(G, E, EC, p, DG, vertical_regions=VR)
        sage: gb
        [[A vertex at (5/2, -5/2), A vertex at (5/2, 5/2), A vertex at (-5/2, 5/2),
          A vertex at (-1/2, 1/2), A vertex at (-1/2, -1/2), A vertex at (1/2, -1/2),
          A vertex at (1/2, 1/2), A vertex at (-1/2, 1/2), A vertex at (-5/2, 5/2),
          A vertex at (5/2, 5/2), A vertex at (5/2, -5/2)],
         [A vertex at (5/2, -5/2), A vertex at (5/2, 5/2), A vertex at (-5/2, 5/2),
          A vertex at (-1/2, 1/2), A vertex at (1/2, 1/2), A vertex at (5/2, 5/2),
          A vertex at (5/2, -5/2)],
         [A vertex at (5/2, -5/2), A vertex at (5/2, 5/2), A vertex at (1/2, 1/2),
          A vertex at (1/2, -1/2), A vertex at (5/2, -5/2)], [A vertex at (5/2, -5/2),
          A vertex at (1/2, -1/2), A vertex at (-1/2, -1/2), A vertex at (-1/2, 1/2),
          A vertex at (-5/2, 5/2), A vertex at (-5/2, -5/2), A vertex at (-1/2, -1/2),
          A vertex at (1/2, -1/2), A vertex at (5/2, -5/2)],
         [A vertex at (5/2, -5/2), A vertex at (1/2, -1/2), A vertex at (-1/2, -1/2),
          A vertex at (-5/2, -5/2), A vertex at (5/2, -5/2)]]
        sage: vd
        {0: 0, 1: 3, 2: 1, 3: 2, 4: 4}
    """
    i = EC0.index(p)
    EC = EC0[i:-1] + EC0[:i + 1]
    # A counterclockwise eulerian circuit on the boundary,
    # starting and ending at p
    if G.size() == E.size():
        if E.is_cycle():
            j = next(dual_graph.vertex_iterator())[0]
            if j in vertical_regions:
                vd = {j: 0}
            else:
                vd = {}
            return [EC], vd
    edges_E = E.edges(sort=True)
    InternalEdges = [e for e in G.edges(sort=True) if e not in edges_E]
    InternalVertices = [v for e in InternalEdges for v in e[:2]]
    Internal = G.subgraph(vertices=InternalVertices, edges=InternalEdges)
    for i, ECi in enumerate(EC):  # q and r are the points we will cut through
        if ECi in Internal:
            EI = [v for v in E if v in
                  Internal.connected_component_containing_vertex(ECi, sort=True)
                  and v != ECi]
            if EI:
                q = ECi
                connecting_path = list(EC[:i])
                break
        if EC[-i] in Internal:
            EI = [v for v in E if v in
                  Internal.connected_component_containing_vertex(EC[-i], sort=True)
                  and v != EC[-i]]
            if EI:
                q = EC[-i]
                connecting_path = list(reversed(EC[-i:]))
                break
    # Precompute distances from q in E and I
    E_dist_q = E.shortest_path_lengths(q)
    I_dist_q = Internal.shortest_path_lengths(q)
    distancequotients = [(E_dist_q[v]**2 / I_dist_q[v], v) for v in EI]
    r = max(distancequotients)[1]
    cutpath = Internal.shortest_path(q, r)
    for i, v in enumerate(cutpath):
        if i > 0 and v in EC:
            r = v
            cutpath = cutpath[:i+1]
            break
    qi = EC.index(q)
    ri = EC.index(r)
    Ecut = copy(E)
    Ecut.delete_vertices([q, r])
    subgraphs = Ecut.connected_components_subgraphs()
    if len(subgraphs) == 2:
        E1, E2 = subgraphs
        if EC[qi + 1] in E2:
            E1, E2 = E2, E1
    elif len(subgraphs) == 1:
        E1 = subgraphs[0]
        E2 = Graph()
        E2.add_edge(q, r, None)
        if r == EC[qi + 1]:
            E1, E2 = E2, E1
    for v in [q, r]:
        for n in E.neighbor_iterator(v):
            if n in E1 and n not in (q, r):
                E1.add_edge(v, n, None)
            if n in E2 and n not in (q, r):
                E2.add_edge(v, n, None)

    for i in range(len(cutpath) - 1):
        E1.add_edge(cutpath[i], cutpath[i + 1], None)
        E2.add_edge(cutpath[i], cutpath[i + 1], None)
    Gd = copy(dual_graph)
    to_delete = [e for e in Gd.edges(sort=True) if e[2][0] in cutpath and
                 e[2][1] in cutpath]
    Gd.delete_edges(to_delete)
    Gd1, Gd2 = Gd.connected_components_subgraphs()
    edges_2 = []
    vertices_2 = []
    for reg in Gd2.vertices(sort=True):
        vertices_2 += reg[1][:-1]
        reg_circuit = reg[1]
        edges_2 += [(v1, reg_circuit[i + 1])
                    for i, v1 in enumerate(reg_circuit[:-1])]
        edges_2 += [(v1, reg_circuit[i - 1])
                    for i, v1 in enumerate(reg_circuit[1:])]
    G2 = G.subgraph(vertices=vertices_2, edges=edges_2)
    edges_1 = []
    vertices_1 = []
    for reg in Gd1.vertices(sort=True):
        vertices_1 += reg[1]
        reg_circuit = reg[1] + (reg[1][0],)
        edges_1 += [(v1, reg_circuit[i + 1])
                    for i, v1 in enumerate(reg_circuit[:-1])]
        edges_1 += [(v1, reg_circuit[i - 1])
                    for i, v1 in enumerate(reg_circuit[1:])]
    G1 = G.subgraph(vertices=vertices_1, edges=edges_1)
    if EC[qi + 1] in G2:
        G1, G2 = G2, G1
        Gd1, Gd2 = Gd2, Gd1

    if qi < ri:
        EC1 = [EC[j] for j in range(qi, ri)] + list(reversed(cutpath))
        EC2 = cutpath + list(EC[ri + 1: -1] + EC[: qi + 1])
    else:
        EC1 = list(EC[qi:] + EC[1:ri]) + list(reversed(cutpath))
        EC2 = cutpath + list(EC[ri + 1:qi + 1])

    gb1, vd1 = geometric_basis(G1, E1, EC1, q, Gd1, vertical_regions=vertical_regions)
    gb2, vd2 = geometric_basis(G2, E2, EC2, q, Gd2, vertical_regions=vertical_regions)

    vd = {j: vd1[j] for j in vd1}
    m = len(gb1)
    for j in vd2.keys():
        vd[j] = vd2[j] + m
    reverse_connecting = list(reversed(connecting_path))
    resul = [connecting_path + path + reverse_connecting
             for path in gb1 + gb2]
    for r in resul:
        i = 0
        while i < len(r) - 2:
            if r[i] == r[i + 2]:
                r.pop(i)
                r.pop(i)
                if i:
                    i -= 1
            else:
                i += 1
    return (resul, vd)


def vertical_lines_in_braidmon(pols) -> list:
    r"""
    Return the vertical lines in ``pols``, unless
    one of the other components has a vertical asymptote.

    INPUT:

    - ``pols`` -- a  list of polynomials with two variables whose
      product equals ``f``

    OUTPUT:

    A list with the indices of the vertical lines in ``flist`` if there is
    no other componnet with vertical asymptote; otherwise it returns an empty
    list.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import vertical_lines_in_braidmon
        sage: R.<x, y> = QQ[]
        sage: flist = [x^2 - y^3, x, x + 3 * y - 5, 1 - x]
        sage: vertical_lines_in_braidmon(flist)
        [1, 3]
        sage: flist += [x * y - 1]
        sage: vertical_lines_in_braidmon(flist)
        []
        sage: vertical_lines_in_braidmon([])
        []
    """
    if not pols:
        return []
    res = []
    for j, f in enumerate(pols):
        C = Curve(f)
        vertical_asymptote = C.has_vertical_asymptote()
        if vertical_asymptote:
            return []
        if C.is_vertical_line():
            res.append(j)
    return res


def strand_components(f, pols, p1):
    r"""
    Compute only the assignment from strands to elements of ``flist``.

    INPUT:

    - ``f`` -- a  reduced polynomial with two variables, over a number field
      with an embedding in the complex numbers

    - ``pols`` -- a  list of polynomials with two variables whose
      product equals ``f``

    - ``p1`` -- a Gauss rational

    OUTPUT:

    - A list and a dictionary.  The first one is an ordered list of pairs
      consisting of ``(z,i)`` where ``z`` is a root of ``f(p_1,y)``
      and `i` is the position of the polynomial in the list whose root
      is ``z``. The second one attaches a number `i` (strand) to a
      number `j` (a polynomial in the list).

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import strand_components
        sage: R.<x, y> = QQ[]
        sage: flist = [x^2 - y^3, x + 3 * y - 5]
        sage: strand_components(prod(flist), flist, 1)
        ([(-0.500000000000000? - 0.866025403784439?*I, 0),
          (-0.500000000000000? + 0.866025403784439?*I, 0),
          (1, 0), (1.333333333333334?, 1)], {0: 0, 1: 0, 2: 0, 3: 1})
    """
    x, y = f.parent().gens()
    F = pols[0].base_ring()
    strands = {}
    roots_base = []
    for i, h in enumerate(pols):
        h0 = h.subs({x: p1})
        h1 = F[y](h0)
        rt = h1.roots(QQbar, multiplicities=False)
        roots_base += [(r, i) for r in rt]
    roots_base.sort()
    strands = {i: par[1] for i, par in enumerate(roots_base)}
    return (roots_base, strands)


def braid_monodromy(f, arrangement=(), vertical=False):
    r"""
    Compute the braid monodromy of a projection of the curve defined by
    a polynomial.

    INPUT:

    - ``f`` -- a polynomial with two variables, over a number field
      with an embedding in the complex numbers

    - ``arrangement`` -- tuple (default: ``()``); an optional tuple
      of polynomials whose product equals ``f``

    - ``vertical`` -- boolean (default: ``False`); if set to ``True``,
      ``arrangements`` contains more than one polynomial, some of them
      are of degree `1` in `x` and degree `0` in `y`, and none of
      the other components have vertical asymptotes, then these
      components are marked as *vertical* and not used for the computation
      of the braid monodromy. The other ones are marked as *horizontal*. If
      a vertical component does not pass through a singular points of the
      projection of the horizontal components a trivial braid is added
      to the list.

    OUTPUT:

    - A list of braids, images by the braid monodromy of a geometric
      basis of the complement of the discriminant of `f` in `\CC`.

    - A dictionary: ``i``, index of a strand is sent to the index of
      the corresponding factor in ``arrangement``.

    - Another dictionary ``dv``, only relevant if ``vertical`` is ``True``.
      If  ``j`` is the index
      of a braid corresponding to a vertical line with index ``i``
      in ``arrangement``, then ``dv[j] = i``.

    - A non-negative integer: the number of strands of the braids,
      only necessary if the list of braids is empty.

    .. NOTE::

        The projection over the `x` axis is used if there are no vertical
        asymptotes. Otherwise, a linear change of variables is done to fall
        into the previous case except if the only vertical asymptotes are lines
        and ``vertical=True``.

    EXAMPLES::

        sage: # needs sirocco
        sage: from sage.schemes.curves.zariski_vankampen import braid_monodromy
        sage: R.<x, y> = QQ[]
        sage: f = (x^2 - y^3) * (x + 3*y - 5)
        sage: bm = braid_monodromy(f); bm
        ([s1*s0*(s1*s2)^2*s0*s2^2*s0^-1*(s2^-1*s1^-1)^2*s0^-1*s1^-1,
          s1*s0*(s1*s2)^2*(s0*s2^-1*s1*s2*s1*s2^-1)^2*(s2^-1*s1^-1)^2*s0^-1*s1^-1,
          s1*s0*(s1*s2)^2*s2*s1^-1*s2^-1*s1^-1*s0^-1*s1^-1,
          s1*s0*s2*s0^-1*s2*s1^-1], {0: 0, 1: 0, 2: 0, 3: 0}, {}, 4)
        sage: flist = (x^2 - y^3, x + 3*y - 5)
        sage: bm1 = braid_monodromy(f, arrangement=flist)
        sage: bm1[0] == bm[0]
        True
        sage: bm1[1]
        {0: 0, 1: 1, 2: 0, 3: 0}
        sage: braid_monodromy(R(1))
        ([], {}, {}, 0)
        sage: braid_monodromy(x*y^2 - 1)
        ([s0*s1*s0^-1*s1*s0*s1^-1*s0^-1, s0*s1*s0^-1, s0], {0: 0, 1: 0, 2: 0}, {}, 3)
        sage: L = [x, y, x - 1, x -y]
        sage: braid_monodromy(prod(L), arrangement=L, vertical=True)
        ([s^2, 1], {0: 1, 1: 3}, {0: 0, 1: 2}, 2)
    """
    global roots_interval_cache
    F = fieldI(f.base_ring())
    I1 = F(QQbar.gen())
    f = f.change_ring(F)
    if not arrangement:
        arrangement1 = (f,)
    else:
        arrangement1 = tuple(g.change_ring(F) for g in arrangement)
    x, y = f.parent().gens()
    if vertical:
        indices_v = vertical_lines_in_braidmon(arrangement1)
    else:
        indices_v = []
    arrangement_h = tuple(f0 for j, f0 in enumerate(arrangement1)
                          if j not in indices_v)
    arrangement_v = tuple(f0 for j, f0 in enumerate(arrangement1)
                          if j in indices_v)
    glist = tuple(fc[0] for f0 in arrangement_h for fc in f0.factor())
    g = f.parent()(prod(glist))
    d = g.degree(y)
    if not arrangement_v:  # change of coordinates only if indices_v is empty
        while not g.coefficient(y**d) in F:
            g = g.subs({x: x + y})
            d = g.degree(y)
            arrangement_h = tuple(f1.subs({x: x + y}) for f1 in arrangement_h)
            arrangement1 = arrangement_h
            glist = tuple(f1.subs({x: x + y}) for f1 in glist)
    if d > 0:
        disc = discrim(glist)
    else:
        disc = []
    vertical_braid = {}
    transversal = {}
    vl = []
    for f0 in arrangement_v:
        pt = [j for j, t in enumerate(disc) if f0.subs({x: t}) == 0]
        if pt:
            vertical_braid[f0] = (pt[0], arrangement1.index(f0))
            vl.append(pt[0])
        else:
            transversal[f0] = arrangement1.index(f0)
    vl.sort()
    vl = frozenset(vl)
    if not disc:
        vertical_braids = {i: transversal[f0]
                           for i, f0 in enumerate(transversal)}
        if d > 1:
            result = [BraidGroup(d).one() for p in transversal]
        else:
            G = FreeGroup(0) / []
            result = [G.one() for p in transversal]
        p1 = F(0)
        if d > 0:
            roots_base, strands = strand_components(g, arrangement_h, p1)
            strands1 = {}
            for j in range(d):
                i = strands[j]
                k = arrangement1.index(arrangement_h[i])
                strands1[j] = k
        else:
            strands1 = {}
        return (result, strands1, vertical_braids, d)
    V = corrected_voronoi_diagram(tuple(disc))
    G, E, p, EC, DG, VR = voronoi_cells(V, vertical_lines=vl)
    p0 = (p[0], p[1])
    p1 = p0[0] + I1 * p0[1]
    roots_base, strands = strand_components(g, arrangement_h, p1)
    strands1 = {}
    for j in range(d):
        i = strands[j]
        k = arrangement1.index(arrangement_h[i])
        strands1[j] = k
    geombasis, vd = geometric_basis(G, E, EC, p, DG, vertical_regions=VR)
    segs = set()
    for p in geombasis:
        for s in zip(p[:-1], p[1:]):
            if (s[1], s[0]) not in segs:
                segs.add((s[0], s[1]))
    I0 = QQbar.gen()
    segs = [(a[0] + I0 * a[1], b[0] + I0 * b[1]) for a, b in segs]
    vertices = list(set(flatten(segs)))
    tocacheverts = tuple([(g, v) for v in vertices])
    populate_roots_interval_cache(tocacheverts)
    end_braid_computation = False
    while not end_braid_computation:
        try:
            braidscomputed = (braid_in_segment([(glist, seg[0], seg[1])
                                                for seg in segs]))
            segsbraids = {}
            for braidcomputed in braidscomputed:
                seg = (braidcomputed[0][0][1], braidcomputed[0][0][2])
                beginseg = (QQ(seg[0].real()), QQ(seg[0].imag()))
                endseg = (QQ(seg[1].real()), QQ(seg[1].imag()))
                b = braidcomputed[1]
                segsbraids[(beginseg, endseg)] = b
                segsbraids[(endseg, beginseg)] = b.inverse()
            end_braid_computation = True
        except ChildProcessError:  # hack to deal with random fails first time
            pass
    B = BraidGroup(d)
    result = []
    for path in geombasis:
        braidpath = B.one()
        for i in range(len(path) - 1):
            x0 = tuple(path[i].vector())
            x1 = tuple(path[i + 1].vector())
            braidpath = braidpath * segsbraids[(x0, x1)]
        result.append(braidpath)
    vertical_braids = {}
    r = len(result)
    t = 0
    for f0 in arrangement_v:
        if f0 in vertical_braid.keys():
            k, j = vertical_braid[f0]
            vertical_braids[vd[k]] = j
        else:
            vertical_braids[r + t] = transversal[f0]
            t += 1
            result.append(B.one())
    return (result, strands1, vertical_braids, d)


def conjugate_positive_form(braid):
    r"""
    For a ``braid`` which is conjugate to a product of *disjoint* positive
    braids a list of such decompositions is given.

    INPUT:

    - ``braid`` -- a braid `\sigma`

    OUTPUT:

    A list of `r` lists. Each such list is another list with two elements, a
    positive braid `\alpha_i` and a list of permutation braids
    `\gamma_{1}^{i},\dots,\gamma_{r}^{n_i}` such that if
    `\gamma_i=\prod_{j=1}^{n_i} \gamma_j^i` then the braids
    `\tau_i=\gamma_i\alpha_i\gamma_i^{-1}` pairwise commute
    and `\alpha=\prod_{i=1}^{r} \tau_i`.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import conjugate_positive_form
        sage: B = BraidGroup(4)
        sage: t = B((1, 3, 2, -3, 1, 1))
        sage: conjugate_positive_form(t)
        [[(s1*s0)^2, [s2]]]
        sage: B = BraidGroup(5)
        sage: t = B((1, 2, 3, 4, -1, -2, 3, 3, 2, -4))
        sage: L = conjugate_positive_form(t); L
        [[s1^2, [s3*s2]], [s1*s2, [s0]]]
        sage: s = B.one()
        sage: for a, l in L:
        ....:   b = prod(l)
        ....:   s *= b * a / b
        sage: s == t
        True
        sage: s1 = B.gen(1)^3
        sage: conjugate_positive_form(s1)
        [[s1^3, []]]
    """
    B = braid.parent()
    d = B.strands()
    braid1 = braid.super_summit_set()[0]
    L1 = braid1.Tietze()
    sg0 = braid.conjugating_braid(braid1)
    gns = set(L1)
    cuts = [j for j in range(d + 1) if j not in gns]
    blocks = []
    for i in range(len(cuts) - 1):
        block = [j for j in L1 if cuts[i] < j < cuts[i + 1]]
        if block:
            blocks.append(block)
    shorts = []
    for a in blocks:
        A = B(a).super_summit_set()
        res = None
        for tau in A:
            sg = (sg0 * B(a) / sg0).conjugating_braid(tau)
            A1 = rightnormalform(sg)
            par = A1[-1][0] % 2
            A1 = [B(a) for a in A1[:-1]]
            b = prod(A1, B.one())
            b1 = len(b.Tietze()) / (len(A1) + 1)
            if res is None or b1 < res[3]:
                res = [tau, A1, par, b1]
        if res[2] == 1:
            r0 = res[0].Tietze()
            res[0] = B([i.sign() * (d - abs(i)) for i in r0])
        res0 = res[:2]
        shorts.append(res0)
    return shorts


@parallel
def conjugate_positive_form_p(braid):
    return conjugate_positive_form(braid)


def braid2rels(L):
    r"""
    Return a minimal set of relations of the group
    ``F / [(b * F([j])) / F([j]) for j in (1..d)]`` where ``F = FreeGroup(d)``
    and ``b`` is a conjugate of a positive braid . One starts from the
    non-trivial relations determined by the positive braid and transform
    them in relations determined by ``b``.

    INPUT:

    - ``L`` -- a tuple whose first element is a positive braid and the second
      element is a list of permutation braids.

    OUTPUT:

    A list of Tietze words for a minimal set of relations of
    ``F / [(g * b) / g for g in F.gens()]``.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import braid2rels
        sage: B.<s0, s1, s2> = BraidGroup(4)
        sage: L = ((s1*s0)^2, [s2])
        sage: braid2rels(L)
        [(4, 1, -2, -1), (2, -4, -2, 1)]
    """
    br = L[0]
    L1 = L[1]
    B = br.parent()
    d = B.strands()
    F = FreeGroup(d)
    T = br.Tietze()
    T1 = set(T)
    m = len(T1) + 1
    k = min(T1) - 1
    B0 = BraidGroup(m)
    F0 = FreeGroup(m)
    br0 = B0([j - k for j in T])
    br0_left = leftnormalform(br0)
    q, r = ZZ(br0_left[0][0]).quo_rem(2)
    br1 = B0.delta()**r * prod(map(B0, br0_left[1:]), B0.one())
    cox = prod(F0.gens())
    U0 = [cox**q * (f0 * br1) / cox**q / f0 for f0 in F0.gens()[:-1]]
    U = [tuple(sign(k1) * (abs(k1) + k) for k1 in br.Tietze()) for br in U0]
    pasos = [B.one()]
    pasos.extend(reversed(L1))
    for C in pasos:
        U = [(F(a) * C.inverse()).Tietze() for a in U]
        ga = F / U
        P = ga.gap().PresentationFpGroup()
        dic = P.TzOptions().sage()
        dic['protected'] = d
        dic['printLevel'] = 0
        P.SetTzOptions(dic)
        P.TzGoGo()
        P.TzGoGo()
        gb = wrap_FpGroup(P.FpGroupPresentation())
        U = [rel.Tietze() for rel in gb.relations()]
    return U


@parallel
def braid2rels_p(L):
    return braid2rels(L)


@parallel
def relation(x, b):
    return x * b / x


def fundamental_group_from_braid_mon(bm, degree=None,
                                     simplified=True, projective=False,
                                     puiseux=True, vertical=[]):
    r"""
    Return a presentation of the fundamental group computed from
    a braid monodromy.

    INPUT:

    - ``bm`` -- a list of braids

    - ``degree`` -- integer (default: ``None``); only needed if the braid
      monodromy is an empty list.

    - ``simplified`` -- boolean (default: ``True``); if set to ``True`` the
      presentation will be simplified (see below)

    - ``projective`` -- boolean (default: ``False``); if set to ``True``,
      the fundamental group of the complement of the projective completion
      of the curve will be computed, otherwise, the fundamental group of
      the complement in the affine plane will be computed

    - ``puiseux`` -- boolean (default: ``True``); if set to ``True``
      a presentation of the fundamental group with the homotopy type
      of the complement of the affine curve will be computed, adding
      one relation if ``projective`` is set to ``True``.

    - ``vertical`` -- list of integers (default: ``[]``); the indices in
      ``[0 .. r - 1]`` of the braids that surround a vertical line

    If ``projective``` is ``False`` and ``puiseux`` is
    ``True``, a Zariski-VanKampen presentation is returned.

    OUTPUT:

    A presentation of the fundamental group of the complement of the
    union of the curve with some vertical lines from its braid monodromy.

    EXAMPLES::

        sage: from sage.schemes.curves.zariski_vankampen import fundamental_group_from_braid_mon
        sage: B.<s0, s1, s2> = BraidGroup(4)
        sage: bm = [s1*s2*s0*s1*s0^-1*s1^-1*s0^-1,
        ....:       s0*s1^2*s0*s2*s1*(s0^-1*s1^-1)^2*s0^-1,
        ....:       (s0*s1)^2]
        sage: g = fundamental_group_from_braid_mon(bm, projective=True); g      # needs sirocco
        Finitely presented group
        < x1, x3 | x3^2*x1^2, x1^-1*x3^-1*x1*x3^-1*x1^-1*x3^-1 >
        sage: print(g.order(), g.abelian_invariants())                         # needs sirocco
        12 (4,)
        sage: B2 = BraidGroup(2)
        sage: bm = [B2(3 * [1])]
        sage: g = fundamental_group_from_braid_mon(bm, vertical=[0]); g         # needs sirocco
        Finitely presented group
        < x0, x1, x2 | x2*x0*x1*x2^-1*x1^-1*x0^-1,
                       x2*x0*x1*x0*x1^-1*x0^-1*x2^-1*x1^-1 >
        sage: fundamental_group_from_braid_mon([]) is None                      # needs sirocco
        True
        sage: fundamental_group_from_braid_mon([], degree=2)                    # needs sirocco
        Finitely presented group < x0, x1 |  >
        sage: fundamental_group_from_braid_mon([SymmetricGroup(1).one()])       # needs sirocco
        Finitely presented group < x |  >
    """
    vertical0 = sorted(vertical)
    v = len(vertical0)
    if not bm:
        d = degree
    elif bm[0].parent().order() == 1:
        d = 1
    else:
        d = bm[0].parent().strands()
    if d is None:
        return None
    F = FreeGroup(d)
    Fv = FreeGroup(d + v)
    if d == 0:
        return Fv / []
    if d == 1:
        return Fv / [(1, j, -1, -j) for j in range(2, d + v + 1)]
    bmh = [br for j, br in enumerate(bm) if j not in vertical0]
    if not puiseux:
        relations_h = (relation([(x, b) for x in F.gens() for b in bmh]))
        rel_h = [r[1] for r in relations_h]
    else:
        conjugate_desc = conjugate_positive_form_p(bmh)
        trenzas_desc = [b1[-1] for b1 in conjugate_desc]
        trenzas_desc_1 = flatten(trenzas_desc, max_level=1)
        relations_h = braid2rels_p(trenzas_desc_1)
        rel_h = [r[1] for r in relations_h]
        rel_h = flatten(rel_h, max_level=1)
    rel_v = []
    for j, k in enumerate(vertical0):
        l1 = d + j + 1
        br = bm[k]
        for gen in F.gens():
            j0 = gen.Tietze()[0]
            rl = (l1,) + (gen * br).Tietze() + (-l1, -j0)
            rel_v.append(rl)
    rel = rel_h + rel_v
    if projective:
        rel.append(prod(Fv.gens()).Tietze())
    G = Fv / rel
    if simplified:
        return G.simplified()
    return G


def fundamental_group(f, simplified=True, projective=False, puiseux=True):
    r"""
    Return a presentation of the fundamental group of the complement of
    the algebraic set defined by the polynomial ``f``.

    INPUT:

    - ``f`` -- a polynomial in two variables, with coefficients in either
      the rationals or a number field with a fixed embedding in `\QQbar`

    - ``simplified`` -- boolean (default: ``True``); if set to ``True`` the
      presentation will be simplified (see below)

    - ``projective`` -- boolean (default: ``False``); if set to ``True``,
      the fundamental group of the complement of the projective completion
      of the curve will be computed, otherwise, the fundamental group of
      the complement in the affine plane will be computed

    - ``puiseux`` -- boolean (default: ``True``); if set to ``True``,
      a presentation of the fundamental group with the homotopy type
      of the complement of the affine curve is computed. If the Euler
      characteristic does not match, the homotopy type is obtained
      with a wedge of 2-spheres. One relation is added if ``projective``
      is set to ``True``.

    If ``projective``` is ``False`` and ``puiseux`` is
    ``True``, a Zariski-VanKampen presentation is returned.

    OUTPUT:

    A presentation of the fundamental group of the complement of the
    curve defined by ``f``.

    EXAMPLES::

        sage: # needs sirocco
        sage: from sage.schemes.curves.zariski_vankampen import fundamental_group, braid_monodromy
        sage: R.<x, y> = QQ[]
        sage: f = x^2 + y^3
        sage: fundamental_group(f)
        Finitely presented group < x0, x1 | x0*x1^-1*x0^-1*x1^-1*x0*x1 >
        sage: fundamental_group(f, simplified=False, puiseux=False).sorted_presentation()
        Finitely presented group < x0, x1, x2 | x2^-1*x1^-1*x0*x1,
                                                x2^-1*x0*x1*x0^-1,
                                                x1^-1*x0^-1*x1^-1*x0*x1*x0 >
        sage: fundamental_group(f, projective=True)
        Finitely presented group < x0 | x0^3 >
        sage: fundamental_group(f).sorted_presentation()
        Finitely presented group < x0, x1 | x1^-1*x0^-1*x1^-1*x0*x1*x0 >

    ::

        sage: # needs sirocco
        sage: from sage.schemes.curves.zariski_vankampen import fundamental_group
        sage: R.<x, y> = QQ[]
        sage: f = y^3 + x^3
        sage: fundamental_group(f).sorted_presentation()
        Finitely presented group < x0, x1, x2 | x2^-1*x1^-1*x0^-1*x2*x0*x1,
                                                x2^-1*x1^-1*x2*x0*x1*x0^-1 >

    It is also possible to have coefficients in a number field with a
    fixed embedding in `\QQbar`::

        sage: from sage.schemes.curves.zariski_vankampen import fundamental_group
        sage: zeta = QQbar['x']('x^2 + x+ 1').roots(multiplicities=False)[0]
        sage: zeta
        -0.50000000000000000? - 0.866025403784439?*I
        sage: F = NumberField(zeta.minpoly(), 'zeta', embedding=zeta)
        sage: F.inject_variables()
        Defining zeta
        sage: R.<x, y> = F[]
        sage: f = y^3 + x^3 + zeta * x + 1
        sage: fundamental_group(f)                                  # needs sirocco
        Finitely presented group < x0 |  >

    We compute the fundamental group of the complement of a
    quartic using the ``puiseux`` option::

        sage: # optional - sirocco
        sage: from sage.schemes.curves.zariski_vankampen import fundamental_group
        sage: R.<x, y> = QQ[]
        sage: f = x^2 * y^2 + x^2 + y^2 - 2 * x * y  * (x + y + 1)
        sage: g = fundamental_group(f); g.sorted_presentation()
        Finitely presented group < x0, x1 | x1^-2*x0^2, (x1^-1*x0)^3 >
        sage: g = fundamental_group(f, projective=True)
        sage: g.order(), g.abelian_invariants()
        (12, (4,))
        sage: fundamental_group(y * (y - 1))
        Finitely presented group < x0, x1 |  >
    """
    g = f
    x, y = g.parent().gens()
    F = g.parent().base_ring()
    d = g.degree(y)
    while not g.coefficient(y**d) in F:
        g = g.subs({x: x + y})
        d = g.degree(y)
    if projective:
        while g.degree(y) < g.degree():
            g = g.subs({x: x + y})
    bm = braid_monodromy(g)[0]
    if not bm:
        d = g.degree(y)
    else:
        d = bm[0].parent().strands()
    return fundamental_group_from_braid_mon(bm, degree=d,
                                            simplified=simplified,
                                            projective=projective,
                                            puiseux=puiseux)


def fundamental_group_arrangement(flist, simplified=True, projective=False,
                                  puiseux=True, vertical=False,
                                  braid_data=None):
    r"""
    Compute the fundamental group of the complement of a curve
    defined by a list of polynomials with the extra information
    about the correspondence of the generators
    and meridians of the elements of the list.

    INPUT:

    - ``flist`` -- a  tuple of polynomial with two variables, over a number
      field with an embedding in the complex numbers

    - ``simplified`` -- boolean (default: ``True``); if set to ``True`` the
      presentation will be simplified (see below)

    - ``projective`` -- boolean (default: ``False``); if set to ``True``,
      the fundamental group of the complement of the projective completion
      of the curve will be computed, otherwise, the fundamental group of
      the complement in the affine plane will be computed

    - ``puiseux`` -- boolean (default: ``True``); if set to ``True``
      a presentation of the fundamental group with the homotopy type
      of the complement of the affine curve will be computed, adding
      one relation if ``projective`` is set to ``True``.

    - ``vertical`` -- boolean (default: ``False``); if set to ``True``,
      whenever no curve has vertical asymptotes the computation of braid
      monodromy is simpler if some lines are vertical

    - ``braid_data`` -- tuple (default: ``None``); if it is not the default
      it is the output of ``fundamental_group_from_braid_mon`` previously
      computed

    OUTPUT:

    - A list of braids. The braids correspond to paths based in the same point;
      each of this paths is the conjugated of a loop around one of the points
      in the discriminant of the projection of ``f``.

    - A dictionary attaching to ``j`` a tuple a list of elements
      of the group  which are meridians of the curve in position ``j``.
      If ``projective`` is ``False`` and the `y`-degree of the horizontal
      components coincide with the total degree, another key is added
      to give a meridian of the line at infinity.

    EXAMPLES::

        sage: # needs sirocco
        sage: from sage.schemes.curves.zariski_vankampen import braid_monodromy
        sage: from sage.schemes.curves.zariski_vankampen import fundamental_group_arrangement
        sage: R.<x, y> = QQ[]
        sage: flist = [x^2 - y^3, x + 3 * y - 5]
        sage: g, dic = fundamental_group_arrangement(flist)
        sage: g.sorted_presentation()
        Finitely presented group
         < x0, x1, x2 | x2^-1*x1^-1*x2*x1, x2^-1*x0^-1*x2^-1*x0*x2*x0,
                        x1^-1*x0^-1*x1*x0 >
        sage: dic
        {0: [x0, x2], 1: [x1], 2: [x0^-1*x2^-1*x1^-1*x0^-1]}
        sage: g, dic = fundamental_group_arrangement(flist, simplified=False, puiseux=False)
        sage: g.sorted_presentation(), dic
        (Finitely presented group
         < x0, x1, x2, x3 | 1, 1, 1, 1, 1, 1, 1,
                            x3^-1*x2^-1*x1^-1*x2*x3*x2^-1*x1*x2,
                            x3^-1*x2^-1*x1^-1*x0^-1*x1*x2*x3*x2,
                            x3^-1*x2^-1*x1^-1*x0^-1*x1*x2*x1^-1*x0*x1*x2,
                            x3^-1*x2^-1*x1^-1*x2*x3*x2^-1*x1*x2,
                            x3^-1*x1^-1*x0*x1,
                            x1^-1*x0^-1*x1*x0, x1^-1*x0^-1*x1*x0,
                            x1^-1*x0^-1*x1*x0, x1^-1*x0^-1*x1*x0 >,
         {0: [x0, x2, x3], 1: [x1], 2: [x3^-1*x2^-1*x1^-1*x0^-1]})
        sage: fundamental_group_arrangement(flist, projective=True)
        (Finitely presented group < x |  >, {0: [x], 1: [x^-3]})
        sage: fundamental_group_arrangement([])
        (Finitely presented group <  |  >, {})
        sage: g, dic = fundamental_group_arrangement([x * y])
        sage: g.sorted_presentation(), dic
        (Finitely presented group < x0, x1 | x1^-1*x0^-1*x1*x0 >,
         {0: [x0, x1], 1: [x1^-1*x0^-1]})
        sage: fundamental_group_arrangement([y + x^2])
        (Finitely presented group < x |  >, {0: [x]})
        sage: fundamental_group_arrangement([y^2 + x], projective=True)
        (Finitely presented group < x | x^2 >, {0: [x]})
        sage: L = [x, y, x - 1, x -y]
        sage: G, dic =fundamental_group_arrangement(L)
        sage: G.sorted_presentation()
        Finitely presented group
        < x0, x1, x2, x3 | x3^-1*x2^-1*x3*x2, x3^-1*x1^-1*x0^-1*x1*x3*x0,
                           x3^-1*x1^-1*x3*x0*x1*x0^-1, x2^-1*x0^-1*x2*x0 >
        sage: dic
        {0: [x1], 1: [x3], 2: [x2], 3: [x0], 4: [x3^-1*x2^-1*x1^-1*x0^-1]}
        sage: fundamental_group_arrangement(L, vertical=True)
        (Finitely presented group
         < x0, x1, x2, x3 | x3*x0*x3^-1*x0^-1, x3*x1*x3^-1*x1^-1,
                            x1*x2*x0*x2^-1*x1^-1*x0^-1,
                            x1*x2*x0*x1^-1*x0^-1*x2^-1 >,
         {0: [x2], 1: [x0], 2: [x3], 3: [x1], 4: [x3^-1*x2^-1*x1^-1*x0^-1]})
    """
    if flist:
        f = prod(flist)
        R = f.parent()
    else:
        R = PolynomialRing(QQ, ('x', 'y'))
        f = R(1)
    x, y = R.gens()
    flist1 = tuple(flist)
    if vertical and vertical_lines_in_braidmon(flist1):
        infinity = all(Curve(g).is_vertical_line() or
                       g.degree(y) == g.degree() for g in flist1)
    else:
        infinity = any(Curve(g).has_vertical_asymptote() or
                       Curve(g).is_vertical_line() for g in flist1)
        if not infinity:
            infinity = all(g.degree(y) == g.degree() for g in flist1)
    if braid_data:
        bm, dic, dv, d1 = braid_data
    elif not flist:
        bm = []
        dic = {}
        dv = {j: j for j, f in flist1}
        d1 = 0
    else:
        bm, dic, dv, d1 = braid_monodromy(f, flist1, vertical=vertical)
    vert_lines = list(dv)
    vert_lines.sort()
    for i, j in enumerate(vert_lines):
        dic[d1 + i] = dv[j]
    g = fundamental_group_from_braid_mon(bm, degree=d1, simplified=False,
                                         projective=projective,
                                         puiseux=puiseux,
                                         vertical=vert_lines)
    if simplified:
        hom = g.simplification_isomorphism()
    else:
        hom = g.hom(codomain=g, im_gens=list(g.gens()), check=False)
    g1 = hom.codomain()
    if not flist:
        return (g1, {})
    dic1 = {}
    for i in range(len(flist1)):
        L = [j1 for j1 in dic if dic[j1] == i]
        dic1[i] = [hom(g.gen(j)) for j in L]
    if not projective and infinity:
        t = prod(hom(a) for a in g.gens()).inverse()
        dic1[len(flist1)] = [t]
    n = g1.ngens()
    rels = [rel.Tietze() for rel in g1.relations()]
    g1 = FreeGroup(n) / rels
    dic1 = {i: list({g1(el.Tietze()) for el in dic1[i]}) for i in dic1}
    return (g1, dic1)