# OK to post homework
# Aurora Hiveley, 2/11/25, Assignment 6

Help := proc(): print(`AllHM2(K), RealSV2(K), UncertaintyContest(K)`): end:

with(combinat):
with(linalg):

# AllHM2(K): inputs a positive integer K and outputs ALL 2 by 2 matrices of the form [[a,b+I*c],[b-I*c,d]]
# where a,b,c,d, are integers such that 1 ≤ a,b,c,d ≤ K (there are K^4 of them)

AllHM2 := proc(K) local a,b,c,d,M:
M := {}:

for a from 1 to K do 
    for b from 1 to K do 
        for c from 1 to K do 
            for d from 1 to K do 
                M := M union {[[a,b+I*c],[b-I*c,d]]}:
            od:
        od:
    od:
od:

M:
end:



# RealSV2(K): inputs a positive integer K and outputs all normalizations of the vectors (i.e. make them unit vectors)
# [a,b] where a and b are integers 1 ≤ a,b ≤ K (there are K2 of them

RealSV2 := proc(K) local a,b,k,U:
U := {}:

for a from 1 to K do 
    for b from 1 to K do 
        k := sqrt(a^2 + b^2):       # normalization factor
        U := U union {[a/k,b/k]}:
    od:
od:

U:
end:



# UncertaintyContest(K): inputs a positive integer K and outpts the triplet [A,B,V] (whose entries are in 1..K, V normalized)
# that make CheckUP(A,B,V) as small as possible and as big as possible (if A and B commute, do not count them).

UncertaintyContest := proc(K) local P,U,p,u,m,M,l:
P := choose(AllHM2(K),2):
U := RealSV2(K):

for p in P do 
    if Com(p[1],p[2])=[[0,0],[0,0]] then 
        P := P minus {p}:           # get rid of pairs which commute 
    fi:
od:

m := [CheckUP(P[1][1], P[1][2],U[1]),[P[1][1],P[1][2],U[1]]]:    # initialize max and min trackers
M := [CheckUP(P[1][1], P[1][2],U[1]),[P[1][1],P[1][2],U[1]]]:
for p in P do
    for u in U do 
        l := CheckUP(p[1],p[2],u):
        
        if l < m[1] then 
            m := [l, [p[1],p[2],u]]:
        fi:

        if l > M[1] then 
            M := [l, [p[1],p[2],u]]:
        fi:
        
    od:
od:

[m[2],M[2]]:
end:


# What are UncertaintyContest(2)?, UncertaintyContest(3)?
# UncertaintyContest(2) outputs:
# smallest: [[[1, 2 + I], [2 - I, 2]], [[2, 2 + I], [2 - I, 1]], [sqrt(2)/2, sqrt(2)/2]]
# biggest: [[[1, 1 + I], [1 - I, 1]], [[1, 1 + 2*I], [1 - 2*I, 2]], [(2*sqrt(5))/5, sqrt(5)/5]]
# UncertaintyContest(3) outputs: 
# smallest: [[[1, 2 + I], [2 - I, 3]], [[3, 3 + I], [3 - I, 1]], [sqrt(2)/2, sqrt(2)/2]]
# biggest: [[[1, 2 + 3*I], [2 - 3*I, 3]], [[2, 3 + 2*I], [3 - 2*I, 2]], [(3*sqrt(13))/13, (2*sqrt(13))/13]]
















### copied from C6.txt 
# C6.txt, 2/10/25 


#Feb. 10, 2025 C6.txt
Help6:=proc(): print(`Mul(A,B), Com(A,B), RSV(n,K), Del(L,A), CheckUP(A,B,V) `): end:

#Mul(A,B): the product of matrix A and B (assuming that it exists)
Mul:=proc(A,B) local i,j,k,n:
n:=nops(A):
[seq([seq( add(A[i][k]*B[k][j]   , k=1..n) , j=1..n)],i=1..n)]:
end:

#Com(A,B: The commutator of A and B
Com:=proc(A,B):    Mul(A,B)-Mul(B,A):end:

#RSV(n,K):
RSV:=proc(n,K) local ra,i,v,a:
ra:=rand(-K..K):
v:=[seq(ra()+I*ra(),i=1..n)]:
a:=sqrt(add(v[i]*conjugate(v[i]),i=1..n)):
[seq(v[i]/a,i=1..n)]:
end:

#Del(L,A): the "avergae error", i.e. standard deviation pf the prob. distribution of
#measuring the observable L when the state is A
Del:=proc(L,A) local n,pd,ev,i,av:
n:=nops(L):
pd:=ProbDist(L,A):
ev:=pd[1]:
av:=pd[3]:
pd:=pd[2]:
sqrt(add(ev[i]^2*pd[i],i=1..n)-av^2):
end:

#CheckUP(A,B,V): the ratio of Del(A,V)*Del(B,V)/Expectation(Com(A,B),V) (it should be >=1/2)
CheckUP:=proc(A,B,V) local C:
C:=Com(A,B):

if abs(ProbDist(C,V)[3])=0 then
 RETURN(FAIL):
fi:

evalf(Del(A,V)*Del(B,V)/abs(ProbDist(C,V)[3])):
end:

#old stuff
#C5.txt, Feb. 6, 2025 Experimental Mathematics (Rutgers, Dr. Z.)
Help5:=proc(): print(` PM(), RHM(n,K) , EV(L), EVC1(L,g), EVC(L), IP(A,B), ProbDist(L,A) `): end:

with(linalg):

#PM(): the three famous Pauli matrices sigma_x,sigma_y, sigma_z
PM:=proc():
[
[[0,1],[1,0]],
[[0,-I],[I,0]],
[[1,0],[0,-1]]
]:
end:

#RHM(n,K): a random n by n Hermitian matrix with entriries from -K to K
RHM:=proc(n,K) local ra,i,j,L:
ra:=rand(-K..K):
for i from 1 to n do
 for j from 1 to i-1 do
  L[i,j]:=ra()+I*ra():
  L[j,i]:=conjugate(L[i,j]):
 od:
L[i,i]:=ra():
od:
[seq([seq(L[i,j],j=1..n)],i=1..n)]:
end:

#EV(L): inputs a Hermitian matrix (say n by n) and outputs the LIST of eigenvalues
EV:=proc(L) local  x,L1,i,n:
n:=nops(L):
L1:=[seq(L[i]-[0$(i-1),x,0$(n-i)],i=1..n)]:
[solve(det(L1),x)]:
end:

#EVC1(L,g): inputs a matrix L and an eigenvalue g finds the NORMALIZED eigenvector
EVC1:=proc(L,g) local v,i,j,var,eq,V,n,L1,V1,va,k:
n:=nops(L):
L1:=[seq(L[i]-[0$(i-1),g,0$(n-i)],i=1..n)]:
if simplify(det(L1))<>0 then
  RETURN(FAIL):
fi:

V:=[seq(v[i],i=1..n)]:
var:={seq(v[i],i=1..n)}:
eq:={seq(add(L1[i,j]*V[j],j=1..n),i=1..n)}:
var:=solve(eq,var):
#V1 the set of unknowns that is arbitrary
V1:={}:
for va in var do
  if op(1,va)=op(2,va) then
  V1:=V1 union {op(1,va)}:
  fi:
od:
if nops(V1)<>1 then
  RETURN(FAIL):
fi:
V:=subs(var,V):
V:=subs(V1[1]=1,V):
k:=sqrt(add(abs(V[i])^2,i=1..n)):
[seq(V[i]/k,i=1..n)]:
end:

#EVC(L): inputs a Hermitian matrix and outputs the pair ListOfEigenvalues, ListOfEigenvectors
EVC:=proc(L) local n,ev,evc,i:
n:=nops(L):
ev:=EV(L):
evc:=[seq(EVC1(L,ev[i]),i=1..n)]:
[ev,evc]:
end:

#IP(A,B): the inner product of state-vectors A and B
IP:=proc(A,B) local n,i:
n:=nops(A):
if nops(B)<>n then
  RETURN(FAIL):
fi:
simplify(add(A[i]*conjugate(B[i]),i=1..n)):
end:

#ProbDist(L,A): Given an obervalble represented by the Hermitian matrix L and a state vector A
#outputs first all the eigenvectors (pure states) and the prob. of winding in the respective pure states
#followed by the expected value of A (i.e. if you make the observation L and the state is A the average in the long run
ProbDist:=proc(L,A) local evc,i,A1,k,PD,n:
n:=nops(L):
k:=sqrt(add(abs(A[i])^2,i=1..n)):
A1:=[seq(A[i]/k,i=1..n)]:
evc:=EVC(L):
PD:=[seq(abs(IP(A1,evc[2][i]))^2,i=1..n)]:
[evc[1],PD,add(PD[i]*evc[1][i],i=1..n)]:
end: