def sum_of_partitions(self, la):
            r"""
            Return the sum over all sets partitions whose shape is ``la``,
            scaled by `\prod_i m_i!` where `m_i` is the multiplicity
            of `i` in ``la``.

            INPUT:

            - ``la`` -- an integer partition

            OUTPUT:

            - an element of ``self``

            EXAMPLES::

                sage: w = SymmetricFunctionsNonCommutingVariables(QQ).dual().w()
                sage: w.sum_of_partitions([2,1,1])
                2*w{{1}, {2}, {3, 4}} + 2*w{{1}, {2, 3}, {4}} + 2*w{{1}, {2, 4}, {3}}
                 + 2*w{{1, 2}, {3}, {4}} + 2*w{{1, 3}, {2}, {4}} + 2*w{{1, 4}, {2}, {3}}
            """
            la = Partition(la)
            c = prod([factorial(_) for _ in la.to_exp()])
            P = SetPartitions()
            return self.sum_of_terms([(P(m), c) for m in SetPartitions(sum(la), la)], distinct=True)

            def is_symmetric(self):
                r"""
                Determine if a `NCSym^*` function, expressed in the
                `\mathbf{w}` basis, is symmetric.

                A function `f` in the `\mathbf{w}` basis is a symmetric
                function if it is in the image of `\chi^*`. That is to say we
                have

                .. MATH::

                    f = \sum_{\lambda} c_{\lambda} \prod_i m_i(\lambda)!
                    \sum_{\lambda(A) = \lambda} \mathbf{w}_A

                where the second sum is over all set partitions `A` whose
                shape `\lambda(A)` is equal to `\lambda` and `m_i(\mu)` is
                the multiplicity of `i` in the partition `\mu`.

                OUTPUT:

                - ``True`` if `\lambda(A)=\lambda(B)` implies the coefficients of
                  `\mathbf{w}_A` and `\mathbf{w}_B` are equal, ``False`` otherwise

                EXAMPLES::

                    sage: w = SymmetricFunctionsNonCommutingVariables(QQ).dual().w()
                    sage: elt = w.sum_of_partitions([2,1,1])
                    sage: elt.is_symmetric()
                    True
                    sage: elt -= 3*w.sum_of_partitions([1,1])
                    sage: elt.is_symmetric()
                    True
                    sage: w = SymmetricFunctionsNonCommutingVariables(ZZ).dual().w()
                    sage: elt = w.sum_of_partitions([2,1,1]) / 2
                    sage: elt.is_symmetric()
                    False
                    sage: elt = w[[1,3],[2]]
                    sage: elt.is_symmetric()
                    False
                    sage: elt = w[[1],[2,3]] + w[[1,2],[3]] + 2*w[[1,3],[2]]
                    sage: elt.is_symmetric()
                    False
                """
                d = {}
                R = self.base_ring()
                for A, coeff in self:
                    la = A.shape()
                    exp = prod([factorial(_) for _ in la.to_exp()])
                    if la not in d:
                        if coeff / exp not in R:
                            return False
                        d[la] = [coeff, 1]
                    else:
                        if d[la][0] != coeff:
                            return False
                        d[la][1] += 1
                # Make sure we've seen each set partition of the shape
                return all(d[la][1] == SetPartitions(la.size(), la).cardinality() for la in d)

def logpp_binom(n, p, p_prec):
    """returns the (integral) power series p^n*(log_p(1+p*z)/log_p(1+p) choose n)"""
    #prod=1+0*z
    L = logpp_gam(p, p_prec)
    ans = prod([(L - j) for j in range(n)])
    #for j in range(0,n):
    #    prod=prod*(L-j)
    ans *= (p ** n) / factorial(n)
    
    return ps_normalize(ans.truncate(p_prec), p, p_prec)

def logp_binom(n, p, p_prec):
    """returns the (integral) power series (log_p(1+z)/log_p(1+p) choose n)"""
    #prod=1+0*z
    if n == 0:
        return PolynomialRing(QQ, 'y')(1)
    L = logp_gam(p, p_prec)
    ans = prod([(L - j) for j in range(n)])
    #for j in range(0,n):
    #    prod=prod*(L-j)
    ans = ans / factorial(n)
    
    return ps_normalize(ans.truncate(p_prec+1), p, p_prec) #Do we need the +1?

 def coeff(p, q):
     ret = QQ.one()
     last = 0
     for val in p:
         count = 0
         s = 0
         while s != val:
             s += q[last+count]
             count += 1
         ret /= factorial(count)
         last += count
     return ret

    def eval_formula(self, n, x):
        """
        Evaluate ``chebyshev_T`` using an explicit formula.
        See [ASHandbook]_ 227 (p. 782) for details for the recurions.
        See also [EffCheby]_ for fast evaluation techniques.

        INPUT:

        - ``n`` -- an integer

        - ``x`` -- a value to evaluate the polynomial at (this can be
          any ring element)

        EXAMPLES::

            sage: chebyshev_T.eval_formula(-1,x)
            x
            sage: chebyshev_T.eval_formula(0,x)
            1
            sage: chebyshev_T.eval_formula(1,x)
            x
            sage: chebyshev_T.eval_formula(2,0.1) == chebyshev_T._evalf_(2,0.1)
            True
            sage: chebyshev_T.eval_formula(10,x)
            512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
            sage: chebyshev_T.eval_algebraic(10,x).expand()
            512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
        """
        if n < 0:
            return self.eval_formula(-n, x)
        elif n == 0:
            return parent(x).one()

        res = parent(x).zero()
        for j in xrange(0, n // 2 + 1):
            f = factorial(n - 1 - j) / factorial(j) / factorial(n - 2 * j)
            res += (-1) ** j * (2 * x) ** (n - 2 * j) * f
        res *= n / 2
        return res

    def __getitem__(self, key):
        """
        EXAMPLES::

            sage: a, b, c, z = var('a b c z')
            sage: hs = HypergeometricSeries([a, b], [c], z)
            sage: for i in range(4): print hs[i]
            1
            a*b*z/c
            1/2*(b + 1)*(a + 1)*a*b*z^2/((c + 1)*c)
            1/6*(b + 1)*(b + 2)*(a + 1)*(a + 2)*a*b*z^3/((c + 1)*(c + 2)*c)

        """
        if key >= 0:
            nominator = PochhammerSymbol(self.list_a, key).evaluate()
            denominator = PochhammerSymbol(self.list_b, key).evaluate()*factorial(key)
            return nominator / denominator * self.z**key
        else:
            return 0

    def expansion_on_basis(self, w):
        r"""
        Return the expansion of `S_w` in words of the shuffle algebra.

        INPUT:

        - ``w`` -- a word

        EXAMPLES::

            sage: S = ShuffleAlgebra(QQ, 'ab').dual_pbw_basis()
            sage: S.expansion_on_basis(Word())
            B[word: ]
            sage: S.expansion_on_basis(Word()).parent()
            Shuffle Algebra on 2 generators ['a', 'b'] over Rational Field
            sage: S.expansion_on_basis(Word('abba'))
            2*B[word: aabb] + B[word: abab] + B[word: abba]
            sage: S.expansion_on_basis(Word())
            B[word: ]
            sage: S.expansion_on_basis(Word('abab'))
            2*B[word: aabb] + B[word: abab]
        """
        from sage.functions.other import factorial
        if len(w) == 0:
            return self._alg.one()
        if len(w) == 1:
           return self._alg.monomial(w)

        if w.is_lyndon():
            W = self.basis().keys()
            letter = W(w[0])
            expansion = self.expansion_on_basis(W(w[1:]))
            return self._alg.sum_of_terms([(letter * i, c) for i,c in expansion])

        lf = w.lyndon_factorization()
        powers = {}
        for i in lf:
            powers[i] = powers.get(i, 0) + 1
        denom = prod(factorial(p) for p in powers.values())
        result = self._alg.prod(self.expansion_on_basis(i) for i in lf)
        return self._alg(result / denom)

    def _eval_(self, s, x):
        r"""
        TESTS::

            sage: hurwitz_zeta(x, 1)
            zeta(x)
            sage: hurwitz_zeta(4, 3)
            1/90*pi^4 - 17/16
            sage: hurwitz_zeta(-4, x)
            -1/5*x^5 + 1/2*x^4 - 1/3*x^3 + 1/30*x
            sage: hurwitz_zeta(3, 0.5)
            8.41439832211716
        """
        if x == 1:
            return zeta(s)
        if s in ZZ and s > 1:
            return ((-1) ** s) * psi(s - 1, x) / factorial(s - 1)
        elif s in ZZ and s < 0:
            return -bernoulli_polynomial(x, -s + 1) / (-s + 1)
        else:
            return

    def order(self):
        r"""
        Returns the number of elements of ``self``.

        EXAMPLES::

            sage: F.<a> = GF(4)
            sage: SemimonomialTransformationGroup(F, 5).order() == (4-1)**5 * factorial(5) * 2
            True
        """
        from sage.functions.other import factorial
        from sage.categories.homset import End
        n = self.degree()
        R = self.base_ring()
        if R.is_field():
            multgroup_size = len(R)-1
            autgroup_size = R.degree()
        else:
            multgroup_size = R.unit_group_order()
            autgroup_size = len([x for x in End(R) if x.is_injective()])
        return multgroup_size**n * factorial(n) * autgroup_size

        def idempotent(self, la):
            """
            Return the idemponent corresponding to the partition ``la``.

            EXAMPLES::

                sage: I = DescentAlgebra(QQ, 4).I()
                sage: E = I.idempotent([3,1]); E
                1/2*I[1, 3] + 1/2*I[3, 1]
                sage: E*E == E
                True
                sage: E2 = I.idempotent([2,1,1]); E2
                1/6*I[1, 1, 2] + 1/6*I[1, 2, 1] + 1/6*I[2, 1, 1]
                sage: E2*E2 == E2
                True
                sage: E*E2 == I.zero()
                True
            """
            from sage.combinat.permutation import Permutations
            k = len(la)
            C = Compositions(self.realization_of()._n)
            return self.sum_of_terms([(C(x), ~QQ(factorial(k))) for x in Permutations(la)])

            def to_symmetric_function(self):
                r"""
                Take a function in the `\mathbf{w}` basis, and return its
                symmetric realization, when possible, expressed in the
                homogeneous basis of symmetric functions.

                OUTPUT:

                - If ``self`` is a symmetric function, then the expansion
                  in the homogeneous basis of the symmetric functions is returned.
                  Otherwise an error is raised.

                EXAMPLES::

                    sage: w = SymmetricFunctionsNonCommutingVariables(QQ).dual().w()
                    sage: elt = w[[1],[2,3]] + w[[1,2],[3]] + w[[1,3],[2]]
                    sage: elt.to_symmetric_function()
                    h[2, 1]
                    sage: elt = w.sum_of_partitions([2,1,1]) / 2
                    sage: elt.to_symmetric_function()
                    1/2*h[2, 1, 1]

                TESTS::

                    sage: w = SymmetricFunctionsNonCommutingVariables(QQ).dual().w()
                    sage: w(0).to_symmetric_function()
                    0
                    sage: w([]).to_symmetric_function()
                    h[]
                    sage: (2*w([])).to_symmetric_function()
                    2*h[]
                """
                if not self.is_symmetric():
                    raise ValueError("not a symmetric function")
                h = SymmetricFunctions(self.parent().base_ring()).homogeneous()
                d = {A.shape(): c for A, c in self}
                return h.sum_of_terms(
                    [(AA, cc / prod([factorial(_) for _ in AA.to_exp()])) for AA, cc in d.items()], distinct=True
                )

def coeff_sp(J,I):
    r"""
    Returns the coefficient `sp_{J,I}` as defined in [NCSF]_.

    INPUT:

    - ``J`` -- a composition
    - ``I`` -- a composition refining ``J``

    OUTPUT:

    - integer

    EXAMPLES::

        sage: from sage.combinat.ncsf_qsym.combinatorics import coeff_sp
        sage: coeff_sp(Composition([1,1,1]), Composition([2,1]))
        2
        sage: coeff_sp(Composition([2,1]), Composition([3]))
        4
    """
    return prod(factorial(len(K))*prod(K) for K in J.refinement_splitting(I))

    def _eval_(self, s, x):
        r"""
        TESTS::

            sage: hurwitz_zeta(x, 1)
            zeta(x)
            sage: hurwitz_zeta(4, 3)
            1/90*pi^4 - 17/16
            sage: hurwitz_zeta(-4, x)
            -1/5*x^5 + 1/2*x^4 - 1/3*x^3 + 1/30*x
            sage: hurwitz_zeta(3, 0.5)
            8.41439832211716
        """
        co = get_coercion_model().canonical_coercion(s, x)[0]
        if is_inexact(co) and not isinstance(co, Expression):
            return self._evalf_(s, x, parent=parent(co))
        if x == 1:
            return zeta(s)
        if s in ZZ and s > 1:
            return ((-1) ** s) * psi(s - 1, x) / factorial(s - 1)
        elif s in ZZ and s < 0:
            return -bernoulli_polynomial(x, -s + 1) / (-s + 1)
        else:
            return

def kontsevich_star_product_terms(K, prime_weights, precision):
    """
    Kontsevich star product terms in KontsevichGraphSums ``K`` up to order ``precision``.

    INPUT:

    - ``K`` - KontsevichGraphSums module.
    - ``prime_weights`` - weights of prime graphs, modulo edge labeling and mirror images.

    EXAMPLES::

        sage: K = KontsevichGraphSums(QQ);
        sage: weights = {}
        sage: weights[KontsevichGraph({'F' : {},'G' : {}}, ground_vertices=('F','G'), immutable=True)] = 1
        sage: weights[KontsevichGraph([(1, 'F', 'L'), (1, 'G', 'R')], ground_vertices=('F','G'), immutable=True)] = 1/2
        sage: weights[KontsevichGraph([(1, 2, 'R'), (1, 'F', 'L'), (2, 1, 'L'), (2, 'G', 'R')], ground_vertices=('F','G'), immutable=True)] = 1/24
        sage: weights[KontsevichGraph([(1, 2, 'R'), (1, 'F', 'L'), (2, 'F', 'L'), (2, 'G', 'R')], ground_vertices=('F','G'), immutable=True)] = 1/12
        sage: S.<h> = KontsevichGraphSeriesRng(K, star_product_terms = kontsevich_star_product_terms(K, weights, 2), default_prec = 2);
        sage: F = S(KontsevichGraph(('F',),immutable=True));
        sage: G = S(KontsevichGraph(('G',),immutable=True));
        sage: H = S(KontsevichGraph(('H',),immutable=True));
        sage: A = (F*G)*H - F*(G*H)
        sage: A.reduce()
        sage: len(A[2])  # three terms in the Jacobi identity
        3
    """
    series_terms = {0 : K([(prime_weights[KontsevichGraph(('F','G'), immutable=True)],
                            KontsevichGraph(('F','G'), immutable=True))])}
    for n in range(1, precision + 1):
        term = K(0)
        for graph in kontsevich_graphs(n, modulo_edge_labeling=True, positive_differential_order=True):
            coeff = kontsevich_weight(graph, prime_weights)*graph.multiplicity()/factorial(len(graph.internal_vertices()))
            if coeff != 0:
                term += K([(coeff, graph)])
        series_terms[n] = term
    return series_terms

#.........这里部分代码省略.........
        [X_1 + X_2, 1/2*X_12, 1/12*X_112 + 1/12*X_122]
        sage: [Z for Z in bch_iterator(X_1 + X_2, X_12)]
        [X_1 + X_2 + X_12, 1/2*X_112 - 1/2*X_122, 0]

    The elements ``X`` and ``Y`` don't need to be elements of the same Lie
    algebra if there is a coercion from one to the other::

        sage: L = LieAlgebra(QQ, 3, step=2)
        sage: L.inject_variables()
        Defining X_1, X_2, X_3, X_12, X_13, X_23
        sage: S = L.subalgebra(X_1, X_2)
        sage: bch1 = [Z for Z in bch_iterator(S(X_1), S(X_2))]; bch1
        [X_1 + X_2, 1/2*X_12]
        sage: bch1[0].parent() == S
        True
        sage: bch2 = [Z for Z in bch_iterator(S(X_1), X_3)]; bch2
        [X_1 + X_3, 1/2*X_13]
        sage: bch2[0].parent() == L
        True

    The BCH formula requires a coercion from the rationals::

        sage: L.<X,Y,Z> = LieAlgebra(ZZ, 2, step=2)
        sage: bch = bch_iterator(X, Y); next(bch)
        Traceback (most recent call last):
        ...
        TypeError: the BCH formula is not well defined since Integer Ring has no coercion from Rational Field

    TESTS:

    Compare to the BCH formula up to degree 5 given by wikipedia::

        sage: from sage.algebras.lie_algebras.bch import bch_iterator
        sage: bch = bch_iterator()
        sage: L.<X,Y> = LieAlgebra(QQ)
        sage: L = L.Lyndon()
        sage: computed_BCH = L.sum(next(bch) for k in range(5))
        sage: wikiBCH = X + Y + 1/2*L[X,Y] + 1/12*(L[X,[X,Y]] + L[Y,[Y,X]])
        sage: wikiBCH += -1/24*L[Y,[X,[X,Y]]]
        sage: wikiBCH += -1/720*(L[Y,[Y,[Y,[Y,X]]]] + L[X,[X,[X,[X,Y]]]])
        sage: wikiBCH += 1/360*(L[X,[Y,[Y,[Y,X]]]] + L[Y,[X,[X,[X,Y]]]])
        sage: wikiBCH += 1/120*(L[Y,[X,[Y,[X,Y]]]] + L[X,[Y,[X,[Y,X]]]])
        sage: computed_BCH == wikiBCH
        True

    ALGORITHM:

    The BCH formula `\log(\exp(X)\exp(Y)) = \sum_k Z_k` is computed starting
    from `Z_1 = X + Y`, by the recursion

    .. MATH::

        (m+1)Z_{m+1} =  \frac{1}{2}[X - Y, Z_m]
        + \sum_{2\leq 2p \leq m}\frac{B_{2p}}{(2p)!}\sum_{k_1+\cdots+k_{2p}=m}
        [Z_{k_1}, [\cdots [Z_{k_{2p}}, X + Y]\cdots],

    where `B_{2p}` are the Bernoulli numbers, see Lemma 2.15.3. in [Var1984]_.

    .. WARNING::

        The time needed to compute each successive term increases exponentially.
        For example on one machine iterating through `Z_{11},...,Z_{18}` for a
        free Lie algebra, computing each successive term took 4-5 times longer,
        going from 0.1s for `Z_{11}` to 21 minutes for `Z_{18}`.
    """
    if X is None or Y is None:
        L = LieAlgebra(QQ, ['X', 'Y']).Lyndon()
        X, Y = L.lie_algebra_generators()
    else:
        X, Y = canonical_coercion(X, Y)
        L = X.parent()

    R = L.base_ring()
    if not R.has_coerce_map_from(QQ):
        raise TypeError("the BCH formula is not well defined since %s "
                        "has no coercion from %s" % (R, QQ))

    xdif = X - Y
    Z = [0, X + Y]  # 1-based indexing for convenience
    m = 1
    yield Z[1]

    while True:
        m += 1
        if L in LieAlgebras.Nilpotent and m > L.step():
            return

        # apply the recursion formula of [Var1984]
        Zm = ~QQ(2 * m) * xdif.bracket(Z[-1])
        for p in range(1, (m - 1) // 2 + 1):
            partitions = IntegerListsLex(m - 1, length=2 * p, min_part=1)
            coeff = bernoulli(2 * p) / QQ(m * factorial(2 * p))
            for kvec in partitions:
                W = Z[1]
                for k in kvec:
                    W = Z[k].bracket(W)
                Zm += coeff * W

        Z.append(Zm)
        yield Zm

def number_of_partial_injection(n, algorithm='binomial'):
    r"""
    Return the number of partial injections on an set of `n` elements
    defined on a subset of `k` elements for each `k` in `0, 1, ..., n`.

    INPUT:

    - ``n`` -- integer
    - ``algorithm`` -- string (default: ``'binomial'``), ``'binomial'``
      or ``'recursive'``. When n>50, the binomial coefficient approach is
      faster (linear time vs quadratic time).

    OUTPUT:

        list

    .. NOTE::

        The recursive code of this function was originally written by
        Vincent Delecroix (Nov 30, 2017) the day after a discussion with
        Pascal Weil and me at LaBRI.

    EXAMPLES::

        sage: from slabbe import number_of_partial_injection
        sage: number_of_partial_injection(0)
        [1]
        sage: number_of_partial_injection(1)
        [1, 1]
        sage: number_of_partial_injection(2)
        [1, 4, 2]
        sage: number_of_partial_injection(3)
        [1, 9, 18, 6]
        sage: number_of_partial_injection(4)
        [1, 16, 72, 96, 24]
        sage: number_of_partial_injection(5)
        [1, 25, 200, 600, 600, 120]
        sage: number_of_partial_injection(6)
        [1, 36, 450, 2400, 5400, 4320, 720]
        sage: number_of_partial_injection(7)
        [1, 49, 882, 7350, 29400, 52920, 35280, 5040]
        sage: number_of_partial_injection(8)
        [1, 64, 1568, 18816, 117600, 376320, 564480, 322560, 40320]

    TESTS::

        sage: number_of_partial_injection(8, algorithm='recursive')
        [1, 64, 1568, 18816, 117600, 376320, 564480, 322560, 40320]

    REFERENCE:

        https://oeis.org/A144084
    """
    if algorithm == 'binomial':
        return [binomial(n,k)**2*factorial(k) for k in range(n+1)]
    elif algorithm == 'recursive':
        L = [ZZ(1)]
        for t in range(1, n+1):
            L.append(t * L[-1])
            for k in range(t-1, 0, -1):
                L[k] = (t * L[k]) // (t-k) + t * L[k-1]
        return L

#.........这里部分代码省略.........
            f0
            sage: f = ScalarField(m, f0, name='f')
            sage: sf = f.hodge_star(g) ; sf 
            4-form '*f' on the 4-dimensional manifold 'M'
            sage: sf.view()
            *f = f0 dt/\dx/\dy/\dz
            sage: ssf = sf.hodge_star(g) ; ssf
            scalar field '**f' on the 4-dimensional manifold 'M'
            sage: ssf.view()
            **f: (t, x, y, z) |--> -f0
            sage: ssf == -f  # must hold for a Lorentzian metric             
            True

        Hodge star of a 1-form in Minkowksi spacetime::
        
            sage: a = OneForm(m, 'A')
            sage: var('At Ax Ay Az')
            (At, Ax, Ay, Az)
            sage: a[:] = (At, Ax, Ay, Az)
            sage: a.view()
            A = At dt + Ax dx + Ay dy + Az dz
            sage: sa = a.hodge_star(g) ; sa
            3-form '*A' on the 4-dimensional manifold 'M'
            sage: sa.view()
            *A = -Az dt/\dx/\dy + Ay dt/\dx/\dz - Ax dt/\dy/\dz - At dx/\dy/\dz
            sage: ssa = sa.hodge_star(g) ; ssa
            1-form '**A' on the 4-dimensional manifold 'M'
            sage: ssa.view()
            **A = At dt + Ax dx + Ay dy + Az dz
            sage: ssa == a  # must hold for a Lorentzian metric in dimension 4
            True

        Hodge star of a 2-form in Minkowksi spacetime::
        
            sage: F = DiffForm(m, 2, 'F')    
            sage: var('Ex Ey Ez Bx By Bz')
            (Ex, Ey, Ez, Bx, By, Bz)
            sage: F[0,1], F[0,2], F[0,3] = -Ex, -Ey, -Ez
            sage: F[1,2], F[1,3], F[2,3] = Bz, -By, Bx
            sage: F[:]
            [  0 -Ex -Ey -Ez]
            [ Ex   0  Bz -By]
            [ Ey -Bz   0  Bx]
            [ Ez  By -Bx   0]
            sage: sF = F.hodge_star(g) ; sF
            2-form '*F' on the 4-dimensional manifold 'M'
            sage: sF[:]
            [  0  Bx  By  Bz]
            [-Bx   0  Ez -Ey]
            [-By -Ez   0  Ex]
            [-Bz  Ey -Ex   0]
            sage: ssF = sF.hodge_star(g) ; ssF
            2-form '**F' on the 4-dimensional manifold 'M'
            sage: ssF[:]   
            [  0  Ex  Ey  Ez]
            [-Ex   0 -Bz  By]
            [-Ey  Bz   0 -Bx]
            [-Ez -By  Bx   0]
            sage: ssF.view()
            **F = Ex dt/\dx + Ey dt/\dy + Ez dt/\dz - Bz dx/\dy + By dx/\dz - Bx dy/\dz
            sage: F.view()
            F = -Ex dt/\dx - Ey dt/\dy - Ez dt/\dz + Bz dx/\dy - By dx/\dz + Bx dy/\dz
            sage: ssF == -F  # must hold for a Lorentzian metric in dimension 4
            True

        Test of the standard identity
        
        .. MATH::
            
            *(A\wedge B) = \epsilon(A^\sharp, B^\sharp, ., .)
            
        where `A` and `B` are any 1-forms and `A^\sharp` and `B^\sharp` the 
        vectors associated to them by the metric `g` (index raising)::

            sage: b = OneForm(m, 'B')
            sage: var('Bt Bx By Bz')
            (Bt, Bx, By, Bz)
            sage: b[:] = (Bt, Bx, By, Bz) ; b.view()
            B = Bt dt + Bx dx + By dy + Bz dz
            sage: epsilon = g.volume_form()
            sage: (a.wedge(b)).hodge_star(g) == epsilon.contract(0, a.up(g), 0).contract(0, b.up(g), 0)
            True

        """
        from sage.functions.other import factorial
        from utilities import format_unop_txt, format_unop_latex
        p = self.rank
        eps = metric.volume_form(p)
        if p == 0:
            resu = self * eps
        else:
            resu = self.contract(0, eps, 0)
            for j in range(1, p):
                resu = resu.self_contract(0, p-j)
            if p > 1:
                resu = resu / factorial(p)
        # Name and LaTeX name of the result:
        resu.name = format_unop_txt('*', self.name)
        resu.latex_name = format_unop_latex(r'\star ', self.latex_name)
        return resu