# hw5 Pablo Blanco
#
read(`C5.txt`): # to be accompanied with Dr. Z's Experimental Math, Spring 2025, C5 code
Randomize():

# outputs a random vector of length n, with complex entries in -K to K
# after being normalized
RandomState:= proc(n,K) local ra,v,i,m:
 ra:=rand(-K..K):
 v:=[seq(ra()+I*ra(),i=1..n)]:
 m:=add( abs(v[i])^2,i=1..n):
 v:=[seq(v[i]/(sqrt(m)),i=1..n)]:
 v:
end:

### Experiments:
### RHM(n,K) with n=3 and K=10
## Test Code (executed 3 times):
# A:=RHM(3,10);
# {seq(Im(EV(A)[i]) ,i=1..3)};  # check eigenvalues are real.
#
## Trial results: {0} in all cases. No imaginary components!
#
# evc := [seq( EVC1(A,simplify(EV(A)[i])),i=1..3)]; # all eigenvectors
            # {IP(evc[1],evc[2]),IP(evc[2],evc[3]), IP(evc[1],evc[3]) };
            ## the code above actually took too long to run. Instead, I ran the code below to approximate the inner product.
# {IP(evalf(evc[1]),evalf(evc[2])),IP(evalf(evc[2]),evalf(evc[3])), IP(evalf(evc[1]),evalf(evc[3])) };
## all inner products were approximately 0.
 

### Test Code (executed 5 times): 
# M:=RHM(3,10);
# A:=RandomState(3,10);
# P:=ProbDist(M,A):
# evalf(P[1]);       ## all components were positive
# add(evalf(P[1]));  ## sum was approximately 1
# 

# assume n is a normalized vector of length 3.
GeneralPauli:= proc(n) local nx,ny,nz:
 if nops(n)<>3 then: # decided not to check that it's normalized because we may have to use approximately normalized vectors
  RETURN(FAIL):
 fi:
 nx:=n[1]:
 ny:=n[2]:
 nz:=n[3]:
 [[nz,nx-I*ny] ,[nx+I*ny,-nz]]:
end:

### Tests (the following indeed return the Pauli matrices):
# matrix(GeneralPauli([1,0,0]));
# matrix(GeneralPauli([0,0,1]));
# matrix(GeneralPauli([0,1,0]));