######################################################################
## QC.txt Save this file as  QC.txt  to use it,                      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `QC.txt` <Enter>                             #
## Then follow the instructions given there                          #
##                                                                   #
## Written by the Experimental Mathematics class taught by Dr. Z.    #
# (Doron Zeilberger), Rutgers University , Spring 2025               #
######################################################################


print(`First Written: March-May tested for Maple 2024 `):
print(`UNDER CONSTRUCTION`):

print(`Version  : March 22, 2025  `):
print():
print(`Written by the Experimental Mathematics class taught by Dr. Z. (Doron Zeilberger), Rutgers University , Spring 2025`):

print(`This is  QC.txt, A Maple package to simulate,in a classic computer, quantum circuits and algorithms`):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/EM25/QC.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "Help();". For specific help type "Help(procedure_name);" `):
print(`For a list of the supporting functions type: Help1();`):
print():
 



Help1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:`):
print(` ApplyT, C1QGnames(), CT, Embed1, ES, MergeV, nBasis, TtoUM `):

else
Help(args):

fi:


end:

Help:=proc()
if args=NULL then

print(` : A Maple package for studying and simulating Quantum circuits and algorithms (under construction) `):

print(`The MAIN procedures are:`):
print(` BV, C1QG, Factorize, IsUM, ModExp, Mul, Ord,  QCtoM, RQC, TP, UM, UMtoT `):


elif nargs=1 and args[1]=ApplyT then
print(`ApplyT(T,u) : Inputs a pair T=[Table,V] where V is the domain of the Table Table, describing the action of the linear transformation,let's call it T1,  on the set of`):
print(`basis elements V.`):
print(`It also inputs u, a linear combination of the members of V, Outputs T1(u), a certain the member of the vector space generated by V`):
print(`[Adapted From Joseph Koutsoutis's hw13 (see https://sites.math.rutgers.edu/~zeilberg/EM25/hw13posted/hw13JosephKoutsoutis.txt)]. Try: `):
print(`T:=UMtoT(C1QG()[7],1,X);ApplyT(T,X[[0]]+4*X[[1]]);`):

elif nargs=1 and args[1]=BV then
print(`BV(n): inputs a non-neg. integer n and outputs the list of length 2^n of all 0-1 vectors in lex. order. Try:`):
print(`BV(3);`):

elif nargs=1 and args[1]=C1QG then
print(`C1QG(): the list of the seven most imporant one-qubit gates (unitary 2 by 2 matrices) in lists-of-lists format`):
print(`in that order: Identity , PauliX , PauliY , PauliZ ,  PhaseS , TPi/8 ,  Hadamard `):
print(`Try:`):
print(`C1QG();`):


elif nargs=1 and args[1]=C1QGnames then
print(`C1QnamesG(): the names the seven most imporant one-qubit gates (unitary 2 by 2 matrices) in lists-of-lists format`):
print(`in that order: Identity , PauliX , PauliY , PauliZ ,  PhaseS , TPi/8 ,  Hadamard `):
print(`Try:`):
print(`C1QGnames();`):

elif nargs=1 and args[1]=CT then
print(`CT(A): the conjugate transpose of A. Try:`):
print(`CT([[a11,a12],[a21,a22]]);`):

elif nargs=1 and args[1]=Embed1 then
print(`Embed1(n,S,M): inputs a pos. integer n and a subset S given as list of length k, and a 2^k by 2^k matrix M, acting on the qubits in S`):
print(`outputs the  correponding 2^n by 2^n matrix where the other qubits are the same. Try:`):
print(`Embed1(4,[1],[[a11,a12],[a21,a22]]);`):

elif nargs=1 and args[1]=ES then
print(`ES(): the four fullly entangled states in the Alice-Bob combined system`):
print(`based on p. 166of Susskind-Friedman's book`):
print(`[singlet, T1,T2,T3]`):
print(`[uu.ud,du,dd]. Try:`):
print(`ES();`):

elif nargs=1 and args[1]=Factorize then
print(`Factorize(n,K): inputs a pos. integer n and outputs the list of prime factors. K is a guessing parameter. Try:`):
print(`Factorize(1001,10);`):

elif nargs=1 and args[1]=IsUM then
Print(`IsUM(A): Is the matrix A (given as a list-of-lists) unitary. Try:`):
print(`IsUM(UM(Pi/3,Pi/4,Pi/5,Pi/3));`):


elif nargs=1 and args[1]=MergeV then
print(`MergeV(v1,v2,S1,S2): Given vector v1 supported on a list S1, and a vector v2 supported on a list S2, let n:=nops(v1)+nops(v2)  such that {op(S1) union op(S2)}={1,.,n}`):
print(`outputs the vector v such that v[S1[i]]=v1[i] and v[S2[i]]:=v2[i]. Try:`):
print(`MergeV([0,1],[1,0],[1,3],[2,4]);`):

elif nargs=1 and args[1]=ModExp then
print(`ModExp(x,r,n): x^r mod n the fast way where r is given in reverse binary [2^0,2^1,2^2,...]. Try:`):
print(`ModExp(2,[1,0,1,0,1],101);`):

elif nargs=1 and args[1]=Mul then
print(`Mul(A,B): the product of matrix A and B (assuming that it exists), using the list-of-lists format. Try:`):
print(`Mul([[a1,a2]],[[b1],[b2]]);`);
print(`Mul([[b1],[b2]],[[a1,a2]]);`);


elif nargs=1 and args[1]=nBasis then
print(`nBasis(n,X): the basis of n-quantum gates circuits X[e1,e2,.., en] corresponds to |e1>|e2>...|en>. Try:`):
print(`nBasis(3,X);`):

elif nargs=1 and args[1]=Ord then
print(`Ord(x,n): inputs a large pos. integer n and x<n m such that gcd(x,n) outputs (by brute force)`):
print(`the order of x mod n, i.e. the smallest pos. integer r such that x^r mod n=1. Try:`):
print(`Ord(2,101);`):

elif nargs=1 and args[1]=QCtoM then
print(`QCtoM(C): inputs  a Quantum circut C=[n,ListOfGates] where there are n qubits and ListOfGates is a list where each entry is a pair [S,Matrix] where S is a subset of`):
print(`{1,...n} indicating the active qubits and M is a 2^nops(S) by 2^nops(S) unitary matrix describing the gate's action. Try:`):
print(`C:=RQC(4,5,2); QCtoM(C);`):


elif nargs=1 and args[1]=RQC then
print(`RQC(n,r,k): generates a random quantum circuit with n qubits of length r, with gates involving at most k qubits at a time`):
print(`obtained by randomly tensor-producting C1QG()[i]. The outputs is a list of length r where S_i is the support of M_i`):
print(`[n,[seq([S_i,M_i],i=1..r)]]. Try:`):
print(`RQC(4,5,2);`):

elif nargs=1 and args[1]=TP then
print(`TP(A,B):the tensor product of matrix A and matrix B (using our data structure`):
print(`of lists-of-lists`):
print(`Let A be an a1 by a2 matrix and `):
print(`Let B be a b1 by b2 matrix`):
print(`TP(A,B): is a1*b1 by a2*b2 matrix. Try`):
print(`TP([[a11,a12],[a21,a22]],[[b11,b12],[b21,b22]]);`):

elif nargs=1 and args[1]=TtoUM then
print(`TtoUM(T,n): inputs a pair T=[Table,Basis], where Basis is a list of basis elements generating the underlying vector space, and Table is a table such that`):
print(`Table[Basis[i]] is the linear combination of the elements of Basis that describes the action of the operator on the basis element Basis[i]. `):
print(`Outputs the corresponding 2^n by 2^n matrix . Try:`):
print(`M:=C1QG()[6]; T:=UMtoT(M,1,X); evalb(TtoUM(T,1)=M);`):
print(`M:=TP(C1QG()[6],C1QG()[4]); T:=UMtoT(M,2,X); evalb(TtoUM(T,2)=M);`):

elif nargs=1 and args[1]=UM then
print(` UM(alpha,beta,gamma,delta):  implements box 1.1. on p. 20 of Nielsen-Chuang `):
print(` Adapted from Auorora Hively's hw4  (See https://sites.math.rutgers.edu/~zeilberg/EM25/hw4posted/hw4AuroraHiveley.txt )`):
print(`note: alpha, beta, delta, gamma are all real valued. Try:`):
print(` UM(a,b,c,d);`):
print(` UM(Pi/4,Pi/3,Pi/2,Pi/5);`):

elif nargs=1 and args[1]=UMtoT then
print(`UMtoT(M,n,X): inputs a unitary matrix M of size 2^n by 2^n (correpsonding to a quantum circuit with n`):
print(`qubits.Outputs the corresponding LINEAR TRANSFORMATION as it acts on the canonical basis, it also returns the list consisting of the canonical basis`):
print(`indexed by X, namely nBasis(n,X); Try:`):
print(`M:=C1QG()[7]; UMtoT(M,1,X);`):
print(`M:=TP(C1QG()[7],C1QG()[5]); UMtoT(M,2,X);`):

else
 print(`There is no such thing as`, args):

fi:


end:




with(combinat):

#MergeV(v1,v2,S1,S2): Given vector v1 supported on a list S1, and a vector v2 supported on a list S2, let n:=nops(v1)+nops(v2)  such that {op(S1) union op(S2)}={1,.,n}
#outputs the vector v such that v[S1[i]]=v1[i] and v[S2[i]]:=v2[i]. Try:
#MergeV([0,1],[1,0],[1,3],[2,4]);
MergeV:=proc(v1,v2,S1,S2) local n,v,i:
n:=nops(v1)+nops(v2):
if {op(S1),op(S2)}<>{seq(i,i=1..n)} then
 RETURN(FAIL):
fi:

if {op(S1)} intersect {op(S2)}<>{} then
 RETURN(FAIL):
fi:

for i from 1 to nops(S1) do
v[S1[i]]:=v1[i]:
od:

for i from 1 to nops(S2) do
v[S2[i]]:=v2[i]:
od:

[seq(v[i],i=1..n)]:
end:

#Embed1(n,S,M): inputs a pos. integer n and an ordered subset S given as list of length k, and a 2^k by 2^k matrix M, acting on the qubits in S
#outputs the  correponding 2^n by 2^n matrix where the other qubits are the same. Try:
#Embed1(4,[1],[[a11,a12],[a21,a22]]);
Embed1:=proc(n,S,M) local k,T,Bk,Bk1, Bn, T1, i, i1, j, j1,S1:
k:=nops(S):

S1:=sort(convert({seq(i,i=1..n)} minus {op(S)},list)):
if not (type(M,list) and nops(M)=2^k and nops(M[1])=2^k) then
 RETURN(FAIL):
fi:

Bk:=BV(k):

Bk1:=BV(n-k):

for i from 1 to nops(Bk) do
for j from 1 to nops(Bk) do
 T1[Bk[i],Bk[j]]:=M[i][j]:
od:
od:

Bn:=BV(n):

for i from 1 to nops(Bn) do
for j from 1 to nops(Bn) do
 T[Bn[i],Bn[j]]:=0:
od:
od:

for i from 1 to nops(Bk1) do
for i1 from 1 to nops(Bk) do
for j1 from 1 to nops(Bk) do
T[MergeV(Bk[i1],Bk1[i],S,S1),MergeV(Bk[j1],Bk1[i],S,S1)]:=T1[Bk[i1],Bk[j1]]:
od:
od:
od:

[seq([seq(T[Bn[i1],Bn[j1]],j1=1..nops(Bn))],i1=1..nops(Bn))]:
end:


#From hw4, Adapted from Auorora Hively's hw4  (See https://sites.math.rutgers.edu/~zeilberg/EM25/hw4posted/hw4AuroraHiveley.txt )
# UM(alpha,beta,gamma,delta):  implements box 1.1. on p. 20 of Nielsen-Chuang 
# note: alpha, beta, delta, gamma are all real valued

UM:=proc(alpha, beta, gamma, delta) local A,X,Y,Z:
A := [[exp(I*alpha),0],[0,exp(I*alpha)]]:
X := [[exp(-I*beta/2),0],[0,exp(I*beta/2)]]:
Y := [[cos(gamma/2), -sin(gamma/2)], [sin(gamma/2), cos(gamma/2)]]:
Z :=[[exp(-I*delta/2),0],[0,exp(I*delta/2)]]:

simplify(evalc(Mul(Mul(A,X),Mul(Y,Z)))):

end:

#ApplyT(T,u) : Inputs a pair T=[Table,V] where V is the domain of the Table Table, describing the action of the linear transformation,let's call it T1,  on the set of
#basis elements V.
#It also inputs u, a linear combination of the members of V, Outputs T1(u), a certain the member of the vector space generated by V
#[Adapted From Joseph Koutsoutis's hw13 (see https://sites.math.rutgers.edu/~zeilberg/EM25/hw13posted/hw13JosephKoutsoutis.txt)]
#T:=UMtoT(C1QG()[7],1,X);
#ApplyT(T,X[[0]]+4*X[[1]]);
ApplyT:=proc(T,u) local T1, V, v, val:
  T1, V := op(T):
  val := 0:
  for v in V do:
    val += coeff(u, v, 1) * T1[v]:
  od:
  expand(simplify(val)):
end:



#RQC(n,r,k): generates a random quantum circuit with n qubits of length r, with gates involving at most k qubits at a time
#obtained by randomly tensor-producting C1QG()[i]. The outputs is a list of length r where S_i is the support of M_i
#[n,[seq([S_i,M_i],i=1..r)]]. Try:
#RQC(4,5,2);
RQC:=proc(n,r,k) local ra,M,k1,S1,M1,i1,i:
ra:=rand(2..7):
M:=[]:
for i from 1 to r do
 k1:=rand(1..k)():
 S1:=convert(randcomb(n,k1),list):
 M1:=C1QG()[ra()]:
 for i1 from 2 to k1 do
  M1:=TP(M1,C1QG()[ra()]):
 od:
M:=[op(M),[S1,M1]]:
od:
[n,M]:
end:



C1QGnames:=proc():[` Identity `, ` PauliX `, `PauliY `, `PauliZ `, ` PhaseS `, `TPi/8` , ` Hadamard `]: end:

#C1QG(): the list of the seven most imporant one-qubit gates (unitary 2 by 2 matrices) in lists-of-lists format
#in that order: Identity `, ` PauliX `, `PauliY `, `PauliZ `, ` PhaseS `, `TPi/8` , ` Hadamard `
#Try:
#C1QG();
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. Try:
#BV(3);
BV:=proc(n) local 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>. Try:
#nBasis(3,X);
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, it also returns the list consisting of the canonical basis
#indexed by X, namely nBasis(n,X); Try:
#M:=C1QG()[7]; UMtoT(M,1,X);
#M:=TP(C1QG()[7],C1QG()[5]); UMtoT(M,2,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(T,n): inputs a pair T=[Table,Basis], where Basis is a list of basis elements generating the underlying vector space, and Table is a table such that
#Table[Basis[i]] is the linear combination of the elements of Basis that describes the action of the operator on the basis element Basis[i]. 
#Outputs the corresponding 2^n by 2^n matrix . Try:
#M:=C1QG()[6]; T:=UMtoT(M,1,X); evalb(TtoUM(T,1)=M);
#M:=TP(C1QG()[6],C1QG()[4]); T:=UMtoT(M,2,X); evalb(TtoUM(T,2)=M);

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:




#Mul(A,B): the product of matrix A and B (assuming that it exists)
Mul:=proc(A,B) local i,j,k:
[seq([seq(expand(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. Try:
#CT([[a11,a12],[a21,a22]]);
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. Try:
#IsUM(UM(Pi/3,Pi/4,Pi/5,Pi/3));
IsUM:=proc(A) local n,i, 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:

#QCtoM(C): inputs  a Quantum circut C=[n,ListOfGates] where there are n qubits and ListOfGates is a list where each entry is a pair [S,Matrix] where S is a subset of
#{1,...n} indicating the active qubits and M is a 2^nops(S) by 2^nops(S) unitary matrix describing the gate's action. Try:
#C:=RQC(4,5,2); QCtoM(C);
QCtoM:=proc(C) local n,M,i,C1:
n:=C[1]:
M:=[seq([0$i,1,0$(2^n-i-1)],i=0..2^n-1)]:
C1:=C[2]:
for i from 1 to nops(C1) do
M:=Mul(Embed1(n,op(C1[i])),M):
od:
M:
end:



#start from C17.txt

#Ord(x,n): inputs a large pos. integer n and x<n m such that gcd(x,n) outputs (by brute force)
#the order of x mod n, i.e. the smallest pos. integer r such that x^r=1. Try:
#Ord(2,101);
Ord:=proc(x,n) local p,r:
if igcd(x,n)<>1 then
 RETURN(FAIL):
fi:
p:=x:
for r from 1 while p mod n<>1 do
p:=p*x mod n:
od:
r:
end:

#FindFact1(n): tries to find a non-trivial factor of n or returns FAIL
FindFact1:=proc(n) local x,r:
x:=rand(2..n-1)():
r:=Ord(x,n):
if r=FAIL then
  RETURN(r):
fi:

if r mod 2=1 then
  RETURN(FAIL):
fi:

if x^(r/2) mod n =n-1  then
  RETURN(FAIL):
fi:
igcd(n,x^(r/2)-1):
end:

#FindFact(n,K): tries to find a factor of n using FindFact1(n) K times or returns FAIL (with prob.
#less than 1/2^K
FindFact:=proc(n,K) local i,p:
for i from 1 to K do
 p:=FindFact1(n):
 if p<>FAIL then
   RETURN(p):
 fi:
od:
FAIL:
end:

#Factorize(n,K): inputs a pos. integer n and outputs the list of prime factors. Try:
#Factorize(1001);
Factorize:=proc(n,K) local p:
if isprime(n) then
  RETURN([n]):
fi:

p:=FindFact(n,K):
if p=FAIL then
 RETURN(FAIL):
fi:
sort([ op(Factorize(p,K)), op(Factorize(n/p,K))    ]):
end:

#ModExp(x,r,n): x^r mod n the fast way where r is given in reverse binary [2^0,2^1,2^2,...]. Try:
#ModExp(2,[1,0,1,0,1],101);
ModExp:=proc(x,r,n) local k,T,p,i:
k:=nops(r)+1:
T[0]:=1:
T[1]:=x:
for i from 1 to k-1 do
 T[i+1]:=T[i]^2 mod n: 
od:
p:=1:
for i from 0 to nops(r)-1 do
if r[i+1]=1 then
  p:=p*T[i+1] mod n:
fi:
od:
p:
end:
#end from C17.txt