#OK to post homework
#Joseph Koutsoutis, 02-09-2025, Assignment 4

read `C4.txt`:

#3

UM := proc(alpha, beta, gamma, delta) local U, i, j, M1, M2, M3:
  M1 := [[exp(-I * beta / 2), 0], [0, exp(I * beta / 2)]]:
  M2 := [[cos(gamma / 2), -sin(gamma/2)], [sin(gamma/2), cos(gamma / 2)]]:
  M3 := [[exp(-I * delta / 2), 0], [0, exp(I * delta / 2)]]:
  U := multiply(M2, M3):
  U := multiply(M1, U):
  [seq( [seq(exp(I * alpha) * U[i,j], j=1..2)] , i=1..2)]:
end:

#TestUM(n, d): Takes as input positive integers n,d and outputs multiply(U, CT(U)) where 
#U is the output of UM(alpha, beta, gamma, delta) with 
#alpha, beta, gamma, delta being randomly generated rational numbers with 
#numerator between -n and n and 
#denominator between 1 and d.
TestUM := proc(n, d) local ra1, ra2, U:
  ra1 := rand(-n..n):
  ra2 := rand(1..d):
  U := UM(ra1() / ra2(), ra1() / ra2(), ra1() / ra2(), ra1() / ra2()):
  evalf(multiply(U, CT(U))):
end:


#I ran TestUM(n,d) a couple of times and it outputted the identity matrix as desired.
#More generally, we compute:
#   multiply(U, CT(U)) = exp(I * alpha)exp(-I * alpha)M1(beta)M2(gamma)M3(delta)CT(M3(delta))CT(M2(gamma))CT(M1(beta))
#                      = M1(beta)M2(gamma) [[exp(-I * delta / 2), 0], [0, exp(I * delta / 2)]][[exp(I * delta / 2), 0], [0, exp(-I * delta / 2)]] CT(M2(gamma))CT(M1(beta))
#                      = M1(beta)M2(gamma) [[exp(-I * delta / 2)exp(I * delta / 2), 0], 
#                                           [0, exp(I * delta / 2)exp(-I * delta / 2)]] CT(M2(gamma))CT(M1(beta))
#                      = M1(beta) [[cos(gamma / 2), -sin(gamma/2)], [sin(gamma/2), cos(gamma / 2)]][[cos(gamma / 2), sin(gamma/2)], [-sin(gamma/2), cos(gamma / 2)]] CT(M1(beta))
#                      = M1(beta) [[cos^2(gamma / 2) + sin^2(gamma/2), cos(gamma/2)sin(gamma/2) - sin(gamma/2)cos(gamma/2)], 
#                                  [sin(gamma/2)cos(gamma/2) - cos(gamma/2)sin(gamma/2), sin^2(gamma / 2) + cos^2(gamma / 2)]] CT(M1(beta))
#                      = [[exp(-I * beta / 2), 0], [0, exp(I * beta / 2)]][[exp(I * beta / 2), 0], [0, exp(-I * beta / 2)]]
#                      = [[exp(-I * beta / 2)exp(I * beta / 2), 0], 
#                         [0, exp(I * beta / 2)exp(-I * beta / 2)]]
#                      = [[1,0],
#                         [0,1]]