from sage.all import *

def is_lyapunov_like(L,K):
    r"""
    Determine whether or not ``L`` is Lyapunov-like on ``K``.

    We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
    L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
    `\left\langle x,s \right\rangle` in the complementarity set of
    ``K``. It is known [Orlitzky]_ that this property need only be
    checked for generators of ``K`` and its dual.

    There are faster ways of checking this property. For example, we
    could compute a `lyapunov_like_basis` of the cone, and then test
    whether or not the given matrix is contained in the span of that
    basis. The value of this function is that it works on symbolic
    matrices.

    INPUT:

    - ``L`` -- A linear transformation or matrix.

    - ``K`` -- A polyhedral closed convex cone.

    OUTPUT:

    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
    and ``False`` otherwise.

    .. WARNING::

        If this function returns ``True``, then ``L`` is Lyapunov-like
        on ``K``. However, if ``False`` is returned, that could mean one
        of two things. The first is that ``L`` is definitely not
        Lyapunov-like on ``K``. The second is more of an "I don't know"
        answer, returned (for example) if we cannot prove that an inner
        product is zero.

    REFERENCES:

    M. Orlitzky. The Lyapunov rank of an improper cone.
    http://www.optimization-online.org/DB_HTML/2015/10/5135.html

    EXAMPLES:

    The identity is always Lyapunov-like in a nontrivial space::

        sage: set_random_seed()
        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
        sage: L = identity_matrix(K.lattice_dim())
        sage: is_lyapunov_like(L,K)
        True

    As is the "zero" transformation::

        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
        sage: R = K.lattice().vector_space().base_ring()
        sage: L = zero_matrix(R, K.lattice_dim())
        sage: is_lyapunov_like(L,K)
        True

        Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
        on ``K``::

        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
        sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
        True

    """
    return all([(L*x).inner_product(s) == 0
                for (x,s) in K.discrete_complementarity_set()])


def positive_operator_gens(K1, K2 = None):
    r"""
    Compute generators of the cone of positive operators on this cone. A
    linear operator on a cone is positive if the image of the cone under
    the operator is a subset of the cone. This concept can be extended
    to two cones, where the image of the first cone under a positive
    operator is a subset of the second cone.

    INPUT:

    - ``K2`` -- (default: ``K1``) the codomain cone; the image of this
                cone under the returned operators is a subset of ``K2``.

    OUTPUT:

    A list of `m`-by-``n`` matrices where ``m == K2.lattice_dim()`` and
    ``n == K1.lattice_dim()``. Each matrix ``P`` in the list should have
    the property that ``P*x`` is an element of ``K2`` whenever ``x`` is
    an element of ``K1``. Moreover, any nonnegative linear combination of
    these matrices shares the same property.

    REFERENCES:

    .. [Orlitzky-Pi-Z]
       M. Orlitzky.
       Positive and Z-operators on closed convex cones.

    .. [Tam]
       B.-S. Tam.
       Some results of polyhedral cones and simplicial cones.
       Linear and Multilinear Algebra, 4:4 (1977) 281--284.

    EXAMPLES:

    Positive operators on the nonnegative orthant are nonnegative matrices::

        sage: K = Cone([(1,)])
        sage: positive_operator_gens(K)
        [[1]]

        sage: K = Cone([(1,0),(0,1)])
        sage: positive_operator_gens(K)
        [
        [1 0]  [0 1]  [0 0]  [0 0]
        [0 0], [0 0], [1 0], [0 1]
        ]

    The trivial cone in a trivial space has no positive operators::

        sage: K = Cone([], ToricLattice(0))
        sage: positive_operator_gens(K)
        []

    Every operator is positive on the trivial cone::

        sage: K = Cone([(0,)])
        sage: positive_operator_gens(K)
        [[1], [-1]]

        sage: K = Cone([(0,0)])
        sage: K.is_trivial()
        True
        sage: positive_operator_gens(K)
        [
        [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
        [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
        ]

    Every operator is positive on the ambient vector space::

        sage: K = Cone([(1,),(-1,)])
        sage: K.is_full_space()
        True
        sage: positive_operator_gens(K)
        [[1], [-1]]

        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
        sage: K.is_full_space()
        True
        sage: positive_operator_gens(K)
        [
        [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
        [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
        ]

    A non-obvious application is to find the positive operators on the
    right half-plane::

        sage: K = Cone([(1,0),(0,1),(0,-1)])
        sage: positive_operator_gens(K)
        [
        [1 0]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
        [0 0], [1 0], [-1  0], [0 1], [ 0 -1]
        ]

    TESTS:

    Each positive operator generator should send the generators of one
    cone into the other cone::

        sage: set_random_seed()
        sage: K1 = random_cone(max_ambient_dim=4)
        sage: K2 = random_cone(max_ambient_dim=4)
        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
        sage: all([ K2.contains(P*x) for P in pi_K1_K2 for x in K1 ])
        True

    Each positive operator generator should send a random element of one
    cone into the other cone::

        sage: set_random_seed()
        sage: K1 = random_cone(max_ambient_dim=4)
        sage: K2 = random_cone(max_ambient_dim=4)
        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
        sage: all([ K2.contains(P*K1.random_element(QQ)) for P in pi_K1_K2 ])
        True

    A random element of the positive operator cone should send the
    generators of one cone into the other cone::

        sage: set_random_seed()
        sage: K1 = random_cone(max_ambient_dim=4)
        sage: K2 = random_cone(max_ambient_dim=4)
        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
        sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim())
        sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ],
        ....:                lattice=L,
        ....:                check=False)
        sage: P = matrix(K2.lattice_dim(),
        ....:            K1.lattice_dim(),
        ....:            pi_cone.random_element(QQ).list())
        sage: all([ K2.contains(P*x) for x in K1 ])
        True

    A random element of the positive operator cone should send a random
    element of one cone into the other cone::

        sage: set_random_seed()
        sage: K1 = random_cone(max_ambient_dim=4)
        sage: K2 = random_cone(max_ambient_dim=4)
        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
        sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim())
        sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ],
        ....:                lattice=L,
        ....:                check=False)
        sage: P = matrix(K2.lattice_dim(),
        ....:            K1.lattice_dim(),
        ....:            pi_cone.random_element(QQ).list())
        sage: K2.contains(P*K1.random_element(ring=QQ))
        True

    The lineality space of the dual of the cone of positive operators
    can be computed from the lineality spaces of the cone and its dual::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: pi_of_K = positive_operator_gens(K)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = pi_cone.dual().linear_subspace()
        sage: U1 = [ vector((s.tensor_product(x)).list())
        ....:        for x in K.lines()
        ....:        for s in K.dual() ]
        sage: U2 = [ vector((s.tensor_product(x)).list())
        ....:        for x in K
        ....:        for s in K.dual().lines() ]
        sage: expected = pi_cone.lattice().vector_space().span(U1 + U2)
        sage: actual == expected
        True

    The lineality of the dual of the cone of positive operators
    is known from its lineality space::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: n = K.lattice_dim()
        sage: m = K.dim()
        sage: l = K.lineality()
        sage: pi_of_K = positive_operator_gens(K)
        sage: L = ToricLattice(n**2)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                 lattice=L,
        ....:                 check=False)
        sage: actual = pi_cone.dual().lineality()
        sage: expected = l*(m - l) + m*(n - m)
        sage: actual == expected
        True

    The dimension of the cone of positive operators is given by the
    corollary in my paper::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: n = K.lattice_dim()
        sage: m = K.dim()
        sage: l = K.lineality()
        sage: pi_of_K = positive_operator_gens(K)
        sage: L = ToricLattice(n**2)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = pi_cone.dim()
        sage: expected = n**2 - l*(m - l) - (n - m)*m
        sage: actual == expected
        True

    The trivial cone, full space, and half-plane all give rise to the
    expected dimensions::

        sage: n = ZZ.random_element().abs()
        sage: K = Cone([[0] * n], ToricLattice(n))
        sage: K.is_trivial()
        True
        sage: L = ToricLattice(n^2)
        sage: pi_of_K = positive_operator_gens(K)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = pi_cone.dim()
        sage: actual == n^2
        True
        sage: K = K.dual()
        sage: K.is_full_space()
        True
        sage: pi_of_K = positive_operator_gens(K)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = pi_cone.dim()
        sage: actual == n^2
        True
        sage: K = Cone([(1,0),(0,1),(0,-1)])
        sage: pi_of_K = positive_operator_gens(K)
        sage: actual = Cone([p.list() for p in pi_of_K], check=False).dim()
        sage: actual == 3
        True

    The lineality of the cone of positive operators follows from the
    description of its generators::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: n = K.lattice_dim()
        sage: pi_of_K = positive_operator_gens(K)
        sage: L = ToricLattice(n**2)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = pi_cone.lineality()
        sage: expected = n**2 - K.dim()*K.dual().dim()
        sage: actual == expected
        True

    The trivial cone, full space, and half-plane all give rise to the
    expected linealities::

        sage: n = ZZ.random_element().abs()
        sage: K = Cone([[0] * n], ToricLattice(n))
        sage: K.is_trivial()
        True
        sage: L = ToricLattice(n^2)
        sage: pi_of_K = positive_operator_gens(K)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = pi_cone.lineality()
        sage: actual == n^2
        True
        sage: K = K.dual()
        sage: K.is_full_space()
        True
        sage: pi_of_K = positive_operator_gens(K)
        sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
        sage: pi_cone.lineality() == n^2
        True
        sage: K = Cone([(1,0),(0,1),(0,-1)])
        sage: pi_of_K = positive_operator_gens(K)
        sage: pi_cone = Cone([p.list() for p in pi_of_K], check=False)
        sage: actual = pi_cone.lineality()
        sage: actual == 2
        True

    A cone is proper if and only if its cone of positive operators
    is proper::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: pi_of_K = positive_operator_gens(K)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: K.is_proper() == pi_cone.is_proper()
        True

    The positive operators of a permuted cone can be obtained by
    conjugation::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
        sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
        sage: pi_of_pK = positive_operator_gens(pK)
        sage: actual = Cone([t.list() for t in pi_of_pK],
        ....:                lattice=L,
        ....:                check=False)
        sage: pi_of_K = positive_operator_gens(K)
        sage: expected = Cone([(p*t*p.inverse()).list() for t in pi_of_K],
        ....:                   lattice=L,
        ....:                   check=False)
        sage: actual.is_equivalent(expected)
        True

    A transformation is positive on a cone if and only if its adjoint is
    positive on the dual of that cone::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: F = K.lattice().vector_space().base_field()
        sage: n = K.lattice_dim()
        sage: L = ToricLattice(n**2)
        sage: W = VectorSpace(F, n**2)
        sage: pi_of_K = positive_operator_gens(K)
        sage: pi_of_K_star = positive_operator_gens(K.dual())
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: pi_star = Cone([p.list() for p in pi_of_K_star],
        ....:                lattice=L,
        ....:                check=False)
        sage: M = MatrixSpace(F, n)
        sage: L = M(pi_cone.random_element(ring=QQ).list())
        sage: pi_star.contains(W(L.transpose().list()))
        True

        sage: L = W.random_element()
        sage: L_star = W(M(L.list()).transpose().list())
        sage: pi_cone.contains(L) ==  pi_star.contains(L_star)
        True

    The Lyapunov rank of the positive operator cone is the product of
    the Lyapunov ranks of the associated cones if they're all proper::

        sage: K1 = random_cone(max_ambient_dim=4,
        ....:                  strictly_convex=True,
        ....:                  solid=True)
        sage: K2 = random_cone(max_ambient_dim=4,
        ....:                  strictly_convex=True,
        ....:                  solid=True)
        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
        sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim())
        sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ],
        ....:                lattice=L,
        ....:                check=False)
        sage: beta1 = K1.lyapunov_rank()
        sage: beta2 = K2.lyapunov_rank()
        sage: pi_cone.lyapunov_rank() == beta1*beta2
        True

    """
    if K2 is None:
        K2 = K1

    # Matrices are not vectors in Sage, so we have to convert them
    # to vectors explicitly before we can find a basis. We need these
    # two values to construct the appropriate "long vector" space.
    F = K1.lattice().base_field()
    n = K1.lattice_dim()
    m = K2.lattice_dim()

    tensor_products = [ s.tensor_product(x) for x in K1 for s in K2.dual() ]

    # Convert those tensor products to long vectors.
    W = VectorSpace(F, n*m)
    vectors = [ W(tp.list()) for tp in tensor_products ]

    check = True
    if K1.is_proper() and K2.is_proper():
        # All of the generators involved are extreme vectors and
        # therefore minimal [Tam]_. If this cone is neither solid nor
        # strictly convex, then the tensor product of ``s`` and ``x``
        # is the same as that of ``-s`` and ``-x``. However, as a
        # /set/, ``tensor_products`` may still be minimal.
        check = False

    # Create the dual cone of the positive operators, expressed as
    # long vectors.
    pi_dual = Cone(vectors, ToricLattice(W.dimension()), check=check)

    # Now compute the desired cone from its dual...
    pi_cone = pi_dual.dual()

    # And finally convert its rays back to matrix representations.
    M = MatrixSpace(F, m, n)
    return [ M(v.list()) for v in pi_cone ]


def Z_operator_gens(K):
    r"""
    Compute generators of the cone of Z-operators on this cone.

    OUTPUT:

    A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
    Each matrix ``L`` in the list should have the property that
    ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element of
    this cone's :meth:`discrete_complementarity_set`. Moreover, any
    conic (nonnegative linear) combination of these matrices shares the
    same property.

    REFERENCES:

    M. Orlitzky.
    Positive and Z-operators on closed convex cones.

    EXAMPLES:

    Z-operators on the nonnegative orthant are just Z-matrices.
    That is, matrices whose off-diagonal elements are nonnegative::

        sage: K = Cone([(1,0),(0,1)])
        sage: Z_operator_gens(K)
        [
        [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
        [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
        ]
        sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
        sage: all([ z[i][j] <= 0 for z in Z_operator_gens(K)
        ....:                    for i in range(z.nrows())
        ....:                    for j in range(z.ncols())
        ....:                    if i != j ])
        True

    The trivial cone in a trivial space has no Z-operators::

        sage: K = Cone([], ToricLattice(0))
        sage: Z_operator_gens(K)
        []

    Every operator is a Z-operator on the ambient vector space::

        sage: K = Cone([(1,),(-1,)])
        sage: K.is_full_space()
        True
        sage: Z_operator_gens(K)
        [[-1], [1]]

        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
        sage: K.is_full_space()
        True
        sage: Z_operator_gens(K)
        [
        [-1  0]  [1 0]  [ 0 -1]  [0 1]  [ 0  0]  [0 0]  [ 0  0]  [0 0]
        [ 0  0], [0 0], [ 0  0], [0 0], [-1  0], [1 0], [ 0 -1], [0 1]
        ]

    A non-obvious application is to find the Z-operators on the
    right half-plane::

        sage: K = Cone([(1,0),(0,1),(0,-1)])
        sage: Z_operator_gens(K)
        [
        [-1  0]  [1 0]  [ 0  0]  [0 0]  [ 0  0]  [0 0]
        [ 0  0], [0 0], [-1  0], [1 0], [ 0 -1], [0 1]
        ]

    Z-operators on a subspace are Lyapunov-like and vice-versa::

        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
        sage: K.is_full_space()
        True
        sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
        sage: zs  = span([ vector(z.list()) for z in Z_operator_gens(K) ])
        sage: zs == lls
        True

    TESTS:

    The Z-property is possessed by every Z-operator::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: Z_of_K = Z_operator_gens(K)
        sage: dcs = K.discrete_complementarity_set()
        sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
        ....:                                  for (x,s) in dcs])
        True

    The lineality space of the cone of Z-operators is the space of
    Lyapunov-like operators::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: Z_cone = Cone([ z.list() for z in Z_operator_gens(K) ],
        ....:               lattice=L,
        ....:               check=False)
        sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
        sage: lls = L.vector_space().span(ll_basis)
        sage: Z_cone.linear_subspace() == lls
        True

    The lineality of the Z-operators on a cone is the Lyapunov
    rank of that cone::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: Z_of_K = Z_operator_gens(K)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: Z_cone  = Cone([ z.list() for z in Z_of_K ],
        ....:                lattice=L,
        ....:                check=False)
        sage: Z_cone.lineality() == K.lyapunov_rank()
        True

    The lineality spaces of the duals of the positive and Z-operator
    cones are equal. From this it follows that the dimensions of the
    Z-operator cone and positive operator cone are equal::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: pi_of_K = positive_operator_gens(K)
        sage: Z_of_K = Z_operator_gens(K)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: pi_cone = Cone([p.list() for p in pi_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: Z_cone = Cone([ z.list() for z in Z_of_K],
        ....:               lattice=L,
        ....:               check=False)
        sage: pi_cone.dim() == Z_cone.dim()
        True
        sage: pi_star = pi_cone.dual()
        sage: z_star = Z_cone.dual()
        sage: pi_star.linear_subspace() == z_star.linear_subspace()
        True

    The trivial cone, full space, and half-plane all give rise to the
    expected dimensions::

        sage: n = ZZ.random_element().abs()
        sage: K = Cone([[0] * n], ToricLattice(n))
        sage: K.is_trivial()
        True
        sage: L = ToricLattice(n^2)
        sage: Z_of_K = Z_operator_gens(K)
        sage: Z_cone = Cone([z.list() for z in Z_of_K],
        ....:               lattice=L,
        ....:               check=False)
        sage: actual = Z_cone.dim()
        sage: actual == n^2
        True
        sage: K = K.dual()
        sage: K.is_full_space()
        True
        sage: Z_of_K = Z_operator_gens(K)
        sage: Z_cone = Cone([z.list() for z in Z_of_K],
        ....:                lattice=L,
        ....:                check=False)
        sage: actual = Z_cone.dim()
        sage: actual == n^2
        True
        sage: K = Cone([(1,0),(0,1),(0,-1)])
        sage: Z_of_K = Z_operator_gens(K)
        sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False)
        sage: Z_cone.dim() == 3
        True

    The Z-operators of a permuted cone can be obtained by conjugation::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: L = ToricLattice(K.lattice_dim()**2)
        sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
        sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
        sage: Z_of_pK = Z_operator_gens(pK)
        sage: actual = Cone([t.list() for t in Z_of_pK],
        ....:                lattice=L,
        ....:                check=False)
        sage: Z_of_K = Z_operator_gens(K)
        sage: expected = Cone([(p*t*p.inverse()).list() for t in Z_of_K],
        ....:                   lattice=L,
        ....:                   check=False)
        sage: actual.is_equivalent(expected)
        True

    An operator is a Z-operator on a cone if and only if its
    adjoint is a Z-operator on the dual of that cone::

        sage: set_random_seed()
        sage: K = random_cone(max_ambient_dim=4)
        sage: F = K.lattice().vector_space().base_field()
        sage: n = K.lattice_dim()
        sage: L = ToricLattice(n**2)
        sage: W = VectorSpace(F, n**2)
        sage: Z_of_K = Z_operator_gens(K)
        sage: Z_of_K_star = Z_operator_gens(K.dual())
        sage: Z_cone = Cone([p.list() for p in Z_of_K],
        ....:               lattice=L,
        ....:               check=False)
        sage: Z_star = Cone([p.list() for p in Z_of_K_star],
        ....:               lattice=L,
        ....:               check=False)
        sage: M = MatrixSpace(F, n)
        sage: L = M(Z_cone.random_element(ring=QQ).list())
        sage: Z_star.contains(W(L.transpose().list()))
        True

        sage: L = W.random_element()
        sage: L_star = W(M(L.list()).transpose().list())
        sage: Z_cone.contains(L) ==  Z_star.contains(L_star)
        True
    """
    # Matrices are not vectors in Sage, so we have to convert them
    # to vectors explicitly before we can find a basis. We need these
    # two values to construct the appropriate "long vector" space.
    F = K.lattice().base_field()
    n = K.lattice_dim()

    # These tensor products contain generators for the dual cone of
    # the cross-positive operators.
    tensor_products = [ s.tensor_product(x)
                        for (x,s) in K.discrete_complementarity_set() ]

    # Turn our matrices into long vectors...
    W = VectorSpace(F, n**2)
    vectors = [ W(m.list()) for m in tensor_products ]

    check = True
    if K.is_proper():
        # All of the generators involved are extreme vectors and
        # therefore minimal. If this cone is neither solid nor
        # strictly convex, then the tensor product of ``s`` and ``x``
        # is the same as that of ``-s`` and ``-x``. However, as a
        # /set/, ``tensor_products`` may still be minimal.
        check = False

    # Create the dual cone of the cross-positive operators,
    # expressed as long vectors.
    Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()), check=check)

    # Now compute the desired cone from its dual...
    Sigma_cone = Sigma_dual.dual()

    # And finally convert its rays back to matrix representations.
    # But first, make them negative, so we get Z-operators and
    # not cross-positive ones.
    M = MatrixSpace(F, n)
    return [ -M(v.list()) for v in Sigma_cone ]


def LL_cone(K):
    gens = K.lyapunov_like_basis()
    L = ToricLattice(K.lattice_dim()**2)
    return Cone([ g.list() for g in gens ], lattice=L, check=False)

def Z_cone(K):
    gens = Z_operator_gens(K)
    L = ToricLattice(K.lattice_dim()**2)
    return Cone([ g.list() for g in gens ], lattice=L, check=False)

def pi_cone(K):
    gens = positive_operator_gens(K)
    L = ToricLattice(K.lattice_dim()**2)
    return Cone([ g.list() for g in gens ], lattice=L, check=False)