def trinomiales_general(A, Vf, Vg, Vhs, Vd, total=True):
    r"""
    Return data necessary to compute the zeta functions of
    a triple singularity defined by

        .. math::
            z^3 -3 \prod_{i=1}^h h_i^{n_i} f z^2
            + \prod_{i=1}^h h_i g,

    where :math:`f, g, h_i` are polynomials in two variables,
    and :math:`n_i \geq 1`.

    We consider an embedded resolution of :math:`fg(f^3 - g^2)` with
    a stratification of the exceptional divisor:

    - :math:`\check{E}_i`: the smooth part of the total divisor

    - :math:`E_i\cap E_j`: the double points of the divisor

    - :math:`E_i\cap \{f=0\}`: the intersections with the strict
      transform of :math:`f`

    - :math:`E_i\cap \{g=0\}`: the intersections with the strict
      transform of :math:`g`

    - :math:`E_i\cap \{f^3=g^2\}`: the intersections with the strict
      transform of :math:`f^3-g^2`

    It is also possible to forget the strict transforms of :math:`f, g`
    and even to consider only resolutions of :math:`f^3 - g^2`.

    INPUT:

    - ``A`` -- the incidence matrix of a resolution

    - ``Vf`` -- the vector with the intersection numbers of the strict
      transform of :math:`f` with the divisors of the resolution

    - ``Vg`` -- the vector with the intersection numbers of the strict
      transform of :math:`g` with the divisors of the resolution

    - ``Vhs`` -- a list of tuple whose elements are tuples `(h_i, n_i)`

    - ``Vd`` -- the vector with the intersection numbers of the strict
      transform of :math:`f^3 - g^2` with the divisors of the resolution

    - ``total`` -- boolean (default: ``True``), whether we consider also
      the strata corresponding to the strict transforms of :math:`f, g`.

    OUTPUT:

    Polynomials :math:`z^3 + M_1 z + M_0` or :math:`z^3 + M_2 z^2 + M_0` are
    assigned to each strata. The monomials are of the form :math:`x^k` for the
    smooth parts of the divisors and :math:`x^k y^l` for the double points.
    Each variable is associated to a divisor or to a strict transform

    - ``n`` -- the number of irreducible components of the exceptional divisor

    - ``Nf`` -- the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`f`

    - ``Ng`` -- the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`g`

    - ``Nhs`` -- a list of tuple whose elements are tuples, the first element
      is the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`h_i` and the second is :math:`n_i`

    - ``Nd`` -- the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`f^3 - g^2`

    - ``euler`` -- the vector in :math:`\mathbb{Z}^n` of the Euler
      characteristic of the smooth part of the strata. If ``total=True`` this
      smooth part does not contain the intersections with total transform
      of :math:`f, g`

    - ``nus`` -- the vector in :math:`\mathbb{Z}^n` of the canonical
      multiplicities

    - ``pols`` -- a dictionary with keys ``(i, )``, ``(i, j)``, ``(i, 'f')``
      (if ``total=True``), ``(i, 'g')`` (if ``total=True``), and ``(i, 'd')``
      assigning the polynomial of the stratum (non-degenerated with respect
      to its Newton polytope), the weights ``(nui, nuj, 1)`` coding
      the form
      :math:`x^{\nu_i-1}y^{\nu_j-1}\operatorname{d}\! x \operatorname{d}\! y \operatorname{d}\! z`

    EXAMPLES::

        sage: R.<x, y, z> = QQ[]
        sage: f, g = (x, y)
        sage: A = Matrix(QQ, 3, [-3, 0, 1, 0, -2, 1, 1, 1, -1])
        sage: Vf = vector(QQ, [1, 0, 0])
        sage: Vg = vector(QQ, [0, 1, 0])
        sage: Vd = vector(QQ, [0, 0, 1])
        sage: trinomiales(A, Vf, Vg, Vd)
        (3, (1, 1, 2), (1, 2, 3), (2, 3, 6), (0, 0, -1), (2, 3, 5),
         {(0, 'f'): (x + x*y*z + z^3, (2, 1, 1), 1),
          (0, 2): (z^3 + x*y^3 + x*y^2*z, (2, 5, 1), 1),
          (1, 'g'): (x*z + x^2*y + z^3, (3, 1, 1), 1),
          (1, 2): (z^3 + x*y^2*z + x^2*y^3, (3, 5, 1), 1),
          (2, 'd'): (x*z^2 + z^3 + x^3*y, (5, 1, 1), 1),
          0: (x + x*z + z^3, (2, 1, 1), 0),
          1: (x^2 + x*z + z^3, (3, 1, 1), 0),
          2: (x^3 + x^2*z + z^3, (5, 1, 1), -1)})
        sage: trinomiales(A, Vf, Vg, Vd, total=False)
        (3, (1, 1, 2), (1, 2, 3), (2, 3, 6), (1, 1, -1), (2, 3, 5),
         {(0, 2): (z^3 + x*y^3 + x*y^2*z, (2, 5, 1), 1),
          (1, 2): (z^3 + x*y^2*z + x^2*y^3, (3, 5, 1), 1),
          (2, 'd'): (x*z^2 + z^3 + x^3*y, (5, 1, 1), 1),
          0: (x + x*z + z^3, (2, 1, 1), 1),
          1: (x^2 + x*z + z^3, (3, 1, 1), 1),
          2: (x^3 + x^2*z + z^3, (5, 1, 1), -1)})
    """
    R = PolynomialRing(QQ, 'x, y, z', order='negdeglex')
    x, y, z = R.gens()
    n = A.nrows()
    A1 = -A ** -1
    Nf, Ng, Nd = [A1 * N0 for N0 in (Vf, Vg, Vd)]
    Nhs = [(A1 * N0, h) for N0, h in Vhs]
    ch = len(Nhs)
    NNf = Nf + sum(j * V0 for V0, j in Nhs)
    NNg = Ng + sum(V0 for V0, j in Nhs)
    NNd = Nd + 2 * sum(V0 for V0, j in Nhs)
    diag = vector(A.diagonal())
    unos = vector(n * [1])
    if total:
        VV = Vf + Vg
    else:
        VV = Vf.parent()(0)
    resta = VV + Vd + sum(V0 for V0, m in Vhs)
    euler = (-A + 2 * identity_matrix(n)) * unos + diag - resta
    nus = A1 * (diag + 2 * unos) + unos
    expos = {}
    for j in range(n):
        a, b = (NNf[j], NNg[j])
        eq = 3 * a == 2 * NNg[j]
        nnd = min(3 * a, 2 * b) == NNd[j]
        if nnd:
            expos[j] = ((a, b), (a / 2, b), nus[j], eq, nnd)
        else:
            expos[j] = ((a, b), (a / 2, NNd[j] - b), nus[j], eq, nnd)
    cadenas = ('d', )
    cadenash = tuple(('h', h) for h in range(ch))
    if total:
        expos['f'] = ((1, 0), (0, 0), 1, False, True)
        expos['g'] = ((0, 1), (0, 0), 1, False, True)
        cadenas = ('f', 'g') + cadenas
    expos['d'] = ((0, 0), (0, 1), 1, True, False)
    for h, V0 in enumerate(Vhs):
        expos['h', h] = ((V0[1], 1), (0, 0), 1, False, True)
    pols = {}
    for j in range(n):
        if expos[j][-1]:
            a, b = expos[j][0]
            pols[j] = (z ** 3 + x ** a * z + x ** b, (expos[j][2], 1, 1),
                       euler[j])
        else:
            a, b = expos[j][1]
            pols[j] = (z ** 3 + x ** a * z ** 2 + x ** b,
                       (expos[j][2], 1, 1), euler[j])
    interdiv = [(i, j) for i, j in Combinations(n, 2) if A[i, j] != 0]
    interstrict = [(i, j) for i in range(n) for j in cadenas
                   if eval('V' + j)[i] != 0]
    interh = [(i, j) for i in range(n) for j in cadenash if Vhs[j[1]][0][i] != 0]
    for i, j in interdiv + interstrict + interh:
        if isinstance(j, str):
            ep = eval('V' + j)[i]
        else:
            ep = 1
        eqi, nndi = expos[i][-2:]
        eqj, nndj = expos[j][-2:]
        nd = eqi and eqj and not (nndi and nndj)
        if not nd:
            ai, bi = expos[i][0]
            aj, bj = expos[j][0]
            pols[i, j] = (z ** 3 + x ** ai * y ** aj * z + x ** bi * y ** bj,
                          (expos[i][2], expos[j][2], 1), ep)
        else:
            ai, bi = expos[i][1]
            aj, bj = expos[j][1]
            pols[i, j] = (z ** 3 + x ** ai * y ** aj * z ** 2 +
                          x ** bi * y ** bj, (expos[i][2], expos[j][2], 1), ep)
    return n, Nf, Ng, Nd, Nhs, euler, nus, pols


def estrato(pol, pesos, euler, info=False):
    r"""
    Return the topological and Igusa zeta functions associated to a stratum

    INPUT:

    - ``pol`` -- a trinomial non degenerated with respect to its Newton
      polyhedron

    - ``pesos`` -- the canonical weights :math:`\nu:i`

    - `euler`` -- the Euler characteristic of the stratum

    - ``info`` -- boolean (default: False), if ``True`` a dictionary giving
      the powers of cyclotomic polynomials for the monodromy zeta function
      is given

    EXAMPLES::

        sage: R.<x, y, z> = QQ[]
        sage: F = z^3 + x^4 * y^6 * z + x^7 * y^9
        sage: estrato(F, (1, 1, 1), 1, info=True)
        ((1/3*s^2 + 1/2*s + 2/9)/(s^3 + 35/18*s^2 + 7/6*s + 2/9), -t^3 + 1,
         {1: 1, 3: 1})
    """
    Z1 = ZetaFunctions(pol)
    zetatop = euler * Z1.topological_zeta(local=True, weights=pesos)
    zetamon = Z1.monodromy_zeta() ** euler
    if info:
        dic = {}
        for p, j in zetamon.factor():
            c = p.is_cyclotomic(certificate=True)
            dic[c] = j
        return zetatop, zetamon, dic
    return zetatop, zetamon


def zetas(n, pols, info=False):
    r"""
    Return the topological and Igusa zeta functions from the data of
    a stratification

    INPUT:

    - ``n`` -- the number of exceptional divisors

    - ``pols`` -- a dictionary from ``trinomiales``

    - ``info`` -- boolean (default: False), if ``True`` a dictionary giving
      the powers of cyclotomic polynomials for the monodromy zeta function
      is given

    EXAMPLES::

        sage: A = Matrix(QQ, 3, [-3, 0, 1, 0, -2, 1, 1, 1, -1])
        sage: Vf = vector(QQ, [1, 0, 0])
        sage: Vg = vector(QQ, [0, 1, 0])
        sage: Vd = vector(QQ, [0, 0, 1])
        sage: n, Nf, Ng, Nd, euler, nus, pols = trinomiales(A, Vf, Vg, Vd)
        sage: zetas(n, pols, info=True)
        (1/(s + 1), -t + 1, {1: 1})
    """
    zetatop = 0
    T = PolynomialRing(QQ, 't')
    t = T.gen(0)
    zetamon = 1 / (1 - t ** 3) ** n
    for k in pols:
        zt, zm = estrato(*pols[k])
        zetatop += zt
        zetamon *= zm
    if info:
        dic = {}
        for p, j in zetamon.factor():
            c = p.is_cyclotomic(certificate=True)
            dic[c] = j
        return zetatop, zetamon, dic
    return zetatop, zetamon


def zetas_components(n, pols, A, Vf, Vg, Vd, Nf, Ng, Nd, info=False, debug=False):
    r"""
    Return the topological and Igusa zeta functions from the data of
    a stratification

    INPUT:

    - ``n`` -- the number of exceptional divisors

    - ``pols`` -- a dictionary from ``trinomiales``

    - ``A`` -- the incidence matrix of the resolution

    - ``Vf`` -- the vector with the intersection numbers of the strict
      transform of :math:`f` with the divisors of the resolution

    - ``Vg`` -- the vector with the intersection numbers of the strict
      transform of :math:`g` with the divisors of the resolution

    - ``Vd`` -- the vector with the intersection numbers of the strict
      transform of :math:`f^3 - g^2` with the divisors of the resolution

    - ``Nf`` -- the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`f`

    - ``Ng`` -- the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`g`

    - ``Nd`` -- the vector in :math:`\mathbb{Z}^n` of the multiplicities for
      the total transform of :math:`f^3 - g^2`

  - ``info`` -- boolean (default: False), if ``True`` a dictionary giving
      the powers of cyclotomic polynomials for the monodromy zeta function
      is given

    EXAMPLES::

        sage: A = Matrix(QQ, 3, [-3, 0, 1, 0, -2, 1, 1, 1, -1])
        sage: Vf = vector(QQ, [1, 0, 0])
        sage: Vg = vector(QQ, [0, 1, 0])
        sage: Vd = vector(QQ, [0, 0, 1])
        sage: n, Nf, Ng, Nd, euler, nus, pols = trinomiales(A, Vf, Vg, Vd)
        sage: zetas(n, pols, info=True)
        (1/(s + 1), -t + 1, {1: 1})
    """
    T = PolynomialRing(QQ, 't')
    t = T.gen(0)
    pares = [a for a in pols if a not in ZZ]
    ceros = [i for i in range(n) if Nd[i] < 3 * Nf[i]]
    Aceros = A.matrix_from_rows_and_columns(ceros, ceros)
    grafoceros = Graph(Aceros)
    grafoceros.relabel(ceros)
    componentesceros = grafoceros.connected_components(sort=False)
    for j, cpt in enumerate(componentesceros):
        cpt2 = [(i, j) for i, j in pares if i in cpt or j in cpt]
        componentesceros[j] = cpt + cpt2
    if debug:
        print('componentes ceros', componentesceros)
    polos = [i for i in range(n) if Nd[i] < 2 * Ng[i]]
    Apolos = A.matrix_from_rows_and_columns(polos, polos)
    grafopolos = Graph(Apolos)
    grafopolos.relabel(polos)
    componentespolos = grafopolos.connected_components(sort=False)
    for j, cpt in enumerate(componentespolos):
        cpt2 = [(i, j) for i, j in pares if i in cpt or j in cpt]
        componentespolos[j] = cpt + cpt2
    if debug:
        print('componentes polos', componentespolos)
    unos = [i for i in range(n) if Nd[i] > 3 * Nf[i]]
    Aunos = A.matrix_from_rows_and_columns(unos, unos)
    grafounos = Graph(Aunos)
    grafounos.relabel(unos)
    componentesunos = grafounos.connected_components(sort=False)
    for j, cpt in enumerate(componentesunos):
        cpt2 = [(i, j) for i, j in pares if i in cpt or j in cpt]
        componentesunos[j] = cpt + cpt2
    resto = [i for i in range(n) if Nd[i] == 2 * Ng[i] == 3 * Nf[i]]
    dic = [i for i in resto if Vd[i] > 0 or [j for j in unos if A[i, j] > 0]]
    componentesunos += [[(i, j)] for i, j in pares if i in dic and j == 'd']
    if debug:
        print('componentes unos', componentesunos)
        print('dicríticos ', dic)
    cte = [i for i in resto if i not in dic]
    Acte = A.matrix_from_rows_and_columns(cte, cte)
    grafocte = Graph(Acte)
    grafocte.relabel(cte)
    componentescte = grafocte.connected_components(sort=False)
    for j, cpt in enumerate(componentescte):
        cpt2 = [(i, j) for i, j in pares if i in cpt or j in cpt]
        componentescte[j] = cpt + cpt2
    componentescte += [(i, j) for i, j in pares if i in dic and j in dic]
    if debug:
        print('componentes constante', componentescte)
    res = {}
    componentes = [('ceros', tuple(cpt)) for cpt in componentesceros]
    componentes += [('polos', tuple(cpt)) for cpt in componentespolos]
    componentes += [('unos', tuple(cpt)) for cpt in componentesunos]
    componentes += [('dic', (cpt,)) for cpt in dic]
    componentes += [('cte', tuple(cpt)) for cpt in componentescte]
    for str, cpt in componentes:
        zetatop = 0
        n0 = len([i for i in cpt if i in ZZ])
        zetamon = 1 / (1 - t ** 3) ** n0
        for k in cpt:
            zt, zm = estrato(*pols[k])
            zetatop += zt
            zetamon *= zm
        if info:
            dic = {}
            for p, j in zetamon.factor():
                c = p.is_cyclotomic(certificate=True)
                dic[c] = j
            res[str, cpt] = (zetatop, zetamon, dic)
        else:
            res[str, cpt] = (zetatop, zetamon)
    return res