Computing conjugacy between two elements of $\mathrm{SL}_2(\mathbb{Z})$ The conjugacy classes of $\mathrm{SL}_2(\mathbb{Z})$ are well characterized (see, e.g., this question). Assuming two matrices $A, B \in \mathrm{SL}_2(\mathbb{Z})$ are conjugate, is there a way to compute their conjugacy, i.e. construct a matrix $C \in \mathrm{SL}_2(\mathbb{Z})$ such that $A C = C B$?
 A: If you have access to a soft as Maple , Mathematica or Sage, you can proceed as follows. Let $A,B\in \mathrm{SL}_2$ be conjugate.
Step 1. Using a Grobner basis library, you solve the equation in $X=\begin{pmatrix}a&b\\c&d\end{pmatrix}$ $XB-AX=0,\det(X)=1$ ($5$ equations in the $4$ unknowns $(a_{i,j}))$.
Since the dimension of the algebraic set of the solutions is $1$, you group the unknowns into 2 blocks, for example $(a,b)$ and $(c,d)$. In this way , you obtain $3$ equations, the first $2$ of degree $1$ and the last one, denoted by (*), is a quadratic function of $(c,d)$ ($uc^2+vd^2+wcd+t=0$). The calculation must be instantaneous with $70$ digits; otherwise change the couples of unknowns.
Step 2. You solve (*); you can do that on line here
https://www.alpertron.com.ar/QUAD.HTM
This type of equations admits an infinity of solutions except for example if we consider the equation $x^2+y^2+1=0$ which, I think so, does not have many solutions.
We obtain one or several primitive solutions and also one or several recursive solutions starting with primitive solutions. With $70$ digits, the time of calculation is $20"$.
Step 3.  The problem is that the first $2$ equations have the form: $pa+qb=r(c,d),p_1a+q_1b=r_1(c,d)$; then the solution obtained in the last equation must kill the denominator $pq_1-p_1q$.
Fortunately, in all the tests I carried out, a primitive solution was suitable. I don't know if this is general or at least generic; however, if not, then we have to use the recurrence formulas.
That follows is an example:



