Louis Lamarche

70 Reputation

6 Badges

14 years, 86 days

MaplePrimes Activity


These are Posts that have been published by Louis Lamarche

This is an application of the previous posts
https://www.mapleprimes.com/posts/209057-Procedure-For-Expanding-Tensor-Product

I have a fourth version of the ExpandQop that will expand automaticaly the power of
quantum tensor product. This is just a minor change to the procedure.

Now here is an application for all this that will help understanding a little about
quantum computing. This is the classical concept of quantum teleportation.

You will need to run the above mentionned file and uncomment the save line in the file
before running the example.

LL
 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument. Manages multiple tensor products   #
#                like A*B*C*F                                        #
#       Version: 3                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

(1)

Setup(unitaryoperators={I,U,X,Y,Z,H,HI,CNOT,CnotI});
Setup(noncommutativeprefix={q,beta,psi});

[unitaryoperators = {CNOT, CnotI, H, HI, I, U, X, Y, Z}]

 

[noncommutativeprefix = {beta, psi, q}]

(2)

Setup(bracketrules= { %Bracket(%Bra(q0), %Ket(q0))=1,
                      %Bracket(%Bra(q1), %Ket(q1))=1,
                      %Bracket(%Bra(q1), %Ket(q0))=0,
                      %Bracket(%Bra(q0), %Ket(q1))=0
                    });

[bracketrules = {%Bracket(%Bra(q0), %Ket(q0)) = 1, %Bracket(%Bra(q0), %Ket(q1)) = 0, %Bracket(%Bra(q1), %Ket(q0)) = 0, %Bracket(%Bra(q1), %Ket(q1)) = 1}]

(3)

####################################################################################
# Load the procedure and set the required global variables
#
read "ExpandQop.m": optp:=op(0,Ket(q0)*Ket(q1)): optpx:= op(0,(Ket(q0)+Ket(q1))^2):
#
####################################################################################

#
# Pauli operators
#
print("Pauli gates");
I:=Ket(q0)*Bra(q0)+Ket(q1)*Bra(q1);        # = sigma[0]
X:=Ket(q1)*Bra(q0)+Ket(q0)*Bra(q1);        # = sigma[1] = sigma[x]
Y:=-i*Ket(q1)*Bra(q0)+i*Ket(q0)*Bra(q1);   # = sigma[2] = sigma[y]
Z:=Ket(q0)*Bra(q0)-Ket(q1)*Bra(q1);        # = sigma[3] = sigma[z]

"Pauli gates"

 

Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q1))

 

Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q1))

 

-I*Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q0))+I*Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q1))

 

Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q0))-Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q1))

(4)

##############################
# Defining the Hadamard gate #
##############################
print("Hadamard gate");
H:= Ket(q0)*Bra(q0)/sqrt(2)+Ket(q0)*Bra(q1)/sqrt(2)+Ket(q1)*Bra(q0)/sqrt(2)-Ket(q1)*Bra(q1)/sqrt(2);

"Hadamard gate"

 

(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q0))+(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q0), Physics:-Bra(q1))+(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q0))-(1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(q1), Physics:-Bra(q1))

(5)

# This is usefull to represent a 2 qubits system
# A more general approach is needed for a n qubit system.
DefineStates:=proc()
    Ket(q00):=Ket(q0)*Ket(q0);  Ket(q01):=Ket(q0)*Ket(q1);
    Ket(q10):=Ket(q1)*Ket(q0);  Ket(q11):=Ket(q1)*Ket(q1);
    Bra(q00):=Dagger(Ket(q00)); Bra(q01):=Dagger(Ket(q01));
    Bra(q10):=Dagger(Ket(q10)); Bra(q11):=Dagger(Ket(q11));
    return;
    end proc:
UndefineStates:=proc()
    Ket(q00):='Ket(q00)'; Ket(q01):='Ket(q01)';
    Ket(q10):='Ket(q10)'; Ket(q11):='Ket(q11)';
    Bra(q00):='Bra(q00)'; Bra(q01):='Bra(q01)';
    Bra(q10):='Bra(q10)'; Bra(q11):='Bra(q11)';
    return;
    end proc:

####################################
# Defining the CNOT gate (2 qubits)
####################################
print("CNOT gate");
CNOT:=Ket(q00)*Bra(q00)+ Ket(q01)*Bra(q01)+ Ket(q11)*Bra(q10)+Ket(q10)*Bra(q11);
DefineStates();
'CNOT'=CNOT;

"CNOT gate"

 

Physics:-`*`(Physics:-Ket(q00), Physics:-Bra(q00))+Physics:-`*`(Physics:-Ket(q01), Physics:-Bra(q01))+Physics:-`*`(Physics:-Ket(q11), Physics:-Bra(q10))+Physics:-`*`(Physics:-Ket(q10), Physics:-Bra(q11))

 

CNOT = Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Bra(q0), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Bra(q1), Physics:-Bra(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Bra(q0), Physics:-Bra(q1))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Bra(q1), Physics:-Bra(q1))

(6)

###########################
# Defining the Bell states
###########################
Ket(beta,x,y)='CNOT.(((H.Ket(x)))*Ket(y))';
Ket(beta00):=CNOT.(Expand((H.Ket(q0)))*Ket(q0));
Ket(beta01):=CNOT.(Expand((H.Ket(q0)))*Ket(q1));
Ket(beta10):=CNOT.(Expand((H.Ket(q1)))*Ket(q0));
Ket(beta11):=CNOT.(Expand((H.Ket(q1)))*Ket(q1));

Physics:-Ket(beta, x, y) = Physics:-`.`(CNOT, Physics:-`*`(Physics:-`.`(H, Physics:-Ket(x)), Physics:-Ket(y)))

 

(1/2)*2^(1/2)*(Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1)))

 

(1/2)*2^(1/2)*(Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0)))

 

-(1/2)*2^(1/2)*(-Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0))+Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1)))

 

(1/2)*2^(1/2)*(Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1))-Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0)))

(7)

##########################################################
# Quantum teleportation
# Reference: Quantum Computation and Quantum Information
#            10th Anniversary Edition
#            Michael A. Nielsen & Isaac L. Chuang
#            Cambridge University Press, Cambridge 2010
#            pp 25-28
##########################################################
print("State to be teleported");
Ket(psi) := a*Ket(q0)+b*Ket(q1);
print("Step 1: Compute the tensor product of the state to be teleported with ", 'Ket(beta00)');
Ket(psi[0])='Ket(psi)'*'Ket(beta00)';
Ket(psi[0]):=Expand(Ket(psi)*Ket(beta00));
print("This is a 3 qubits state");
#######
print("Step 2: Pass these 3 qubits through a  CNOT*I  operator");
'CnotI'='CNOT*I';
CnotI:=ExpandQop(Expand(CNOT*I)):
#
# To see what the CNOTI operator looks like
#
# print("CNOTI=");
# print(op(1,CNOTI)+op(2,CNOTI)+op(3,CNOTI)+op(4,CNOTI));
# print(op(5,CNOTI)+op(6,CNOTI)+op(7,CNOTI)+op(8,CNOTI));
'Ket(psi[1])'='CnotI.Ket(psi[0])';
Ket(psi[1]):=Expand(CnotI.Ket(psi[0]));
#######
print("Step 3: Pass these 3 qubits through an Haldamard*I  operator");
'HalI'='H*I';
HalI:=ExpandQop(Expand(H*I)):
#
# To see what the Haldamard*I operator looks like
#
# print("HalI=");
# print(op(1,HalI)+op(2,HalI)+op(3,HalI)+op(4,HalI));
# print(op(5,HalI)+op(6,HalI)+op(7,HalI)+op(8,HalI));
'Ket(psi[2])'='HalI.Ket(psi[1])';
Ket(psi[2]):=Expand(HalI.Ket(psi[1]));
 

"State to be teleported"

 

a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

 

"Step 1: Compute the tensor product of the state to be teleported with ", Physics:-Ket(beta00)

 

Physics:-Ket(psi[0]) = Physics:-`*`(Physics:-Ket(psi), Physics:-Ket(beta00))

 

(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q1))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q1))

 

"This is a 3 qubits state"

 

"Step 2: Pass these 3 qubits through a  CNOT*I  operator"

 

CnotI = Physics:-`*`(CNOT, I)

 

Physics:-Ket(psi[1]) = Physics:-`.`(CnotI, Physics:-Ket(psi[0]))

 

(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*2^(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q1))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q0))+(1/2)*2^(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q1))

 

"Step 3: Pass these 3 qubits through an Haldamard*I  operator"

 

HalI = Physics:-`*`(H, I)

 

Physics:-Ket(psi[2]) = Physics:-`.`(HalI, Physics:-Ket(psi[1]))

 

(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q1))+(1/2)*b*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q0), Physics:-Ket(q1))+(1/2)*b*Physics:-`*`(Physics:-Ket(q0), Physics:-Ket(q1), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q1))-(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q0), Physics:-Ket(q1))-(1/2)*b*Physics:-`*`(Physics:-Ket(q1), Physics:-Ket(q1), Physics:-Ket(q0))

(8)

UndefineStates();
print("Using contracted names for the first two qubits");
Ket(q00)*Bra(q0)*Bra(q0)='I';
Ket(q01)*Bra(q0)*Bra(q1)='I';
Ket(q10)*Bra(q1)*Bra(q0)='I';
Ket(q11)*Bra(q1)*Bra(q1)='I';
'Ket(psi[2])'=Ket(q00)*Bra(q0)*Bra(q0).Ket(psi[2])+
              Ket(q01)*Bra(q0)*Bra(q1).Ket(psi[2])+
              Ket(q10)*Bra(q1)*Bra(q0).Ket(psi[2])+
              Ket(q11)*Bra(q1)*Bra(q1).Ket(psi[2]);

"Using contracted names for the first two qubits"

 

Physics:-`*`(Physics:-Ket(q00), Physics:-Bra(q0), Physics:-Bra(q0)) = I

 

Physics:-`*`(Physics:-Ket(q01), Physics:-Bra(q0), Physics:-Bra(q1)) = I

 

Physics:-`*`(Physics:-Ket(q10), Physics:-Bra(q1), Physics:-Bra(q0)) = I

 

Physics:-`*`(Physics:-Ket(q11), Physics:-Bra(q1), Physics:-Bra(q1)) = I

 

Physics:-Ket(psi[2]) = (1/2)*a*Physics:-`*`(Physics:-Ket(q00), Physics:-Ket(q0))+(1/2)*b*Physics:-`*`(Physics:-Ket(q00), Physics:-Ket(q1))+(1/2)*a*Physics:-`*`(Physics:-Ket(q01), Physics:-Ket(q0))-(1/2)*b*Physics:-`*`(Physics:-Ket(q01), Physics:-Ket(q1))+(1/2)*a*Physics:-`*`(Physics:-Ket(q10), Physics:-Ket(q1))+(1/2)*b*Physics:-`*`(Physics:-Ket(q10), Physics:-Ket(q0))+(1/2)*a*Physics:-`*`(Physics:-Ket(q11), Physics:-Ket(q1))-(1/2)*b*Physics:-`*`(Physics:-Ket(q11), Physics:-Ket(q0))

(9)

print("Rewriting this result by hand");
'Ket(psi[2])'=(Ket(q00)*(a*Ket(q0)+b*Ket(q1))+
               Ket(q01)*(a*Ket(q0)-b*Ket(q1))+
               Ket(q10)*(a*Ket(q1)+b*Ket(q0))+
               Ket(q11)*(a*Ket(q1)-b*Ket(q0)))/2;

"Rewriting this result by hand"

 

Physics:-Ket(psi[2]) = (1/2)*Physics:-`*`(Physics:-Ket(q00), a*Physics:-Ket(q0)+b*Physics:-Ket(q1))+(1/2)*Physics:-`*`(Physics:-Ket(q01), a*Physics:-Ket(q0)-b*Physics:-Ket(q1))+(1/2)*Physics:-`*`(Physics:-Ket(q10), a*Physics:-Ket(q1)+b*Physics:-Ket(q0))+(1/2)*Physics:-`*`(Physics:-Ket(q11), a*Physics:-Ket(q1)-b*Physics:-Ket(q0))

(10)

DefineStates();
print("If Alice measures 00 Bob does noting");
''I'.   '2*Bra(q00).Ket(psi[2])'' =  I.   2*Bra(q00).Ket(psi[2]);
print("If Alice measures 01 Bob applies the X gate");
''X'.   '2*Bra(q01).Ket(psi[2])'' =  X.   2*Bra(q01).Ket(psi[2]);
print("If Alice measures 10 Bob applies the Z gate");
''Z'.   '2*Bra(q10).Ket(psi[2])'' =  Z.   2*Bra(q10).Ket(psi[2]);
print("If Alice measures 11 Bob applies the X gate and then the Z gate");
''Z'.'X'. '2*Bra(q11).Ket(psi[2])'' =  Z.X. 2*Bra(q11).Ket(psi[2]);

"If Alice measures 00 Bob does noting"

 

Physics:-`.`('I', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q00)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

 

"If Alice measures 01 Bob applies the X gate"

 

Physics:-`.`('X', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q01)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

 

"If Alice measures 10 Bob applies the Z gate"

 

Physics:-`.`('Z', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q10)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

 

"If Alice measures 11 Bob applies the X gate and then the Z gate"

 

Physics:-`.`('Z', 'X', 'Physics:-`.`(Physics:-`*`(2, Physics:-Bra(q11)), Physics:-Ket(psi[2]))') = a*Physics:-Ket(q0)+b*Physics:-Ket(q1)

(11)

 


 

Download QuantumTeleportation.mw

 

 

Version 2 do not enable to expand multiple product like A*A*B*E
Version 3 will now do that.
I just forgot to add this feature.

LL.
 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument. Manages multiple tensor products   #
#                like A*B*C*F                                        #
#       Version: 3                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

(1)

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b,q});

[quantumoperators = {A, B, C, Cn}]

 

[noncommutativeprefix = {a, b, q}]

(2)

opexp:= op(0,10^x):            # exponentiation id
opnp := op(0,10*x):            # normal product id
optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
opdiv:= `Fraction`:            # fraction       id          
opsym:= op(0,x):               # symbol         id
opint:= op(0,10):              # integer        id
opflt:= op(0,10.0):            # float          id
opcpx:= op(0,i):               # complex        id
opbra:= op(0,Bra(q)):          # bra            id
opket:= op(0,Ket(q)):          # ket            id
opmod:= op(0, a mod b):        # mod            id
ExpandQop:=proc(x)
    local nx,ret,j,lkb,cbk,rkb,no,lop,success;
    lop:=op(0,x);
    no:=nops(x);
    if lop = opsym or lop = opint or lop = opflt or
       lop = opbra or lop = opket or lop = opcpx then
         ret:=x;
    else
    if lop = opexp then
        ret:=x;
    else       
    if lop = opnp then
        ret:=1;
        for j from 1 to no do
            ret:=ret*ExpandQop(op(j,x));
        end do;        
    else
    if lop = `+` then
        ret:=0;
        for j from 1 to no do
            ret:=ret+ExpandQop(op(j,x));
        end do;
    else
    if lop = `-` then
        ret:=0;
        for j from 1 to no do
            ret:=ret-ExpandQop(op(j,x));
        end do;
    else
    if lop = opdiv then
       ret:=1;
       for j from 1 to no do
           ret:=ret/ExpandQop(op(j,x));
       end do;
    else
    if lop = opmod then
       ret:=x;
    else
    if lop = optp then
       if (no > 3 ) then
           success:=false;
           nx:=x;
           while not success do
             lkb:=0; cbk:=0; rkb:=0;ret:=1;
             for j from 1 to no do
                 if (j>1) then
                      if(lkb=0) then
                          if( type(op(j-1,nx),Ket) and
                              type(op(j,nx),Bra) ) then lkb:=j-1; fi;
                      else
                          if( type(op(j-1,nx),Ket) and
                              type(op(j,nx),Bra) ) then rkb:=j;   fi;
                      fi;
                      if( type(op(j-1,nx),Bra) and type(op(j,nx),Ket) )
                                                   then cbk:=j;   fi;
                 fi;
             end do;
             if ( (lkb < cbk) and (cbk<rkb) ) then
                 for j from 1     to lkb   do ret := ret*op(j,nx); end do;
                 for j from cbk   to no    do ret := ret*op(j,nx); end do;
                 for j from lkb+1 to cbk-1 do ret := ret*op(j,nx); end do;
             else
               ret:=nx;
             fi;
             
             if nx = ret then
                success := true;
             else
                nx := ret;
             fi
           end do;
       else
           ret:=x;
       fi;
    else ret:=x;
    fi; # optp
    fi; # opmod
    fi; # opdiv
    fi; # `-`
    fi; # `+`
    fi; # `opnp
    fi; # `opexp`
    fi; # opsym, opint, opflt, opbra, opket, opcpx

    return ret;
end proc:

# For saving
# save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"Maple do not compute the tensor product of operators,"

 

"C=A*B gives:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"ExpandQop(C) gives the expected result:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

(3)

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix elements computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix elements computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

 

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

"Matrix elements computed with C appears curious"

 

Physics:-`.`(bras3, C, kets3) = "..."

 

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

 

"Matrix elements computed with Cn as expected"

 

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

(4)

print("Example");
En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

"Example"

 

-(5*I)*2^(1/2)*(-1+x-y-z)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

(5)

print("Another example");
'F'='A*B/sqrt(2)+B*A/sqrt(2)';
F:=A*B/sqrt(2)+B*A/sqrt(2):
'op(1,F)'=op(1,F);
'op(2,F)'=op(2,F);

'Fn'='ExpandQop(F)';
Fn:=ExpandQop(F):
'op(1,Fn)'=op(1,Fn);
'op(2,Fn)'=op(2,Fn);

"Another example"

 

F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

 

op(1, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

Fn = ExpandQop(F)

 

op(1, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

(6)

print("Final example, multiple products");
G:=B*B*B;
'G'=ExpandQop(G);

"Final example, multiple products"

 

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

G = Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

(7)

 


 

Download ExpandQopV3.mw

Here is a major upgrade of the procedure I submitted in february.

https://www.mapleprimes.com/posts/209030-Procedure-For-Computing-The-Tensor-Product

There is a line after the procedure to save it in the file "ExpandQop.m"
In future post I will use it in order to minimize the size of the examples.

Louis Lamarche
 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument                                     #
#       Version: 2                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

(1)

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b,q});

[quantumoperators = {A, B, C, Cn}]

 

[noncommutativeprefix = {a, b, q}]

(2)

opexp:= op(0,10^x):            # exponentiation id
opnp := op(0,10*x):            # normal product id
optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
opdiv:= `Fraction`:            # fraction       id          
opsym:= op(0,x):               # symbol         id
opint:= op(0,10):              # integer        id
opflt:= op(0,10.0):            # float          id
opcpx:= op(0,i):               # complex        id
opbra:= op(0,Bra(q)):          # bra            id
opket:= op(0,Ket(q)):          # ket            id
opmod:= op(0, a mod b):        # mod            id
ExpandQop:=proc(x)
    local ret,j,lkb,cbk,rkb,no,lop;
    lkb:=0; cbk:=0; rkb:=0;
    lop:=op(0,x);
    no:=nops(x);
    if lop = opsym or lop = opint or lop = opflt or
       lop = opbra or lop = opket or lop = opcpx then
         ret:=x;
    else
    if lop = opexp then
        ret:=x;
    else       
    if lop = opnp then
        ret:=1;
        for j from 1 to no do
            ret:=ret*ExpandQop(op(j,x));
        end do;        
    else
    if lop = `+` then
        ret:=0;
        for j from 1 to no do
            ret:=ret+ExpandQop(op(j,x));
        end do;
    else
    if lop = `-` then
        ret:=0;
        for j from 1 to no do
            ret:=ret-ExpandQop(op(j,x));
        end do;
    else
    if lop = opdiv then
       ret:=1;
       for j from 1 to no do
           ret:=ret/ExpandQop(op(j,x));
       end do;
    else
    if lop = opmod then
       ret:=x;
    else
    if lop = optp then
        ret:=1;
       if (no > 3 ) then
           for j from 1 to no do
               if (j>1) then
                    if(lkb=0) then
                        if( type(op(j-1,x),Ket) and
                            type(op(j,x),Bra) ) then lkb:=j-1; fi;
                    else
                        if( type(op(j-1,x),Ket) and
                            type(op(j,x),Bra) ) then rkb:=j;   fi;
                    fi;
                    if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                                then cbk:=j;   fi;
               fi;
           end do;
           if ( (lkb < cbk) and (cbk<rkb) ) then
               for j from 1     to lkb   do ret := ret*op(j,x); end do;
               for j from cbk   to no    do ret := ret*op(j,x); end do;
               for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
           else
               ret:=x;
           fi;
       else
           ret:=x;
       fi;
    else ret:=x;
    fi; # optp
    fi; # opmod
    fi; # opdiv
    fi; # `-`
    fi; # `+`
    fi; # `opnp
    fi; # `opexp`
    fi; # opsym, opint, opflt, opbra, opket, opcpx

    return ret;
end proc:

# For saving
# save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"Maple do not compute the tensor product of operators,"

 

"C=A*B gives:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"ExpandQop(C) gives the expected result:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

(3)

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix elements computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix elements computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

 

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

"Matrix elements computed with C appears curious"

 

Physics:-`.`(bras3, C, kets3) = "..."

 

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

 

"Matrix elements computed with Cn as expected"

 

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

(4)

print("Example");
En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

"Example"

 

-(5*I)*2^(1/2)*(-1+x-y-z)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

(5)

print("Another example");
'F'='A*B/sqrt(2)+B*A/sqrt(2)';
F:=A*B/sqrt(2)+B*A/sqrt(2):
'op(1,F)'=op(1,F);
'op(2,F)'=op(2,F);

'Fn'='ExpandQop(F)';
Fn:=ExpandQop(F):
'op(1,Fn)'=op(1,Fn);
'op(2,Fn)'=op(2,Fn);

"Another example"

 

F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

 

op(1, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

Fn = ExpandQop(F)

 

op(1, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

(6)

 


 

Download ExpandQopV2.mw

I wanted to use MAPLE to preform symbolic quantum computations. The role
of quantum operators and their tensor product is very important in simplying
the understating of such new calculus at least for the beginners. For instance,
(using "o" for the tensor product and "." for the scalar product, H being the Hadamard
operator on a qubit, I the identity operator, and CNOT the 2 qubit controled not
operator)
1) generating the Bells states |Bxy> two stages of operators are needed
     (CNOT) .  (H o  I)  . |x> o |y>

2) performing quantum teleportation of |psi>
     (H o I o I) . (CNOT o I ) . |psi>o |B00>
    followed by a measurements on the first two qubits for driving the application of
    quantum gates to the third qubit.

All these tensor products of operators can be easily written in MAPLE.

Here is a first version of the ExpandQop procedure that will be usefull the purpose of
expanding correctly the tensor product of two quantum operator expressed in Dirac notation.

I hope this is usefull.

LL 

 

######################################################################
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018/02/20                                                 #
#                                                                    #
#         Function name: ExpandQop (x)                               #
#               Purpose: Compute the tensor product of two quantum   #
#                        operators in Dirac notations                #
#              Argument: x - a simple quantum operator               #
#   Future improvements: Manage all +, -, *, /, ^, mod  operations   #
#                        in the argument                             #
#               Version: 1.0                                         #
######################################################################
restart;

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

(1)

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b});

[quantumoperators = {A, B, C, Cn}]

 

[noncommutativeprefix = {a, b}]

(2)

ExpandQop:=proc(x)
    local ret,j,lkb,cbk,rkb,no;
    ret:=1; lkb:=0; cbk:=0; rkb:=0; no:=nops(x);
    if (no > 3 ) then
        for j from 1 to no do
            if (j>1) then
                 if(lkb=0) then
                     if( type(op(j-1,x),Ket) and
                         type(op(j,x),Bra) ) then lkb:=j-1; fi;
                 else
                     if( type(op(j-1,x),Ket) and
                         type(op(j,x),Bra) ) then rkb:=j;   fi;
                 fi;
                 if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                             then cbk:=j;   fi;
            fi;
        end do;
        if ( (lkb < cbk) and (cbk<rkb) ) then
            for j from 1     to lkb   do ret := ret*op(j,x); end do;
            for j from cbk   to no    do ret := ret*op(j,x); end do;
            for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
        else
            ret:=x;
        fi;
    else
        ret:=x;
    fi;
    return ret;
end proc:

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"Maple do not compute the tensor product of operators,"

 

"C=A*B gives:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"ExpandQop(C) gives the expected result:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

(3)

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix element computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix element computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

 

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

"Matrix element computed with C"

 

Physics:-`.`(bras3, C, kets3) = "..."

 

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

 

"Matrix element computed with Cn as expected"

 

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

(4)

 


 

Download ExpandQop.mw

 

 

Page 1 of 1