# 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 #: 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: