from sympy import *
init_printing(use_unicode=True)

r1,r2,r3,s1,s2,s3 = symbols('r_1 r_2 r_3 s_1 s_2 s_3')
I = Matrix([[1,0],[0,1]]); X = Matrix([[0,1],[1,0]]); Y = Matrix([[0,-1j],[1j,0]]); Z = Matrix([[1,0],[0,-1]])
rho = (1/2)*(I+r1*X+r2*Y+r3*Z); sigma = (1/2)*(I+s1*X+s2*Y+s3*Z)
#rho, sigma

def frho(r1,r2,r3):
    return (1/2)*(I+r1*X+r2*Y+r3*Z)
def fsigma(s1,s2,s3):
    return (1/2)*(I+s1*X+s2*Y+s3*Z)

A = frho(r1,0,r2)*fsigma(0,s2,0)
simplify(A.eigenvals())

A = frho(0,0,r3); B = fsigma(s1,s2,0)
simplify(A*(B**2)*A - B*(A**2)*B)

M = A*B; simplify(M)

simplify(M.eigenvals()) # parecem ser positivos o que esta na raiz e o autovalores









