# OK to post homework # Aurora Hiveley, 3/7/25, Assignment 13 Help := proc(): print(`quadruples(), quintuples(), divisor(C,A), ApplyT(T,X,n,u)`): end: with(combinat): with(linalg): ### HXH = Z # Find all triples 2 ≤ i,j,k ≤ 7 (there are 6^3=216 such triples) such that # Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k]) is equal to a member of C1QG(), say C1QG()[r]. # example: Mul(Mul(C1QG()[7],C1QG()[2]),C1QG()[7])=C1QG()[4] # Output the set of all such quadruples [i,j,k,r] (in particular [7,2,7,4] should be a member of that set). quadruples := proc() local C,S,Q,i,j,k,s,r,A: S := {seq(seq(seq([i,j,k],i=2..7),j=2..7),k=2..7)}; # set of all possible triples Q := {}: # initialize set of quadruples C := C1QG(): for s in S do A := Mul(Mul(C[s[1]],C[s[2]]),C[s[3]]): if A in C then member(A,C,'r'): Q := Q union {[op(s), r]}: fi: od: Q: end: # How big is that set? # Q has 54 quadruples: # {[2, 2, 2, 2], [2, 2, 3, 3], [2, 2, 4, 4], [2, 2, 5, 5], [2, 2, 6, 6], [2, 2, 7, 7], [2, 3, 3, 2], [2, 4, 4, 2], [2, 7, 4, 7], [2, 7, 7, 2], # [3, 2, 2, 3], [3, 3, 2, 2], [3, 3, 3, 3], [3, 3, 4, 4], [3, 3, 5, 5], [3, 3, 6, 6], [3, 3, 7, 7], [3, 4, 4, 3], [3, 5, 2, 5], [3, 7, 7, 3], # [4, 2, 2, 4], [4, 3, 3, 4], [4, 4, 2, 2], [4, 4, 3, 3], [4, 4, 4, 4], [4, 4, 5, 5], [4, 4, 6, 6], [4, 4, 7, 7], [4, 5, 4, 5], [4, 5, 5, 1], # [4, 6, 4, 6], [4, 7, 2, 7], [4, 7, 7, 4], [5, 2, 2, 5], [5, 3, 3, 5], [5, 4, 4, 5], [5, 4, 5, 1], [5, 5, 4, 1], [5, 7, 7, 5], [6, 2, 2, 6], # [6, 3, 3, 6], [6, 4, 4, 6], [6, 7, 7, 6], [7, 2, 2, 7], [7, 2, 7, 4], [7, 3, 3, 7], [7, 4, 4, 7], [7, 4, 7, 2], [7, 7, 2, 2], [7, 7, 3, 3], # [7, 7, 4, 4], [7, 7, 5, 5], [7, 7, 6, 6], [7, 7, 7, 7]} # adapted to find Mul(Mul(C1QG()[i],C1QG()[j]),C1QG()[k])=c*C1QG()[r] quintuples := proc() local C,S,Q,i,j,k,s,c,d,r,A: S := {seq(seq(seq([i,j,k],i=2..7),j=2..7),k=2..7)}; # set of all possible triples Q := {}: # initialize set of quintuples C := C1QG(): for s in S do A := Mul(Mul(C[s[1]],C[s[2]]),C[s[3]]): for r from 2 to 7 do c := C[r]: d := divisor(c,A): # helper function, returns common divisor if d<>FAIL then Q := Q union {[op(s), d,r]}: fi: od: od: Q: end: # there are 248 such quintuples # if A = mC for some constant m, returns m. else returns fail divisor := proc(C,A) local ind,div,i: # written specifically for two 2x2 matrices ind := [[1,1],[1,2],[2,1],[1,2]]: div := {}: for i in ind do if C[i[1]][i[2]] = 0 and A[i[1]][i[2]]<>0 then RETURN(FAIL): elif C[i[1]][i[2]]<>0 and A[i[1]][i[2]] = 0 then RETURN(FAIL): elif C[i[1]][i[2]] = 0 and A[i[1]][i[2]] = 0 then div := div: else div := div union {simplify(A[i[1]][i[2]] / C[i[1]][i[2]])}: fi: od: if nops(div)<>1 then RETURN(FAIL): fi: op(div): end: # ApplyT(T,X,n,u): has T=[T1,V], T1 is a table and V is nBasis(n,X), and u is an arbitrary linear combination of of members of nBasis(n,X). # uses the linearity of the table, and outputs the result of applying the linear transformation T (given as a table defined on the canonical # basis), as a linear combination of the basis elements. ApplyT := proc(T,X,n,u) local T1,V,v,i: T1 := T[1]: # transformation V := T[2]: # basis v := u: # initialize output vector pre-substitution # replace the basis elements in u with dummy variables to avoid double substitution for i from 1 to 2^n do v := subs(V[i]=a[i],v): od: # apply the transformation of each basis element for i from 1 to 2^n do v := subs(a[i] = T1[T[2][i]],v): od: simplify(v): end: # If T=UMtoT(TP(H,Z)), what is # ApplyT(T,X,2,X[[0,0]]+I*X[[0,1]]+2*I*X[[1,0]]+5*X[[1,1]]); ? # output is: # (((5 - I)*X[[1, 1]] - (5 + I)*X[[0, 1]] + (1 - 2*I)*X[[1, 0]] + (1 + 2*I)*X[[0, 0]])*sqrt(2))/2 ### copied from class #March 6, 2025 Computing Quantum gates Help13:=proc(): print(`C1QG(), C1QGnames(), BV(n), nBasis(n,X), UMtoT(M,n,X) , TtoUM(T,n) `): end: C1QGnames:=proc():[` Identity `, ` PauliX `, `PauliY `, `PauliZ `, ` PhaseS `, `TPi/8` , ` Hadamard `]: end: C1QG:=proc(): [ [[1,0],[0,1]], [[0,1],[1,0]], [[0,-I],[I,0]], [[1,0],[0,-1]], [[1,0],[0,I]], [[1,0],[0,sqrt(2)/2+sqrt(2)/2*I]], expand([[1,1],[1,-1]]/sqrt(2)) ]: end: #BV(n): inputs a non-neg. integer n and outputs the list of length 2^n of all 0-1 vectors in lex. order BV:=proc(n) local V,v,i: option remember: if n=0 then RETURN([[]]): fi: V:=BV(n-1): [ seq([0,op(V[i])],i=1..nops(V)),seq([1,op(V[i])],i=1..nops(V))]: end: #nBasis(n,X): the basis of n-quantum gates circuits X[e1,e2,.., en] corresponds to |e1>|e2>...|en> nBasis:=proc(n,X) local V,i: V:=BV(n): [seq(X[V[i]],i=1..nops(V))]: end: #UMtoT(M,n,X): inputs a unitary matrix M of size 2^n by 2^n (correpsonding to a quantum circuit with n #qbits outputs the corresponding LINEAR TRANSFORMATION as it acts on the canonical basis #nBasis(n,X) UMtoT:=proc(M,n,X) local T,V,i,j: if not IsUM(M) then print(M, `is not unitary `): RETURN(FAIL): fi: if nops(M)<>2^n then RETURN(FAIL): fi: V:=nBasis(n,X): for j from 1 to nops(V) do T[V[j]]:=add( V[i]*M[i][j],i=1..nops(M)): od: [op(T),V]: end: TtoUM:=proc(T,n) local T1,B,i,j: B:=T[2]: T1:=T[1]: if nops(B)<>2^n then RETURN(FAIL): fi: [seq ([ seq ( coeff(T1[B[j]],B[i],1),j=1..2^n)],i=1..2^n)]: end: ####from previous classes HelpOld:=proc(): print(`CT(A), IsUM(A), Mul(A,B), ES(), TP(A,B)`):end: #Mul(A,B): the product of matrix A and B (assuming that it exists) Mul:=proc(A,B) local i,j,k: [seq([seq(add(A[i][k]*B[k][j],k=1..nops(A[i])),j=1..nops(B[1]))],i=1..nops(A))]: end: #CT(A): the conjugate transpose of A CT:=proc(A): local i,j,n:n:=nops(A):[seq([seq(conjugate(A[j][i]),j=1..n)],i=1..n)]:end: Con:=proc(A): local i,j,n:n:=nops(A):[seq([seq(A[j][i],j=1..n)],i=1..n)]:end: #IsUM(A): Is the matrix A unitary IsUM:=proc(A) local n,i,j, P: n:=nops(A): P:=simplify(Mul(A,CT(A))): evalb(P=[seq([0$(i-1),1,0$(n-i)],i=1..n)]): end: #ES(): the four fullly entangled states in the Alice-Bob combined system #based on p. 166of 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 B (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 #TP(A,B): is a1*b1 by a2*b2 matrix TP:=proc(A,B) local a1,a2,b1,b2,AB,i,i1,j,j1,AB1: a1:=nops(A): a2:=nops(A[1]): b1:=nops(B): b2:=nops(B[1]): #The rows of TP(A,B) are called [i,i1] where 1<=i<=a1 and 1<=i1<=b1 #The columns of TP(A,B) are called [j,j1] where 1<=j<=a2 and 1<=j1<=b2 AB:=[]: 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: ##End Old stuff