# Please do not post
# Matthew Esaia, 2/9/2025, Homework 5

with(linalg):

Help:=proc(): print(`RandomState(n,K), GeneralPauli(n), TestPaul(), P2(n,K)`); end:

## PROBLEM 1
# RandomState(n,K), outputs a random state vector with n complex components from -K to K, normalized
RandomState:=proc(n,K) local ra,v,norm,vnorm,i:
	ra:=rand(-K..K):
	v:= [seq(ra() + I*ra(), 1..n)]:
	
	norm:=sqrt(add(abs(v[i])^2,i=1..n)):
	vnorm:=[seq(simplify(v[i]/norm),i=1..n)]:
end:

## PROBLEM 2

# RHM(n,K): a random n by n Hermitian matrix with entries 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:

# 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:
	add(A[i]*conjugate(B[i]),i=1..n):
end:

P2 := proc(n,K) local L,i,vals,vecs,num,vi,vj:
	for i from 1 to 3 do:
		L:=RHM(n,K):
		
		vals := eigenvalues(L):
		vecs := eigenvectors(L):
		
		print(vals);

		for vi in vecs do:
			for vj in vecs do:
				if vi <> vj then
					if IP(vi, vj) = 0 then
						print(`TRUE`);
					else
						print(`FALSE`):
					fi:
				fi:
			od:
		od:
	od:
end:


## PROBLEM 4
# GeneralPauli(n), inputs a normalized vector, constructs 2x2 matrix and finds eigenvalues + eigenvectors
GeneralPauli:=proc(n) local sigma, x, y, z, vals, vecs:
	x := n[1]:
	y := n[2]:
	z := n[3]:
	sigma := [[z, x - I * y], [x + I * y, -z]];
	
	vals := eigenvalues(sigma);
	vecs := eigenvectors(sigma);
	
	print(sigma);
	print(vals);
	print(vecs);
end:

TestPauli := proc() local ra, i, v, norm, vnorm:
	ra:=rand(-100..100):
	v := [seq(ra(), 1..3)]:
	
	norm:=sqrt(add(abs(v[i])^2,i=1..3)):
	vnorm:=[seq(simplify(v[i]/norm),i=1..3)]:
	
	GeneralPauli(vnorm):
end:

## The eigenvalues are always 1 and -1, with eigenvectors [1, 1] and [-1, 1]