{{{id=1|
# given two matrices, checks for conjugacy
def are_conjugate(A,B):
    if A.charpoly() != B.charpoly():
        return False
    else:
        d = -( A.charpoly()[0] )
        K.<u> = QuadraticField(d)
        IA = K.ideal( A[0][1] , u - A[0][0] )
        IB = K.ideal( B[0][1] , u - B[0][0] )
        if (IA/IB).is_principal():
            return True
        else: return False
///
}}}

{{{id=2|
# given two conjugate matrices, runs through U in GL_2(CC) satisfying A*U=U*B, until U is in GL_2(ZZ)
def matrix_U(A,B):
    if are_conjugate(A,B):
        u00 = -500
        while u00 < 1000:
            def u10plus(n):
                return (-1*A[0][0]*n + A[1][1]*n - sqrt( (n*(A[1][1]-A[0][0]))^2 + 4*A[0][1]*(-B[1][0]+A[1][0]*n^2) ) )/(2*A[0][1])
            def u10minus(n):
                return (-1*A[0][0]*n + A[1][1]*n - sqrt( (n*(A[1][1]-A[0][0]))^2 + 4*A[0][1]*(B[1][0]+A[1][0]*n^2) ))/(2*A[0][1])
            if u10plus(u00) in ZZ:
                u10 = u10plus(u00)
                u01 = (A[0][0]*u00 - B[0][0]*u00 + A[0][1]*u10)/B[1][0]
                u11 = (A[1][0]*u00 + A[1][1]*u10 - B[0][0]*u10)/B[1][0]
                if u01 in ZZ and u11 in ZZ:
                    return matrix(ZZ,2,[u00,u01,u10,u11])
            if u10minus(u00) in ZZ:
                u10 = u10minus(u00)
                u01 = (A[0][0]*u00 - B[0][0]*u00 + A[0][1]*u10)/B[1][0]
                u11 = (A[1][0]*u00 + A[1][1]*u10 - B[0][0]*u10)/B[1][0]
                if u01 in ZZ and u11 in ZZ:
                    return matrix(ZZ,2,[u00,u01,u10,u11])
            u00 += 1
    else:
        print "Not conjugate"
///
}}}

{{{id=3|
# generates a "random" element of SL_n(ZZ)
def rand_U():
    a = ZZ.random_element(-100,100)
    while a==0:
        a = ZZ.random_element(-100,100)
    b = ZZ.random_element(-100,100)
    while gcd(a,b)!=1 or b==0:
        b = ZZ.random_element(-100,100)
    d = inverse_mod(a,b)
    c = (a*d-1)/b
    return matrix(ZZ,2,[a,b,c,d])
///
}}}

{{{id=4|
# generates a "random" matrix with characteristic polynomial x^2 - w
def rand_matrix(w):
    a = ZZ.random_element(-100,100)
    bc = w - a^2
    if bc == 0:
        b = ZZ.random_element(-100,100)
        c = 0
    elif bc == 1:
        b = 1
        c = 1
    elif bc == -1:
        b = -1
        c = 1
    else:
        fac = factor(bc)
        exps = [ ZZ.random_element(0,x[1]+1) for x in fac ]
        b = 1
        for i in range(len(fac)):
            b *= (fac[i][0])**exps[i]
        c = bc/b
    return matrix(ZZ,2,[a,b,c,-a])
///
}}}

{{{id=5|
# test: conjugate a matrix A by a "randomly" chosen U to get B, and evaluate are_conjugate(A,B)
A = matrix(ZZ,2,[0,-19,1,0])
U = rand_U()
B = U*A*(U.inverse())
print B
print ""
are_conjugate(A,B)
///
[ 17329 -29383]
[ 10220 -17329]

True
}}}

{{{id=6|
# example: check whether two matrices with characteristic polynomial x^2 + 13 are conjugate over ZZ
A = rand_matrix(-13)
B = rand_matrix(-13)
print A
print ""
print B
print ""
are_conjugate(A,B)
///
[  35 1238]
[  -1  -35]

[ 12 157]
[ -1 -12]

True
}}}

{{{id=7|
# if possible, calculate U in GL_2(ZZ) such that A = U*B*(U^(-1))
if are_conjugate(A,B):
    U = matrix_U(A,B)
    print U
    print ""
    A == U*B*(U.inverse())
///
[-1 23]
[ 0 -1]

True
}}}