I have a bunch of matrices with integer coefficients that simultaneously commute. I know that there is a basis that simultaneously diagonalizes all of them, and I want to find it exactly so that I can recover all the corresponding eigenvalues as algebraic numbers.
I've tried casting to QQbar and using eigenvectors, but this occasionally tries to divide by zero for no reason I can discern. Any ideas?