#1: RandomState(n,K)

#that inputs a pos. integer n and a pos. integer K and outputs a #random state vector with n complex components whose real and #imaginary parts are from -K to K, then normalize it.

RandomState:=proc(n,K):
vector:=[]:
for i from 1 to n do
 vector:=[op(vector), ra(-K..K)+I*ra(-K..K)]:
od:
vector:
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:


Experiment with three different random random Hermitian matrices using RHM(n,K) with n=3 and K=10 and check that the eigenvalues are real and that the eigenvectors (pure states) are all orhogonal to each other.

# I wrote a proc below to do this checking

with(LinearAlgebra) 

OrthogCheck:=(n,K) local matrices, i, M, m, j, k, dotproduct:

matrices:=[]:
for i from 1 to 3 do 
matrices:=[op(matrices), RHM(n,k)]:
i:=i+1:
od:

for M in matrices do
#vector of eigenvalues
V:=Eigenvectors(M)[1]:
#matrix of eigenvectors
E:=Eigenvectors(M)[2]:
for j from 1 to n do 
 if not Im(E[j])=0 then
	print("There is an eigenvalue with nonzero imaginary part"):
fi:
od:

for m from 1 to n do
 for p from m+1 to n do
  if not simplify(DotProduct(Column(V,j), Column(V,k)))=0 then 
   print("Eigenvectors are not orthogonal"):
  fi:
 od:
od:

end:



#__________________________________________________________



Experiment with five different random Hermitian matrices using RHM(n,K) with n=3 and K=10 and corresponding five random state vector and apply ProbDist(L,A) to them. For each case check that the proba. distribution consists of non-negative numbers that add up to 1.

#First we want to get 5 hermitian matrices and 5 random state vectors


matrices:=[]:
for i from 1 to 5 do 
M:=RHM(3, 10):
matrices:=[op(matrices), M]:
i:=i+1:
od:

vectors:=[]:
for i from 1 to 5 do
vector:=RandomState(3,10):
vectors:=[op(vectors), vector]:
i:=i+1:
od:

for L in matrices do
 for A in vectors do
#for n=3 the matrix is nxn so there are 3 eigenvectors and 3 elts in #probdistribution 
  SUM:=add(ProbDist(L,A)[1][1], ProbDist(L,A)[1][2], ProbDist(L,A)[1][3]):
  if SUM <> 1 then return FAIL:
  fi:
  if SUM = 1 then print("probabilities add to one"):
  fi:
  for i from 1 to 3 do 
   if ProbDist(L,A)[1][i] <0 return FAIL:
   fi:
   if ProbDist(L,A)[1][i] >= 0 print("probs all nonnegative") :
  od:
od:
od:


#__________________________________________________________


Write a procedure
GeneralPauli(n)

that inputs a NORMALIZED vector n=[nx,ny,nz]

that constructs the 2 by 2 matrix. at the bottm of p. 349 of the Susskind-Frieman Quantum Mechanics book book.

#general pauli matrices are of the form [nz, nx-iny][nx+iny, nz]
#couldn't find on given page but looked up online

Find its eigenvalues and eigenvectors.

proc(n):=GeneralPauli(n)
nx:=n[1]:
ny:=n[2]:
nz:=n[3]:

N:=Matrix([[nz, nx-iny],[nx+iny, nz]]):

E:=Eigenvectors(N)[1]:
V:=Eigenvectors(N)[2]:

end: