/* QCOMP.MAC VERSION 2.0
 * Copyright (c) Akira SaiToh (2011)
 **************************************************************************
 * qcomp.mac -- A function package for Maxima to simulate quantum computing
 * Version 1.0 was a public-domain file written by Akira SaiToh in 2005.
 * Version 2.0 (this file, updated on 4 July 2011) is an extended version
 * written by Akira SaiToh in 2011, which is released under the GPL ver. 2.
 **************************************************************************
*/
/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
/*
********************************************************************
*This package contains the following functions. Arguments are kept
*unchanged by the functions.
******Note: qubit labels start from 1.*****
*  mtrace(M) -- return the trace of a matrix M.
*  cjt(M) -- return the conjugate transpose of a matrix M.
*  tensorprod2m(A,B) -- return the tensor product (Kronecker product)
*            of matrices A and B.
*  tensorprod([b]) -- return the tensor product of matrices listed in [b].
*            Usage: tensorprod(A, B, C, D);
*  vnentropy(M) -- return the von Neumann entropy of the density matrix
*            M.
*  traceoutAofAB(M,dA,dB) -- for a given dAdB x dAdB matrix M,
*            return a dB x dB (reduced) matrix by tracing-out the
*            subsystem A.
*  traceoutBofAB(M,dA,dB) -- for a given dAdB x dAdB matrix M,
*            return a dA x dA (reduced) matrix by tracing-out the
*            subsystem B.
*  traceoutBofABC(M,dA,dB,dC) -- for a given dAdBdC x dAdBdC matrix M,
*            return a dAdC x dAdC (reduced) matrix by tracing-out the
*            subsystem B.
*  traceoutq(M,i) -- assuming that M is a density matrix of qubits,
*            return the reduced density matrix of the qubits remaining
*            after tracing-out the i-th qubit.
*  traceoutqubits(M,[u]) -- return the reduced density matrix of the
*            qubits remaining after tracing out the qubits labeled by
*            integers found in [u].
*            Usage: traceoutqubits(M, 2, 5, 7);
*  transposeAofAB(M,dA,dB) -- for a given dAdB x dAdB matrix M,
*            return its partial transposition with respect to the
*            subsystem A.
*  transposeBofAB(M,dA,dB) -- for a given dAdB x dAdB matrix M,
*            return its partial transposition with respect to the
*            subsystem B.
*  transposeBofABC(M,dA,dB,dC) -- for a given dAdBdC x dAdBdC matrix
*            M, return its partial transposition with respect to the
*            subsystem B.
*  transposeq(M,i) -- for a given square matrix M, reutrn its partial
*            transposition with respect to the i-th qubit.
*  transposequbits(M,[u]) -- for a given square matrix M, return its
*            partial transposition with respect to the qubits listed in
*            [u].
*            Usage: transposequbits(M, 2, 5, 7);
*
*Example:
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*    load(qcomp);
*    A:matrix([0.7,0],[0,0.3]); B:matrix([0.6,0],[0,0.4]);
*    C:tensorprod(A,B);
*                          [ 0.42   0     0     0   ]
*                          [                        ]
*                          [  0    0.28   0     0   ]
*(%o6)                     [                        ]
*                          [  0     0    0.18   0   ]
*                          [                        ]
*                          [  0     0     0    0.12 ]
*(%i7) traceoutqubits(C,2);
*                                 [ 0.7   0  ]
*(%o7)                            [          ]
*                                 [  0   0.3 ]
*(%i8) traceoutAofAB(C,2,2);
*                                 [ 0.6   0  ]
*(%o8)                            [          ]
*                                 [  0   0.4 ]
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*Plese read the following code carefully before use.
********************************************************************
*/

/* Trace of a square matrix */
mtrace(M) := block([t,i],
	t:0,
	if not (matrixp(M) and length(M)=length(M[1])) then
		error ("Error: the argument is not a square matrix for mtrace().")
	else
	for i:1 thru length(M) do
		t:t+M[i,i],
	t
);

/* Conjugate transpose of a matrix */
cjt(M) := block(
	if not (matrixp(M) or listp(M)) then
		error ("Error: the argument is not a matrix for cjt().")
	else
		conj(transpose(M))
);

/* Trace out the space of a single qubit */
traceoutq(M,i) := block([P,N,Nm,Nn,Nk,Nl,m,n,k,l],
	if not (matrixp(M) and integerp(i)) then
		error("Error: invalid arguments to traceoutq(); this function is to trace out the space of the i-th qubit from a given density matrix M.")
        else if not (length(M)=length(M[1])) then
		error ("Error: the first argument is not a square matrix for traceoutq().")
	else
	(
		N:length(M)/2,
		Nm: 2^(i-1), Nn:N/Nm,
		Nk: Nm, Nl: Nn,
		P:zeromatrix(N,N),
		for m:1 thru Nm do
		for n:1 thru Nn do
		for k:1 thru Nk do
		for l:1 thru Nl do
	P[(m-1)*Nn+n,(k-1)*Nl+l]:M[(m-1)*2*Nn+n,(k-1)*2*Nl+l]+M[(m-1)*2*Nn+Nn+n,(k-1)*2*Nl+Nl+l],
		P
	)
);

/* Trace out the space of specified qubits */
traceoutqubits(M,[u]) := block([P,i,j],
	if not (length(M)=length(M[1])) then
                error ("Error: the first argument is not a square matrix for traceoutqubits().")
	else
	(
		P:M,
		for i:1 thru length(u) do
		(
			P:traceoutq(P,u[i]),
			for j:i+1 thru length(u) do
				if (u[j]>u[i]) then
					u[j]:u[j]-1
				else if (u[j]=u[i]) then
					error("Error: cannot traceout the same space twice.")
		),
		P
	)
);

/* Trace out the space of A of system AB */
traceoutAofAB(M,dA,dB) := block([P,k,l,m],
	if not (matrixp(M) and (integerp(dA) and integerp(dB))) then
		error("Error: invalid arguments to traceoutAofAB().")
	elseif not (length(M) = dA * dB and length(M)=length(M[1])) then
		error("Error: for traceoutAofAB(), M should be a dAdB x dAdB matrix.")
	else
	(
		P:zeromatrix(dB,dB),
		for k:1 thru dA do
		for l:1 thru dB do
		for m:1 thru dB do
	P[l,m]:P[l,m]+M[(k-1)*dB+l,(k-1)*dB+m],
		P
	)
);

/* Trace out the space of B of system AB */
traceoutBofAB(M,dA,dB) := block([P,k,l,m],
	if not (matrixp(M) and (integerp(dA) and integerp(dB))) then
		error("Error: invalid arguments to traceoutBofAB().")
	elseif not (length(M) = dA * dB and length(M)=length(M[1])) then
		error("Error: for traceoutBofAB(), M should be a dAdB x dAdB matrix.")
	else
	(
		P:zeromatrix(dA,dA),
		for k:1 thru dA do
		for l:1 thru dA do
		for m:1 thru dB do
	P[k,l]:P[k,l]+M[(k-1)*dB+m,(l-1)*dB+m],
		P
	)
);

/* Trace out the space of B of system ABC */
traceoutBofABC(M,dA,dB,dC) := block([P,N,i,j,k,l,m],
        if not (matrixp(M) and (integerp(dA) and integerp(dB) and integerp(dC))) then
                error("Error: invalid arguments to traceoutBofABC().")
        elseif not (length(M) = dA * dB * dC and length(M)=length(M[1])) then
                error("Error: for traceoutBofABC(), M should be a dAdBdC x dAdBdC matrix.")
        else
        (
		N:dA * dC,
		P:zeromatrix(N,N),
		for i:1 thru dA do
		for j:1 thru dA do
		for k:1 thru dB do
		for l:1 thru dC do
		for m:1 thru dC do
	P[(i-1)*dC+l,(j-1)*dC+m]:P[(i-1)*dC+l,(j-1)*dC+m]+M[(i-1)*dB*dC+(k-1)*dC+l,(j-1)*dB*dC+(k-1)*dC+m],
		P 
	)
);

/* Partial transpose on A of system AB */
transposeAofAB(M,dA,dB) := block([P,N,i,j,k,l],
	if not (matrixp(M) and (integerp(dA) and integerp(dB))) then
                error("Error: invalid arguments to transposeAofAB().")
        elseif not (length(M) = dA * dB and length(M)=length(M[1])) then
                error("Error: for transposeAofAB(), M should be a dAdB x dAdB matrix.")
        else
	(
		N:dA * dB,
		P:zeromatrix(N,N),
		for i:1 thru dA do
		for j:1 thru dA do
		for k:1 thru dB do
		for l:1 thru dB do
	P[(i-1)*dB+k,(j-1)*dB+l]:M[(j-1)*dB+k,(i-1)*dB+l],
		P
	)
);

/* Partial transpose on B of system AB */
transposeBofAB(M,dA,dB) := block([P,N,i,j,k,l],
	if not (matrixp(M) and (integerp(dA) and integerp(dB))) then
                error("Error: invalid arguments to transposeBofAB().")
        elseif not (length(M) = dA * dB and length(M)=length(M[1])) then
                error("Error: for transposeBofAB(), M should be a dAdB x dAdB matrix.")
        else
	(
		N:dA * dB,
		P:zeromatrix(N,N),
		for i:1 thru dA do
		for j:1 thru dA do
		for k:1 thru dB do
		for l:1 thru dB do
	P[(i-1)*dB+k,(j-1)*dB+l]:M[(i-1)*dB+l,(j-1)*dB+k],
		P
	)
);

/* Partial transpose on B of system ABC */
transposeBofABC(M,dA,dB,dC) := block([P,N,i,j,k,l,m,n],
        if not (matrixp(M) and (integerp(dA) and integerp(dB) and integerp(dC))) then
                error("Error: invalid arguments to transposeBofABC().")
        elseif not (length(M) = dA * dB * dC and length(M)=length(M[1])) then
                error("Error: for transposeBofABC(), M should be a dAdBdC x dAdBdC matrix.")
        else
        (
		N:dA * dB * dC,
		P:zeromatrix(N,N),
		for i:1 thru dA do
		for j:1 thru dA do
		for k:1 thru dB do
		for l:1 thru dB do
		for m:1 thru dC do
		for n:1 thru dC do
	P[(i-1)*dB*dC+(k-1)*dC+m,(j-1)*dB*dC+(l-1)*dC+n]:M[(i-1)*dB*dC+(l-1)*dC+m,(j-1)*dB*dC+(k-1)*dC+n],
		P
	)
);

/* Partial transpose on a single qubit */
transposeq(M,i) := block([W,dA,dB,dC],
	if not (matrixp(M) and integerp(i) and i>0) then
		error("Error: invalid arguments to transposeq().")
	elseif not (length(M)=length(M[1])) then
		error("Error: the first argument is not a square matrix for transposeq().")
	else
	(
		W:M,
		dA:2^(i-1),
		dB:2,
		dC:(length(W)/dA/dB),
		if not (integerp(dC) and dC>0) then
			error("Invalid dimension in transposeq().")
		else
		(
			W:transposeBofABC(W,dA,dB,dC),
			W
		)
	)
);

/* Partial transpose on specified qubits */
transposequbits(M,[u]) := block([W,i],
	if not (matrixp(M) and length(M)=length(M[1])) then
		error("Error: M should be a square matrix for transposequbits().")
	else
	(
		W:M,
		for i:1 thru length(u) do
		(
			if not (integerp(u[i])) then
				error("Error: there is a non-integer specified in transposequbits().")
			else
			(
				W:transposeq(W,u[i])
			)
		),
		W
	)
);

/* Direct product of two matrices */
tensorprod2m(A,B) := block([P,i,j,k,l],
        if not (matrixp(A) and matrixp(B)) then
                error("Error: invalid arguments to tensorprod2m()")
        else
        (
                P:zeromatrix(length(A)*length(B),length(A[1])*length(B[1])),
                for i:1 thru length(A) do
                        for j:1 thru length(A[1]) do
                                for k:1 thru length(B) do
                                        for l:1 thru length(B[1]) do
	P[(i-1)*length(B)+k,(j-1)*length(B[1])+l]:A[i,j]*B[k,l],
                P
        )
);

/* Direct product of matrices in the given order (b1 x b2 x b3 x ...) */
tensorprod([b]) := block([P,i],
	if not (matrixp(b[1])) then
		error("Error: there is a non-matrix entry for tensorprod().")
	else
	(
		P:b[1],
		for i:2 thru length(b) do
		(
			if not (matrixp(b[i])) then
				error("Error: there is a non-matrix entry for tensorprod().")
			else
			(
				P:tensorprod2m(P,b[i])
			)
		),
		P
	)
);

/* Load package `eigen' */
load(eigen);

/* von Neumann entropy (the logarithm is in base 2) */
vnentropy(M) := block([s,L,i],
	if not (matrixp(M) and length(M)=length(M[1])) then
		error("Error: M is not a square matrix for vnenntropy()"),
	
	L:eigenvalues(M), s:0,
	for i:1 thru length(L[1]) do
		if not (L[1][i]=0) then
			s:s+L[2][i]*L[1][i]*log(L[1][i]),/*Note: multiplicity L[2][i] should be taken into account.*/
	s:-s/log(2),
	s
);