/********************************************** * tensor.mac * * * * general tensor operations on Maxima * * * * (c) Kenrou Adachi, 2007/1/12- * * Version 2.0 2007/3/7 * * * **********************************************/ load(eigen)$ /* tensortype, tensorrank, zerotensor, tensorp sametensortypep */ 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)) )$ 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 nil )$ zerotensor(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 )$ 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 */ pseudoproduct(a, b) := block( if not listp(first(a)) then maplist(lambda([x], x * b), a) else maplist(lambda([x], pseudoproduct(x, b)), a) )$ tensorproduct(a, b) := block( if tensorp(a) and tensorp(b) then pseudoproduct(a, b) else nil )$ infix("@x@")$ a @x@ b := tensorproduct(a, b)$ /* 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 )$ exchnext(a, i) := block([n], n : tensorrank(a), if (i > n - 1) or (i < 1) then ( print("Index is out of range"), FALSE ) else if i = 1 then exchinext(a) else maplist(lambda([x], exchnext(x, i - 1)), a) )$ exchangetensortype(a, i, j) := block([b, sw], if not tensorp(a) then ( print("Object is not a tensor"), FALSE ) else if i = j then ( print("Invalid indices"), FALSE ) else if i > j then exchengetensortype(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 */ exchttype(a, i, j) := exchangetensortype(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"), FALSE ) 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"), FALSE ) else if (i = j) or (i < 1) or (i > n) or (j < 1) or (j > n) then ( print("Invalid indices"), FALSE ) 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"), FALSE ) 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) )$ /* 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"), FALSE ) else tensorcontract(a @x@ b, i, n + j) )$ /* alias */ ttcontract(a, b, i, j) := tensortensorcontract(a, b, i, j)$ /* end of tensor.mac */