/*A function package for Maxima to simulate a quantum computing*/
/*This is a public-domain file written by Akira SaiToh in 2005.*/
/*
********************************************************************
*This package contains the following functions:
*  mtrace
*  traceoutq
*  traceoutqubits
*  tensorprod2m
*  tensorprod
*  vnentropy
*Among them, the following functions should be used mainly:
*  mtrace, the trace of a matrix A
*          Usage: mtrace(A);
*  traceoutqubits, trace out specified qubits
*          Usage: traceoutqubits(A, 2, 3, 4);
*  tensorprod, tensor product of matrices
*          Usage: tensorprod(A, B, C, D);
*  vnentropy, the von Neumann entropy (base 2) of a matrix A
*          Usage: vnentropy(A);
*
*Note: qubit labels start from 1.
*
*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 ]
*    traceoutqubits(C,2);
*                                 [ 0.7   0  ]
*(%o7)                            [          ]
*                                 [  0   0.3 ]
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*Plese read the following code carefully before use.
********************************************************************
*/

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


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

/* Trace out the space of specified qubits */
traceoutqubits(A,[u]) := ([P],
	P:A,
	for i 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("Cannot traceout the same space twice.")
	),
	P
);


/* Direct product of two matrices */
tensorprod2m(A,B):=([P],
        if not (matrixp(A) and matrixp(B)) then
                error("Invalid arguments to tensorprod()")
        else
        (
                P:zeromatrix(length(A)*length(B),length(A[1])*length(B[1])),
                for i thru length(A) do
                        for j thru length(A[1]) do
                                for k thru length(B) do
                                        for l 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]):=([P],
	P:b[1],
	for i:2 thru length(b) do
		P:tensorprod2m(P,b[i]),
	P
);

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

/* von Neumann entropy (the logarithm is in base 2) */
vnentropy(A):=([s],
	if not (matrixp(A) and length(A)=length(A[1])) then
		error("A is not a square matrix for vnenntropy()"),
	
	L:eigenvalues(A), s:0,
	for i thru length(L[1]) do
		if not (L[1][i]=0) then
			s:s+L[2][i]*L[1][i]*log(L[1][i]),
	s:-s/log(2),
	s
);