# Ok to post homework
# Lucy Martinez, 02-27-2025, Assignment 10

with(combinat):

######################
#Problem 1: Write a progmram SimulateNAND(n,k, K)
# that for K times does the following
# picks a random P=RSLP(n,k)
# computes nops(EvalSP(P,n)) (a number from 0 to 2n)
# adds all of these numbers up
# divides by K

# Whad did you get for:
# evalf(SimulateNAND(10,5,1000)) = 696.3840000
# evalf(SimulateNAND(10,10,1000))= 661.8560000
# evalf(SimulateNAND(10,20,1000))= 644.7960000
# evalf(SimulateNAND(10,50,1000))= 632.6170000

SimulateNAND:=proc(n,k, K) local co,i,P,len:
co:=0:
for i from 1 to K do
 P:=RSLP(n,k):
 len:=nops(EvalSP(P,n)):
 co:=co+len:
od:
co/K:
end:






#Problem 2: Adapt RSLP(n,k) to write a procedure
# Monotone(n,k,p) where n and k are positive integers and
# p is a rational number between 0 and 1, that has no NOT,
# i.e. only using AND and OR, and a typical member is
# [i,j,AND] or [i,j,OR] (now i <= j, i=j is not allowed)
# and for each non-input gate, the probability of AND is p and
# the probability of OR is 1-p)

# then write SimulateMonotone(n,k, p, K)

# What did you get for:
# evalf(SimulateMonotone(10,5, 1/2,1000))=   710.5920000
# evalf(SimulateMonotone(10,10, 1/2, 1000))= 680.1040000
# evalf(SimulateMonotone(10,20, 1/2, 1000))= 648.7980000
# evalf(SimulateMonotone(10,50, 1/2, 1000))= 640.2320000
# evalf(SimulateMonotone(10,5, 3/4,1000))=   715.0400000
# evalf(SimulateMonotone(10,10, 3/4, 1000))= 675.3840000
# evalf(SimulateMonotone(10,20, 3/4, 1000))= 647.1540000
# evalf(SimulateMonotone(10,50, 3/4, 1000))= 629.9420000
# evalf(SimulateMonotone(10,5, 1,1000))=     714.3680000
# evalf(SimulateMonotone(10,10, 1, 1000))=   675.7200000  
# evalf(SimulateMonotone(10,20, 1, 1000))=   650.8400000
# evalf(SimulateMonotone(10,50,1, 1000))=    636.5570000

# Hint: In order to simulate a "loaded die" that
# returns true with probabality p (where p is a rational number between 0 and 1)
# you can use
# LD:=proc(p):evalb(rand(1..denom(p))()<=numer(p)):end:


Monotone:=proc(n,k,p) local P,i,i1,j1:
P:=[seq(i,i=1..n)]:
for i from n+1 to k+n do
 i1:=rand(1..nops(P))():
 j1:=rand(1..nops(P))():
 if i1<>j1 then
  if LD(p) then
   P:=[op(P),[i1,j1,AND]]:
  else
   P:=[op(P),[i1,j1,OR]]:
  fi:
 fi:
od:
P:
end:


SimulateMonotone:=proc(n,k, p, K) local co,i,P,len:
co:=0:
for i from 1 to K do
 P:=Monotone(n,k,p):
 len:=nops(EvalSP(P,n)):
 co:=co+len:
od:
co/K:
end:


LD:=proc(p):
evalb(rand(1..denom(p))()<=numer(p)):
end:

####################################previous classes
#C10.txt, Feb. 24, 2025
Help10:=proc():
print(`NAND(x,y),NOT(x),AND(x,y),OR(x,y),RSLP(n,k),TF(n)`):
print(`EvalSP1(P,IN),EvalSP(P,n)`):
end:


#NAND(x,y): not (x and y)= not x or not y 
NAND:=proc(x,y):
not (x and y):
end:

#NOT(x): uses NAND to get the boolean operation not
NOT:=proc(x):
#not (x and x)= not x or not x=not x
NAND(x,x):
end:

#AND(x,y): uses NAND to get the boolean operation x AND y
AND:=proc(x,y):
#NAND(x,y)= not x or not x
#NAND(NAND(x,y),NAND(x,y))=not (not x or not y) or not (not x or not y)
#  = (x and y) or (x and y)= x and y
NAND(NAND(x,y),NAND(x,y)):
end:

#x OR y= NOT( Not x AND Not y)
#OR(x,y): uses NAND to get the boolean operation x OR y
OR:=proc(x,y):
#NAND(x,x)=not x or not x=not x
#NAND(y,y)=not y or not y=not y
#NAND(NAND(x,x),NAND(y,y))=not(not x) or not(not y)=x or y
NAND(NAND(x,x),NAND(y,y)):
end:

#Def: Straight-line NAND program has n input gates labeled [1,...,n],
# and k (say) non-input gates of the form [i,j]
# where 1<=i<=j<=k-1 are previous gates
# [1,2,3,[1,1],[3,4],[1,3],[4,4]]
# To compute such a program you evaluate each gate in turn, and the last
# gate is the value of the Boolean function


#RSLP(n,k): A random straight-line progran with NAND only
RSLP:=proc(n,k) local P,i,i1,j1:
P:=[seq(i,i=1..n)]:

# this should output k+n-(n+1)+1=k gates
for i from n+1 to k+n do
 i1:=rand(1..nops(P))():
 j1:=rand(1..nops(P))():
 P:=[op(P),[i1,j1]]:
od:
P:
end:


#TF(n): the set of all true-false vectors of length n
TF:=proc(n) local V,v:
option remember:
if n=0 then
 return({[]}):
fi:

V:=TF(n-1):
[seq([op(v),false],v in V),seq([op(v),true],v in V)]:
end:


#EvalSP1(P,IN): inputs a straight-line NAND program and an input 
# true-false vector IN, outputs the value of the last gate
EvalSP1:=proc(P,IN) local P1,n,i,i1,j1:
P1:=IN:
n:=nops(IN):
for i from n+1 to nops(P) do
 i1:=P[i][1]:
 j1:=P[i][2]:
 P1:=[op(P1),NAND(P1[i1],P1[j1])]:
od:
P1[-1]:
end:


#EvalSP(P,n): Inputs a NAND straight-line program (SLP) and outputs
# the set of vectors of length n that yields TRUE
EvalSP:=proc(P,n) local V,v,T:
V:=TF(n):
T:={}:
for v in V do
 if EvalSP1(P,v) then
  T:=T union {v}:
 fi:
od:
T:
end:


###################################from previous classes

#C8.txt, Feb. 17, 2025
Help8:=proc():
print(`ES(),TP(A,B),ExpMV(M,V)`):
end:


#ES(): the four fully entangled states in the Alice-Bob combined system
# based on pg. 166 of Susskind-Friedman's book
# [singlet, T1,T2,T3]
# [uu, ud, du, dd]
ES:=proc():
[ 
 [0, 1/sqrt(2), -1/sqrt(2), 0],
 [0, 1/sqrt(2), 1/sqrt(2), 0],
 [1/sqrt(2),0,0,1/sqrt(2)],
 [1/sqrt(2),0,0,-1/sqrt(2)]
]:
end:


#TP(A,B): the tensor product of matrix A and matrix Bob
# (using our data structure of lists-of-lists)
# Let A be an a1 by a2 matrix and
# Let B be a b1 by b2 matrix
# the resulting matrix (after applying tensor prod) is dim a1*b1 by a2*b2
TP:=proc(A,B) local a1,a2,b1,b2,i,j,i1,j1,AB1,AB:
a1:=nops(A):
a2:=nops(A[1]):
b1:=nops(B):
b2:=nops(B[1]):

AB:=[]:
#the rows of TP(A,B) are called [i,i1] where 1<=i<=a1 and 1<=i1<=b1
#the cols of TP(A,B) are called [j,j1] where 1<=j<=a2 and 1<=j1<=b2
for i from 1 to a1 do
 for i1 from 1 to b1 do
  #AB1 is the [i,i1]-row of TP(A,B)
  AB1:=[]:
  for j from 1 to a2 do
   for j1 from 1 to b2 do
    AB1:=[op(AB1),A[i][j]*B[i1][j1]]:
   od:
  od:
  AB:=[op(AB),AB1]:
 od:
od:
AB:
end:


#Added after class:
#ExpMV(M,V): inputs  a Hermitian matrix M and a state vector V finds the expected value of the observable M in state V
#<conjuate(V),M,V>: 
ExpMV:=proc(M,V) local i,n,j:
add(add(conjugate(V[i])*M[i][j]*V[j],j=1..nops(V)),i=1..nops(V)):
end:

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

#Mul(A,B): the product of matrix A and B (assumming 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): A (normalized) random state vector with n complex components
# whose real and imaginary parts are from -K to K
RSV:=proc(n,K) local ra,v,i,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:


#SOME NOTES:
# mu:=Sum(a*Pr(X=a),a in A):
# av=average
# sigma^2=Sum((a-av)^2*Pr(X=a),a in A) = Sum(a^2*Pr(), a in A)
# The standard deviation of a random variable is the square root of its variance

#NOTE: Com(PM()[1],PM()[2])=2*I*PM()[3]

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


#Here, we will change notation by making 
#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:


#C5.txt, Feb. 6, 2025
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]],
[[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 eigenvalues, 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 outputs the pair ListOfEigenvalues and ListOfEigenvectors
evc:=EVC(L):
PD:=[seq(abs(IP(A1,evc[2][i]))^2,i=1..n)]:
#evc[1] is the list of eigenvalues
[evc[1],PD,add(PD[i]*evc[1][i],i=1..n)]:
end: