Vehicle-cpp/hkpc/software/sympy/matrices/tests/test_eigen.py

from sympy.core.evalf import N
from sympy.core.numbers import (Float, I, Rational)
from sympy.core.symbol import (Symbol, symbols)
from sympy.functions.elementary.complexes import Abs
from sympy.functions.elementary.exponential import exp
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.trigonometric import (cos, sin)
from sympy.matrices import eye, Matrix
from sympy.core.singleton import S
from sympy.testing.pytest import raises, XFAIL
from sympy.matrices.matrices import NonSquareMatrixError, MatrixError
from sympy.matrices.expressions.fourier import DFT
from sympy.simplify.simplify import simplify
from sympy.matrices.immutable import ImmutableMatrix
from sympy.testing.pytest import slow
from sympy.testing.matrices import allclose


def test_eigen():
    R = Rational
    M = Matrix.eye(3)
    assert M.eigenvals(multiple=False) == {S.One: 3}
    assert M.eigenvals(multiple=True) == [1, 1, 1]

    assert M.eigenvects() == (
        [(1, 3, [Matrix([1, 0, 0]),
                 Matrix([0, 1, 0]),
                 Matrix([0, 0, 1])])])

    assert M.left_eigenvects() == (
        [(1, 3, [Matrix([[1, 0, 0]]),
                 Matrix([[0, 1, 0]]),
                 Matrix([[0, 0, 1]])])])

    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])

    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    assert M.eigenvects() == (
        [
            (-1, 1, [Matrix([-1, 1, 0])]),
            ( 0, 1, [Matrix([0, -1, 1])]),
            ( 2, 1, [Matrix([R(2, 3), R(1, 3), 1])])
        ])

    assert M.left_eigenvects() == (
        [
            (-1, 1, [Matrix([[-2, 1, 1]])]),
            (0, 1, [Matrix([[-1, -1, 1]])]),
            (2, 1, [Matrix([[1, 1, 1]])])
        ])

    a = Symbol('a')
    M = Matrix([[a, 0],
                [0, 1]])

    assert M.eigenvals() == {a: 1, S.One: 1}

    M = Matrix([[1, -1],
                [1,  3]])
    assert M.eigenvects() == ([(2, 2, [Matrix(2, 1, [-1, 1])])])
    assert M.left_eigenvects() == ([(2, 2, [Matrix([[1, 1]])])])

    M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    a = R(15, 2)
    b = 3*33**R(1, 2)
    c = R(13, 2)
    d = (R(33, 8) + 3*b/8)
    e = (R(33, 8) - 3*b/8)

    def NS(e, n):
        return str(N(e, n))
    r = [
        (a - b/2, 1, [Matrix([(12 + 24/(c - b/2))/((c - b/2)*e) + 3/(c - b/2),
                              (6 + 12/(c - b/2))/e, 1])]),
        (      0, 1, [Matrix([1, -2, 1])]),
        (a + b/2, 1, [Matrix([(12 + 24/(c + b/2))/((c + b/2)*d) + 3/(c + b/2),
                              (6 + 12/(c + b/2))/d, 1])]),
    ]
    r1 = [(NS(r[i][0], 2), NS(r[i][1], 2),
        [NS(j, 2) for j in r[i][2][0]]) for i in range(len(r))]
    r = M.eigenvects()
    r2 = [(NS(r[i][0], 2), NS(r[i][1], 2),
        [NS(j, 2) for j in r[i][2][0]]) for i in range(len(r))]
    assert sorted(r1) == sorted(r2)

    eps = Symbol('eps', real=True)

    M = Matrix([[abs(eps), I*eps    ],
                [-I*eps,   abs(eps) ]])

    assert M.eigenvects() == (
        [
            ( 0, 1, [Matrix([[-I*eps/abs(eps)], [1]])]),
            ( 2*abs(eps), 1, [ Matrix([[I*eps/abs(eps)], [1]]) ] ),
        ])

    assert M.left_eigenvects() == (
        [
            (0, 1, [Matrix([[I*eps/Abs(eps), 1]])]),
            (2*Abs(eps), 1, [Matrix([[-I*eps/Abs(eps), 1]])])
        ])

    M = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
    M._eigenvects = M.eigenvects(simplify=False)
    assert max(i.q for i in M._eigenvects[0][2][0]) > 1
    M._eigenvects = M.eigenvects(simplify=True)
    assert max(i.q for i in M._eigenvects[0][2][0]) == 1

    M = Matrix([[Rational(1, 4), 1], [1, 1]])
    assert M.eigenvects() == [
        (Rational(5, 8) - sqrt(73)/8, 1, [Matrix([[-sqrt(73)/8 - Rational(3, 8)], [1]])]),
        (Rational(5, 8) + sqrt(73)/8, 1, [Matrix([[Rational(-3, 8) + sqrt(73)/8], [1]])])]

    # issue 10719
    assert Matrix([]).eigenvals() == {}
    assert Matrix([]).eigenvals(multiple=True) == []
    assert Matrix([]).eigenvects() == []

    # issue 15119
    raises(NonSquareMatrixError,
           lambda: Matrix([[1, 2], [0, 4], [0, 0]]).eigenvals())
    raises(NonSquareMatrixError,
           lambda: Matrix([[1, 0], [3, 4], [5, 6]]).eigenvals())
    raises(NonSquareMatrixError,
           lambda: Matrix([[1, 2, 3], [0, 5, 6]]).eigenvals())
    raises(NonSquareMatrixError,
           lambda: Matrix([[1, 0, 0], [4, 5, 0]]).eigenvals())
    raises(NonSquareMatrixError,
           lambda: Matrix([[1, 2, 3], [0, 5, 6]]).eigenvals(
               error_when_incomplete = False))
    raises(NonSquareMatrixError,
           lambda: Matrix([[1, 0, 0], [4, 5, 0]]).eigenvals(
               error_when_incomplete = False))

    m = Matrix([[1, 2], [3, 4]])
    assert isinstance(m.eigenvals(simplify=True, multiple=False), dict)
    assert isinstance(m.eigenvals(simplify=True, multiple=True), list)
    assert isinstance(m.eigenvals(simplify=lambda x: x, multiple=False), dict)
    assert isinstance(m.eigenvals(simplify=lambda x: x, multiple=True), list)


@slow
def test_eigen_slow():
    # issue 15125
    from sympy.core.function import count_ops
    q = Symbol("q", positive = True)
    m = Matrix([[-2, exp(-q), 1], [exp(q), -2, 1], [1, 1, -2]])
    assert count_ops(m.eigenvals(simplify=False)) > \
        count_ops(m.eigenvals(simplify=True))
    assert count_ops(m.eigenvals(simplify=lambda x: x)) > \
        count_ops(m.eigenvals(simplify=True))


def test_float_eigenvals():
    m = Matrix([[1, .6, .6], [.6, .9, .9], [.9, .6, .6]])
    evals = [
        Rational(5, 4) - sqrt(385)/20,
        sqrt(385)/20 + Rational(5, 4),
        S.Zero]

    n_evals = m.eigenvals(rational=True, multiple=True)
    n_evals = sorted(n_evals)
    s_evals = [x.evalf() for x in evals]
    s_evals = sorted(s_evals)

    for x, y in zip(n_evals, s_evals):
        assert abs(x-y) < 10**-9


@XFAIL
def test_eigen_vects():
    m = Matrix(2, 2, [1, 0, 0, I])
    raises(NotImplementedError, lambda: m.is_diagonalizable(True))
    # !!! bug because of eigenvects() or roots(x**2 + (-1 - I)*x + I, x)
    # see issue 5292
    assert not m.is_diagonalizable(True)
    raises(MatrixError, lambda: m.diagonalize(True))
    (P, D) = m.diagonalize(True)

def test_issue_8240():
    # Eigenvalues of large triangular matrices
    x, y = symbols('x y')
    n = 200

    diagonal_variables = [Symbol('x%s' % i) for i in range(n)]
    M = [[0 for i in range(n)] for j in range(n)]
    for i in range(n):
        M[i][i] = diagonal_variables[i]
    M = Matrix(M)

    eigenvals = M.eigenvals()
    assert len(eigenvals) == n
    for i in range(n):
        assert eigenvals[diagonal_variables[i]] == 1

    eigenvals = M.eigenvals(multiple=True)
    assert set(eigenvals) == set(diagonal_variables)

    # with multiplicity
    M = Matrix([[x, 0, 0], [1, y, 0], [2, 3, x]])
    eigenvals = M.eigenvals()
    assert eigenvals == {x: 2, y: 1}

    eigenvals = M.eigenvals(multiple=True)
    assert len(eigenvals) == 3
    assert eigenvals.count(x) == 2
    assert eigenvals.count(y) == 1


def test_eigenvals():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    m = Matrix([
        [3,  0,  0, 0, -3],
        [0, -3, -3, 0,  3],
        [0,  3,  0, 3,  0],
        [0,  0,  3, 0,  3],
        [3,  0,  0, 3,  0]])

    # XXX Used dry-run test because arbitrary symbol that appears in
    # CRootOf may not be unique.
    assert m.eigenvals()


def test_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert M*vec_list[0] == val*vec_list[0]


def test_left_eigenvects():
    M = Matrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    vecs = M.left_eigenvects()
    for val, mult, vec_list in vecs:
        assert len(vec_list) == 1
        assert vec_list[0]*M == val*vec_list[0]


@slow
def test_bidiagonalize():
    M = Matrix([[1, 0, 0],
                [0, 1, 0],
                [0, 0, 1]])
    assert M.bidiagonalize() == M
    assert M.bidiagonalize(upper=False) == M
    assert M.bidiagonalize() == M
    assert M.bidiagonal_decomposition() == (M, M, M)
    assert M.bidiagonal_decomposition(upper=False) == (M, M, M)
    assert M.bidiagonalize() == M

    import random
    #Real Tests
    for real_test in range(2):
        test_values = []
        row = 2
        col = 2
        for _ in range(row * col):
            value = random.randint(-1000000000, 1000000000)
            test_values = test_values + [value]
        # L     -> Lower Bidiagonalization
        # M     -> Mutable Matrix
        # N     -> Immutable Matrix
        # 0     -> Bidiagonalized form
        # 1,2,3 -> Bidiagonal_decomposition matrices
        # 4     -> Product of 1 2 3
        M = Matrix(row, col, test_values)
        N = ImmutableMatrix(M)

        N1, N2, N3 = N.bidiagonal_decomposition()
        M1, M2, M3 = M.bidiagonal_decomposition()
        M0 = M.bidiagonalize()
        N0 = N.bidiagonalize()

        N4 = N1 * N2 * N3
        M4 = M1 * M2 * M3

        N2.simplify()
        N4.simplify()
        N0.simplify()

        M0.simplify()
        M2.simplify()
        M4.simplify()

        LM0 = M.bidiagonalize(upper=False)
        LM1, LM2, LM3 = M.bidiagonal_decomposition(upper=False)
        LN0 = N.bidiagonalize(upper=False)
        LN1, LN2, LN3 = N.bidiagonal_decomposition(upper=False)

        LN4 = LN1 * LN2 * LN3
        LM4 = LM1 * LM2 * LM3

        LN2.simplify()
        LN4.simplify()
        LN0.simplify()

        LM0.simplify()
        LM2.simplify()
        LM4.simplify()

        assert M == M4
        assert M2 == M0
        assert N == N4
        assert N2 == N0
        assert M == LM4
        assert LM2 == LM0
        assert N == LN4
        assert LN2 == LN0

    #Complex Tests
    for complex_test in range(2):
        test_values = []
        size = 2
        for _ in range(size * size):
            real = random.randint(-1000000000, 1000000000)
            comp = random.randint(-1000000000, 1000000000)
            value = real + comp * I
            test_values = test_values + [value]
        M = Matrix(size, size, test_values)
        N = ImmutableMatrix(M)
        # L     -> Lower Bidiagonalization
        # M     -> Mutable Matrix
        # N     -> Immutable Matrix
        # 0     -> Bidiagonalized form
        # 1,2,3 -> Bidiagonal_decomposition matrices
        # 4     -> Product of 1 2 3
        N1, N2, N3 = N.bidiagonal_decomposition()
        M1, M2, M3 = M.bidiagonal_decomposition()
        M0 = M.bidiagonalize()
        N0 = N.bidiagonalize()

        N4 = N1 * N2 * N3
        M4 = M1 * M2 * M3

        N2.simplify()
        N4.simplify()
        N0.simplify()

        M0.simplify()
        M2.simplify()
        M4.simplify()

        LM0 = M.bidiagonalize(upper=False)
        LM1, LM2, LM3 = M.bidiagonal_decomposition(upper=False)
        LN0 = N.bidiagonalize(upper=False)
        LN1, LN2, LN3 = N.bidiagonal_decomposition(upper=False)

        LN4 = LN1 * LN2 * LN3
        LM4 = LM1 * LM2 * LM3

        LN2.simplify()
        LN4.simplify()
        LN0.simplify()

        LM0.simplify()
        LM2.simplify()
        LM4.simplify()

        assert M == M4
        assert M2 == M0
        assert N == N4
        assert N2 == N0
        assert M == LM4
        assert LM2 == LM0
        assert N == LN4
        assert LN2 == LN0

    M = Matrix(18, 8, range(1, 145))
    M = M.applyfunc(lambda i: Float(i))
    assert M.bidiagonal_decomposition()[1] == M.bidiagonalize()
    assert M.bidiagonal_decomposition(upper=False)[1] == M.bidiagonalize(upper=False)
    a, b, c = M.bidiagonal_decomposition()
    diff = a * b * c - M
    assert abs(max(diff)) < 10**-12


def test_diagonalize():
    m = Matrix(2, 2, [0, -1, 1, 0])
    raises(MatrixError, lambda: m.diagonalize(reals_only=True))
    P, D = m.diagonalize()
    assert D.is_diagonal()
    assert D == Matrix([
                 [-I, 0],
                 [ 0, I]])

    # make sure we use floats out if floats are passed in
    m = Matrix(2, 2, [0, .5, .5, 0])
    P, D = m.diagonalize()
    assert all(isinstance(e, Float) for e in D.values())
    assert all(isinstance(e, Float) for e in P.values())

    _, D2 = m.diagonalize(reals_only=True)
    assert D == D2

    m = Matrix(
        [[0, 1, 0, 0], [1, 0, 0, 0.002], [0.002, 0, 0, 1], [0, 0, 1, 0]])
    P, D = m.diagonalize()
    assert allclose(P*D, m*P)


def test_is_diagonalizable():
    a, b, c = symbols('a b c')
    m = Matrix(2, 2, [a, c, c, b])
    assert m.is_symmetric()
    assert m.is_diagonalizable()
    assert not Matrix(2, 2, [1, 1, 0, 1]).is_diagonalizable()

    m = Matrix(2, 2, [0, -1, 1, 0])
    assert m.is_diagonalizable()
    assert not m.is_diagonalizable(reals_only=True)


def test_jordan_form():
    m = Matrix(3, 2, [-3, 1, -3, 20, 3, 10])
    raises(NonSquareMatrixError, lambda: m.jordan_form())

    # the next two tests test the cases where the old
    # algorithm failed due to the fact that the block structure can
    # *NOT* be determined  from algebraic and geometric multiplicity alone
    # This can be seen most easily when one lets compute the J.c.f. of a matrix that
    # is in J.c.f already.
    m = Matrix(4, 4, [2, 1, 0, 0,
                    0, 2, 1, 0,
                    0, 0, 2, 0,
                    0, 0, 0, 2
    ])
    P, J = m.jordan_form()
    assert m == J

    m = Matrix(4, 4, [2, 1, 0, 0,
                    0, 2, 0, 0,
                    0, 0, 2, 1,
                    0, 0, 0, 2
    ])
    P, J = m.jordan_form()
    assert m == J

    A = Matrix([[ 2,  4,  1,  0],
                [-4,  2,  0,  1],
                [ 0,  0,  2,  4],
                [ 0,  0, -4,  2]])
    P, J = A.jordan_form()
    assert simplify(P*J*P.inv()) == A

    assert Matrix(1, 1, [1]).jordan_form() == (Matrix([1]), Matrix([1]))
    assert Matrix(1, 1, [1]).jordan_form(calc_transform=False) == Matrix([1])

    # If we have eigenvalues in CRootOf form, raise errors
    m = Matrix([[3, 0, 0, 0, -3], [0, -3, -3, 0, 3], [0, 3, 0, 3, 0], [0, 0, 3, 0, 3], [3, 0, 0, 3, 0]])
    raises(MatrixError, lambda: m.jordan_form())

    # make sure that if the input has floats, the output does too
    m = Matrix([
        [                0.6875, 0.125 + 0.1875*sqrt(3)],
        [0.125 + 0.1875*sqrt(3),                 0.3125]])
    P, J = m.jordan_form()
    assert all(isinstance(x, Float) or x == 0 for x in P)
    assert all(isinstance(x, Float) or x == 0 for x in J)


def test_singular_values():
    x = Symbol('x', real=True)

    A = Matrix([[0, 1*I], [2, 0]])
    # if singular values can be sorted, they should be in decreasing order
    assert A.singular_values() == [2, 1]

    A = eye(3)
    A[1, 1] = x
    A[2, 2] = 5
    vals = A.singular_values()
    # since Abs(x) cannot be sorted, test set equality
    assert set(vals) == {5, 1, Abs(x)}

    A = Matrix([[sin(x), cos(x)], [-cos(x), sin(x)]])
    vals = [sv.trigsimp() for sv in A.singular_values()]
    assert vals == [S.One, S.One]

    A = Matrix([
        [2, 4],
        [1, 3],
        [0, 0],
        [0, 0]
        ])
    assert A.singular_values() == \
        [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221))]
    assert A.T.singular_values() == \
        [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221)), 0, 0]

def test___eq__():
    assert (Matrix(
        [[0, 1, 1],
        [1, 0, 0],
        [1, 1, 1]]) == {}) is False


def test_definite():
    # Examples from Gilbert Strang, "Introduction to Linear Algebra"
    # Positive definite matrices
    m = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
    assert m.is_positive_definite == True
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    m = Matrix([[5, 4], [4, 5]])
    assert m.is_positive_definite == True
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    # Positive semidefinite matrices
    m = Matrix([[2, -1, -1], [-1, 2, -1], [-1, -1, 2]])
    assert m.is_positive_definite == False
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    m = Matrix([[1, 2], [2, 4]])
    assert m.is_positive_definite == False
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    # Examples from Mathematica documentation
    # Non-hermitian positive definite matrices
    m = Matrix([[2, 3], [4, 8]])
    assert m.is_positive_definite == True
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    # Hermetian matrices
    m = Matrix([[1, 2*I], [-I, 4]])
    assert m.is_positive_definite == True
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    # Symbolic matrices examples
    a = Symbol('a', positive=True)
    b = Symbol('b', negative=True)
    m = Matrix([[a, 0, 0], [0, a, 0], [0, 0, a]])
    assert m.is_positive_definite == True
    assert m.is_positive_semidefinite == True
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == False

    m = Matrix([[b, 0, 0], [0, b, 0], [0, 0, b]])
    assert m.is_positive_definite == False
    assert m.is_positive_semidefinite == False
    assert m.is_negative_definite == True
    assert m.is_negative_semidefinite == True
    assert m.is_indefinite == False

    m = Matrix([[a, 0], [0, b]])
    assert m.is_positive_definite == False
    assert m.is_positive_semidefinite == False
    assert m.is_negative_definite == False
    assert m.is_negative_semidefinite == False
    assert m.is_indefinite == True

    m = Matrix([
        [0.0228202735623867, 0.00518748979085398,
         -0.0743036351048907, -0.00709135324903921],
        [0.00518748979085398, 0.0349045359786350,
         0.0830317991056637, 0.00233147902806909],
        [-0.0743036351048907, 0.0830317991056637,
         1.15859676366277, 0.340359081555988],
        [-0.00709135324903921, 0.00233147902806909,
         0.340359081555988, 0.928147644848199]
    ])
    assert m.is_positive_definite == True
    assert m.is_positive_semidefinite == True
    assert m.is_indefinite == False

    # test for issue 19547: https://github.com/sympy/sympy/issues/19547
    m = Matrix([
        [0, 0, 0],
        [0, 1, 2],
        [0, 2, 1]
    ])
    assert not m.is_positive_definite
    assert not m.is_positive_semidefinite


def test_positive_semidefinite_cholesky():
    from sympy.matrices.eigen import _is_positive_semidefinite_cholesky

    m = Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
    assert _is_positive_semidefinite_cholesky(m) == True
    m = Matrix([[0, 0, 0], [0, 5, -10*I], [0, 10*I, 5]])
    assert _is_positive_semidefinite_cholesky(m) == False
    m = Matrix([[1, 0, 0], [0, 0, 0], [0, 0, -1]])
    assert _is_positive_semidefinite_cholesky(m) == False
    m = Matrix([[0, 1], [1, 0]])
    assert _is_positive_semidefinite_cholesky(m) == False

    # https://www.value-at-risk.net/cholesky-factorization/
    m = Matrix([[4, -2, -6], [-2, 10, 9], [-6, 9, 14]])
    assert _is_positive_semidefinite_cholesky(m) == True
    m = Matrix([[9, -3, 3], [-3, 2, 1], [3, 1, 6]])
    assert _is_positive_semidefinite_cholesky(m) == True
    m = Matrix([[4, -2, 2], [-2, 1, -1], [2, -1, 5]])
    assert _is_positive_semidefinite_cholesky(m) == True
    m = Matrix([[1, 2, -1], [2, 5, 1], [-1, 1, 9]])
    assert _is_positive_semidefinite_cholesky(m) == False


def test_issue_20582():
    A = Matrix([
        [5, -5, -3, 2, -7],
        [-2, -5, 0, 2, 1],
        [-2, -7, -5, -2, -6],
        [7, 10, 3, 9, -2],
        [4, -10, 3, -8, -4]
    ])
    # XXX Used dry-run test because arbitrary symbol that appears in
    # CRootOf may not be unique.
    assert A.eigenvects()

def test_issue_19210():
    t = Symbol('t')
    H = Matrix([[3, 0, 0, 0], [0, 1 , 2, 0], [0, 2, 2, 0], [0, 0, 0, 4]])
    A = (-I * H * t).jordan_form()
    assert A == (Matrix([
                    [0, 1,                  0,                0],
                    [0, 0, -4/(-1 + sqrt(17)), 4/(1 + sqrt(17))],
                    [0, 0,                  1,                1],
                    [1, 0,                  0,                0]]), Matrix([
                    [-4*I*t,      0,                         0,                         0],
                    [     0, -3*I*t,                         0,                         0],
                    [     0,      0, t*(-3*I/2 + sqrt(17)*I/2),                         0],
                    [     0,      0,                         0, t*(-sqrt(17)*I/2 - 3*I/2)]]))


def test_issue_20275():
    # XXX We use complex expansions because complex exponentials are not
    # recognized by polys.domains
    A = DFT(3).as_explicit().expand(complex=True)
    eigenvects = A.eigenvects()
    assert eigenvects[0] == (
        -1, 1,
        [Matrix([[1 - sqrt(3)], [1], [1]])]
    )
    assert eigenvects[1] == (
        1, 1,
        [Matrix([[1 + sqrt(3)], [1], [1]])]
    )
    assert eigenvects[2] == (
        -I, 1,
        [Matrix([[0], [-1], [1]])]
    )

    A = DFT(4).as_explicit().expand(complex=True)
    eigenvects = A.eigenvects()
    assert eigenvects[0] == (
        -1, 1,
        [Matrix([[-1], [1], [1], [1]])]
    )
    assert eigenvects[1] == (
        1, 2,
        [Matrix([[1], [0], [1], [0]]), Matrix([[2], [1], [0], [1]])]
    )
    assert eigenvects[2] == (
        -I, 1,
        [Matrix([[0], [-1], [0], [1]])]
    )

    # XXX We skip test for some parts of eigenvectors which are very
    # complicated and fragile under expression tree changes
    A = DFT(5).as_explicit().expand(complex=True)
    eigenvects = A.eigenvects()
    assert eigenvects[0] == (
        -1, 1,
        [Matrix([[1 - sqrt(5)], [1], [1], [1], [1]])]
    )
    assert eigenvects[1] == (
        1, 2,
        [Matrix([[S(1)/2 + sqrt(5)/2], [0], [1], [1], [0]]),
         Matrix([[S(1)/2 + sqrt(5)/2], [1], [0], [0], [1]])]
    )


def test_issue_20752():
    b = symbols('b', nonzero=True)
    m = Matrix([[0, 0, 0], [0, b, 0], [0, 0, b]])
    assert m.is_positive_semidefinite is None