Similarity Of Matrices
I don't get very well the point of the requirement in the problem of a randomized algorithm. Or maybe I am misinterpreting what is the notion of problem size with respect to which the running time is required to be polynomially bounded. With respect to the number of entries of $A$ and $B$ there is the following deterministic algorithm.
The requirement of $C$ having rational entries, or more generally entries in the same field as the entries of $A$ and $B$ points to using the Frobenius/Rational normal form (FNF). The Jordan normal form is inconvenient because, in general, it requires working in a field extension in which the characteristic polynomials of $A$ and $B$ factor completely (split).
Just like in the case of Jordan normal form, there are invertible matrices $P_1,P_2$ such that $P_1AP_1^{-1}=A_f$ and $P_2BP_2^{-1}=B_f$ and $A_f, B_f$ are the Frobenius normal forms of $A,B$, respectively. The additional property that we mentioned before is that all these matrices have entries in the same field as the entries of $A$ and $B$.So, if your matrices have integer entries, then these matrices will have rational entries.
Now, the matrices $A$ and $B$ are similar if and only if $A_f=B_f$. When this happens we can compute the required $C$ from $P_1AP_1^{-1}=A_f=B_f=P_2BP_2^{-1}$, since we get $(P_2^{-1}P_1)A(P_2^{-1}P_1)^{-1}=B$. So, we can take $C=P_2^{-1}P_1$.
The FNF can be computed (without emphasis on efficiency) by considering the polynomial matrices $A-xI$ and $B-xI$ and reducing them to their Smith normal form. This consists in a combination of Gauss elimination (row/column reduction, if you call it that way) combined with Euclid algorithm to compute some polynomial GCDs that are required. The Smith normal form tells you the invariant factors, and then the Frobenius normal form is the block-diagonal matrix in which the blocks are the companion matrices of the invariant factors. The sequence of steps to reduce $A-xI$ and $B-xI$ to their Smith normal forms give you the $P_1$ and $P_2$ required above.
Justifying all these claims requires a long exposition of the linear algebra (and either explicitly or implicitly the structure theorem of finitely generated modules over a principal ideal domain). This can be found in many linear algebra books. One that I know has it is Dummit and Foote's Abstract Algebra.

A guess of how randomization could be used is by trying random vectors $v$ and finding the maximal Krylov subspaces that it generates. This is, the space generated by $v,Av,A^2v,A^3v,...$, which ought to be generated by finitely many of those vectors. Each one of those give you a cyclic invariant subspace of $A$, which tell you a block of the Frobenius normal form. Then pass to the quotient space by this subspace and repeat until you find all blocks of the FNF. Understanding the theory is required again to see why.