Revision history [back]

One can employ matrix square root to significantly speed up things. The following implementation illustrates this idea and supports only matrices with distinct eigenvalues:

def sqrt_pm1(M):
    n = M.nrows()

    d, L = (M/3).eigenmatrix_right()
    d = d.diagonal()
    assert len(set(d)) == n          # only this case is supported
    L = L.change_ring(QQbar)  

    pm1 = {QQbar(-1),QQbar(1)}
    return [R for signs in Tuples((-1,1),n) if set((R := L * diagonal_matrix([sqrt(e)*s for e,s in zip(d,signs)]) * L.inverse()).list()).issubset(pm1)]

For example,

sqrt_pm1(Matrix([[0,-6],[6,0]]))

gives the following list of two $\pm 1$-matrices that are square roots of Matrix([[0,-6],[6,0]]):

[
[-1  1]  [ 1 -1]
[-1 -1], [ 1  1]
]

One can employ matrix square root to significantly speed up things. The following implementation illustrates this idea and supports only matrices with distinct eigenvalues:

def sqrt_pm1(M):
    n = M.nrows()

    d, L = (M/3).eigenmatrix_right()
    d = d.diagonal()
    assert len(set(d)) == n          # only this case is supported
    L = L.change_ring(QQbar)  

    pm1 = {QQbar(-1),QQbar(1)}
    return [R for signs in Tuples((-1,1),n) if set((R := L * diagonal_matrix([sqrt(e)*s for e,s in zip(d,signs)]) * L.inverse()).list()).issubset(pm1)]

For example,

sqrt_pm1(Matrix([[0,-6],[6,0]]))

gives the following list of two $\pm 1$-matrices that are square roots of Matrix([[0,-6],[6,0]])Matrix([[0,-6],[6,0]])/3:

[
[-1  1]  [ 1 -1]
[-1 -1], [ 1  1]
]

One can employ matrix square root to significantly speed up things. The following implementation illustrates this idea and supports only matrices with distinct eigenvalues:

def sqrt_pm1(M):
    if not M.det().is_square():    # quick test for insolubility
        return []

    n = M.nrows()

    d, L = (M/3).eigenmatrix_right()
M.eigenmatrix_right()
    d = d.diagonal()
    assert len(set(d)) == n          # only this case is supported
    L = L.change_ring(QQbar)  

    pm1 = {QQbar(-1),QQbar(1)}
    return [R for signs in Tuples((-1,1),n) if set((R := L * diagonal_matrix([sqrt(e)*s for e,s in zip(d,signs)]) * L.inverse()).list()).issubset(pm1)]

For example,

sqrt_pm1(Matrix([[0,-6],[6,0]]))
sqrt_pm1(Matrix([[0,-6],[6,0]]) / 3)

gives the following list of two $\pm 1$-matrices that are square roots of Matrix([[0,-6],[6,0]])/3Matrix([[0,-6],[6,0]]) / 3:

[
[-1  1]  [ 1 -1]
[-1 -1], [ 1  1]
]