Do you mean that it may be limited by approximate accuracy, say, an exactly non-zero value been evaluated to zero internally, during the computation process?

I assume that JordanDecomposition in Mathematica works with arbitrary complex numbers, exact or approximate, not just rationals or algebraics:

JordanDecomposition[{{E, Pi}, {ArcTan[2], Sin[I]}}]

I've noticed the following technical explanation of the Jordan form based method implemented in sage math:

Warning

When the two matrices are similar, this routine may fail to find the similarity transformation. A technical explanation follows.

The similarity check is accomplished with rational form, which will be successful for any pair of matrices over the same field. However, the computation of rational form does not provide a transformation. So we instead compute Jordan form, which does provide a transformation. But Jordan form will require that the eigenvalues of the matrix can be represented within Sage, requiring the existence of the appropriate extension field. When this is not possible, a RuntimeError is raised, as demonstrated in an example below.

I wonder if the same limitation applies for the method discussed here.

Regards, HZ

The differences? Gianlucas method use 2 Jordandecomposition, my method translates your problem to a sort of Gauss-elimination, which after all could be described as similar. Pros and cons? I don't know, I just wanted to see if this method works.

You could run each method (my could be a bit shorter: don't calculate the determinats of mats) a 1000 (or so) times and look which is faster.

Hi Daniel Lichtblau,

Exactly, the problem's background is here:

Suppose I already have a unitary irreducible complex representation for a group. In this case, can I find the conjugate/similar transformation matrices connecting the two unitary irreducible complex representations, i.e., the one used by GAP and the one I already have, of the same group?

You may try with JordanDecomposition:

a = {{-I, 0, 0, 0}, {0, I, 0, 0}, {0, 0, 0, -I}, {0, 0, -I, 0}};
b = {{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, I}, {0, 0, I, 0}};
{a0, a1} = JordanDecomposition[a];
{b0, b1} = JordanDecomposition[b];
a1 == b1
b == b0 . b1 . Inverse[b0] == b0 . Inverse[a0] . a . a0 . Inverse[b0]

It seems that the is_similar function provided by sagemath is exactly for this purpose, as described here. But I'm not sure if there is an equivalent implementation in Wolfram.

Regards, HZ