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

read `C5.txt`:

#1

RandomState := proc(n,K) local ra, k, V, i:
  ra := rand(-K..K):
  V := [seq(ra() + I * ra(), i=1..n)]:
  k := sqrt(add(abs(V[i])^2, i=1..n)):
  [seq(V[i] / k, i=1..n)]:
end:


#2

TestRHMEigen := proc() local i, L, ev, evc, v, V:
  for i from 1 to 3 do:
    L := RHM(3,10):
    ev, evc := op(EVC(L)):
    printf("\#i = %d, eigenvalues: [%Zf, %Zf, %Zf], inner products: [%Zf, %Zf, %Zf]\n", 
        i, evalf(ev[1]), evalf(ev[2]), evalf(ev[3]), 
        evalf(IP(evc[1], evc[2])), evalf(IP(evc[1], evc[3])), evalf(IP(evc[2], evc[3]))):
  od:
end:

#This outputted the following:

#i = 1, eigenvalues: [14.688338-0.000000I, -10.481598+0.000000I, 3.793261+0.000000I], inner products: [0.000000+0.000000I, 0.000000+0.000000I, 0.000000+0.000000I]
#i = 2, eigenvalues: [4.544407-0.000000I, -14.954649-0.000000I, -2.589758+0.000000I], inner products: [0.000000+0.000000I, 0.000000+0.000000I, 0.000000+0.000000I]
#i = 3, eigenvalues: [19.743843-0.000000I, -11.265375-0.000000I, 0.521532+0.000000I], inner products: [0.000000+0.000000I, 0.000000+0.000000I, 0.000000+0.000000I]


#3

TestRHMProbDist := proc() local i, L, A, pr:
  for i from 1 to 5 do:
    L := RHM(3,10):
    A := RandomState(3, 10): 
    pr := evalf(ProbDist(L, A)):
    printf("\#i = %d, ProbDist outputted: [%Zf, %Zf, %Zf], summing to: %Zf\n", i, pr[1], pr[2], pr[3], add(pr)):
  od:
end:

#This outputted the following:

#i = 1, ProbDist outputted: [0.000396+0.000000I, 0.949161+0.000000I, 0.050442+0.000000I], summing to: 1.000000+0.000000I
#i = 2, ProbDist outputted: [0.372660+0.000000I, 0.119462+0.000000I, 0.507878+0.000000I], summing to: 1.000000+0.000000I
#i = 3, ProbDist outputted: [0.497107+0.000000I, 0.310073+0.000000I, 0.192819+0.000000I], summing to: 1.000000+0.000000I
#i = 4, ProbDist outputted: [0.365543+0.000000I, 0.006527+0.000000I, 0.627930+0.000000I], summing to: 1.000000+0.000000I
#i = 5, ProbDist outputted: [0.814111+0.000000I, 0.013435+0.000000I, 0.172453+0.000000I], summing to: 1.000000+0.000000I


#5

GeneralPauli := proc(n) local nx, ny, nz:
  nx, ny, nz := op(n):
  [[nz, (nx - I * ny)],
   [(nx + I*ny), -nz]]:
end:


#simplify(subs(nx^2+ny^2+nz^2 = 1, EVC(GeneralPauli([nx, ny, nz])) )) outputs:

#[[1, -1], [[(1 + nz)*abs(nx + ny*I)/(sqrt(abs(1 + nz)^2 + abs(nx + ny*I)^2)*(nx + ny*I)), abs(nx + ny*I)/sqrt(abs(1 + nz)^2 + abs(nx + ny*I)^2)], 
#           [(-1 + nz)*abs(nx + ny*I)/(sqrt(abs(nx + ny*I)^2 + abs(-1 + nz)^2)*(nx + ny*I)), abs(nx + ny*I)/sqrt(abs(nx + ny*I)^2 + abs(-1 + nz)^2)]]]