/********************************************** * tensor.mac * * * * general tensor operations on Maxima * * * * (c) Kenrou Adachi, 2007/1/12- * * Version 2.0 2007/3/7 * * Version 2.3 2014/7/12 * * for Maxima-5.9.0, 5.9.1 * * * **********************************************/ load(eigen)$ /***************************************************************/ /* tensortype, tensorrank, zerotensor, tensorp sametensortypep */ /***************************************************************/ /***** prototype of tensortype *******************/ pseudotype(f) := block([ttyp, ttmp], ttmp : copylist(f), ttyp : cons(length(ttmp), []), rectyp(typ, tmp) := block( tmp : first(tmp), if listp(tmp) then ( typ : cons(length(tmp), typ), rectyp(typ, tmp) ) else typ ), reverse(rectyp(ttyp, ttmp)) )$ /*************************************************/ /***** old version zerotensor ********************/ zerotensor2(ttype) := block([n, rtyp, rlst, xero], n : length(ttype), rtyp : reverse(ttype), rlst : first(rtyp), zerolist(m) := block( if m = 0 then [] else cons(0, zerolist(m - 1)) ), rec(m, rz) := block( if m = 0 then [] else cons(rz, rec(m - 1, rz)) ), for i : 1 thru n do if i = 1 then ( xero : zerolist(rlst), rtyp : rest(rtyp) ) else ( rlst : first(rtyp), xero : rec(rlst, xero), rtyp : rest(rtyp) ), xero )$ /****************************************************/ /*** prototype of tensorproduct *********************/ pseudoproduct(a, b) := block( if not listp(first(a)) then maplist(lambda([x], x * b), a) else maplist(lambda([x], pseudoproduct(x, b)), a) )$ /****************************************************/ /***** new version zerotensor ************************/ zerolist(n) := if n = 0 then [] else cons(0, zerolist(n - 1))$ zerotensor(ttyp) := block( if (not listp(ttyp)) or (ttyp = []) then ( print("the argument is not a type of tensor"), [] ) else if rest(ttyp) = [] then zerolist(first(ttyp)) else pseudoproduct(zerolist(first(ttyp)), zerotensor(rest(ttyp))) )$ tensortype(f) := block([typ, za, ta], /* depending on zerotensor below */ typ : pseudotype(f), za : zerotensor(typ), ta : 0 * f, if za = ta then typ else [] )$ tensorp(a) := block([za, ta], za : zerotensor(pseudotype(a)), ta : 0 * a, if za = ta then TRUE else FALSE )$ sametensortypep(a, b) := block([za : 0 * a, zb : 0 * b], if za = zb then TRUE else FALSE )$ tensorrank(a) := length(tensortype(a))$ /*********************************************************/ /******************** tensorproduct **********************/ /*********************************************************/ tensorproduct(a, b) := block( if tensorp(a) and tensorp(b) then pseudoproduct(a, b) else [] )$ infix("@x@")$ a @x@ b := tensorproduct(a, b)$ /***********************************************************/ /******************** contraction **************************/ /***********************************************************/ /* exchangetensortype */ exchinext(a) := block([b, c, ni, nip], ni : length(a), nip : length(a[1]), b : [], for j : 1 thru nip do ( c : [], for k : 1 thru ni do c : endcons(a[k][j], c), b : endcons(c, b) ), b )$ matrixlikep(a) := block( if tensorrank(a) = 2 then TRUE else FALSE )$ matrixtranspose(a) := block( if matrixlikep(a) then exchinext(a) else [] )$ exchnext(a, i) := block([n], n : tensorrank(a), if (i > n - 1) or (i < 1) then ( print("Index is out of range"), [] ) else if i = 1 then exchinext(a) else maplist(lambda([x], exchnext(x, i - 1)), a) )$ tensortranspose(a, i, j) := block([b, sw], if not tensorp(a) then ( print("Object is not a tensor"), [] ) else if i = j then ( print("Invalid indices"), [] ) else if i > j then tensortranspose(a, j, i) else if j = i + 1 then exchnext(a, i) else ( b : copylist(a), for k : i thru j - 1 do b : exchnext(b, k), for k : j - 2 thru i step -1 do b : exchnext(b, k), b ) )$ /* ailias */ exchangetensortype(a, i, j) := tensortranspose(a, i, j)$ exchttype(a, i, j) := tensortranspose(a, i, j)$ /* tensorcontract, contractlast2, tensorvectorcontract */ squarematrixlikep(a) := block( if listp(a[1][1]) then FALSE else if not length(a) = length(a[1]) then FALSE else TRUE )$ matrixtrace(a) := block([n, tra], if not squarematrixlikep(a) then ( print("Indices are not square-matrix-like"), [] ) else ( n : length(a), tra : sum(a[i][i], i, 1, n), tra ) )$ contractlast2(a) := block( if squarematrixlikep(a) then matrixtrace(a) else maplist(lambda([x], contractlast2(x)), a) )$ tensorcontract(a, i, j) := block([n, b], n : tensorrank(a), if n = 0 then ( print("Object is not a tensor"), [] ) else if (i = j) or (i < 1) or (i > n) or (j < 1) or (j > n) then ( print("Invalid indices"), [] ) else if i > j then tensorcontract(a, j, i) else if i = n - 1 then contractlast2(a) else ( b : copylist(a), if not j = n then for p : j thru n - 1 do b : exchnext(b, p), if not i = n - 1 then for p : i thru n - 2 do b : exchnext(b, p), b : contractlast2(b), b ) )$ /* alias */ cntrct(a, i, j) := tensorcontract(a, i, j)$ /* tensortensorcontract, tensorvectorcontract */ simplistp(f) := block( if not listp(first(f)) then TRUE else FALSE )$ /* vector = tensorrank 1 tensor */ /* Warning: the followings don't check Consistency of Indices */ tensorvectorcontract(a, b) := block( if (not tensorrank(b) = 1) or (not tensorp(a)) then ( print("Objects have invalid types"), [] ) /************************************************************* * else if simplistp(a) then * * inprod(a, b) * * else if simplistp(first(a)) then * * maplist(lambda([x], inprod(x, b)), a) * * else * * maplist(lambda([x], tensorvectorcontract(x, b)), a) * *************************************************************/ /* to avoid using inprod */ else if first(reverse(tensortype(a))) = length(b) then contractlast2(tensorproduct(a, b)) else [] )$ /* alias */ tvcontract(a, b) := tensorvectorcontract(a, b)$ infix("@tvc")$ a @tvc b := tensorvectorcontract(a, b)$ tensortensorcontract(a, b, i, j) := block([n], n : tensorrank(a), if (not n = 0) or (not tensorp(b)) then ( print("Objects may not be tensor"), [] ) else tensorcontract(a @x@ b, i, n + j) )$ /* alias */ ttcontract(a, b, i, j) := tensortensorcontract(a, b, i, j)$ /**********************************************************/ /***************** Appendix **********************/ /**********************************************************/ /* simple version of tensortranspose */ ttranspose(a, i, j) := block([b], if i > j then ttranspose(a, j, i) else if j = i+1 then exchnext(a, i) else ( b : copylist(a), for k : i thru j-1 do b : exchnext(b, k), for k : j-2 thru i step -1 do b : exchnext(b, k), b ) )$ /***********************************************************/ /* a list utility from the book ''On Lisp'' by Paul Graham */ /***********************************************************/ flatten(x) := block( rec(xx, acc) := block( if xx = [] then acc else if not listp(xx) then cons(xx, acc) else rec(first(xx), rec(rest(xx), acc)) ), rec(x, []) )$ /**********************************************/ /*** Remark: The newer maxima has "flatten" ***/ /**********************************************/ flatten1(a) := maplist(lambda([x], flatten(x)), a)$ demote(x) := block([tx, tt], tt : [], tx : ttranspose(copylist(x), 2, 3), loop, if not first(tx) = [] then ( tt : append(tt, flatten1(first(tx))), tx : rest(tx) ), if not tx = [] then go(loop), tt )$ pauli(g1, g2) := demote(g1 @x@ g2)$ pauli4(g1, g2, g3, g4) := pauli(pauli(g1, g2), pauli(g3, g4))$ pauli8(g1, g2, g3, g4, g5, g6, g7, g8) := pauli(pauli4(g1, g2, g3, g4), pauli4(g5, g6, g7, g8))$ listtomatrix(a) := block([tv, m, b], if listp(a) then ( b : copylist(a), m : length(a), tv : matrix(first(b)), b : rest(b), for i : 1 thru m - 1 do tv : addrow(tv, b[i]), tv ) else nil )$ matrixtolist(a) := block([tl, m], if matrixp(a) then ( tl : [], m : length(a), for i : m thru 1 step -1 do tl : cons(a[i], tl), tl ) else nil )$ /* alias */ ltom(a) := listtomatrix(a)$ mtol(a) := matrixtolist(a)$ gamma0 : matrix([1,0],[0,1])$ gamma1 : matrix([0,1],[1,0])$ gamma2 : matrix([0,-%i],[%i,0])$ gamma3 : matrix([1,0],[0,-1])$ gm0 : [[1,0],[0,1]]$ gm1 : [[0,1],[1,0]]$ gm2 : [[0,-%i],[%i,0]]$ gm3 : [[1,0],[0,-1]]$ pm1 : pauli(gm1,gm1)$ pm2 : pauli(gm1,gm2)$ pm3 : pauli(gm1,gm3)$ pm0 : pauli(gm0,gm0)$ pdirac1 : ltom(pm1)$ pdirac2 : ltom(pm2)$ pdirac3 : ltom(pm3)$ pdirac0 : ltom(pm0)$ /* end of tensor.mac */