/* tensor.v */

rectyp(typ, tmp) := block(
    tmp : first(tmp),
    if listp(tmp) then
    (
        typ : cons(length(tmp), typ),
        rectyp(typ, tmp)
    )
    else
        typ
)$


pseudotype(f) := block([ttyp, ttmp],
    mode_declare(ttyp, list, ttmp, list),
    ttmp : copylist(f),
    ttyp : cons(length(ttmp), []),
    reverse(rectyp(ttyp, ttmp))
)$


pseudoproduct(a, b) := block(
    if not listp(first(a)) then
        maplist(lambda([x], x * b), a)
    else
        maplist(lambda([x], pseudoproduct(x, b)), a)
)$


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],
    mode_declare(typ, list, za, list, ta, list),
    typ : pseudotype(f),
    za : zerotensor(typ),
    ta : 0 * f,
    if za = ta then
        typ
    else
        []
)$


tensorp(a) := block([za, ta],
    mode_declare(za, list, ta, list),
    za : zerotensor(pseudotype(a)),
    ta : 0 * a,
    if za = ta then
        true
    else
        false
)$


sametensortypep(a, b) := block([za : 0 * a, zb : 0 * b],
    mode_declare(za, list, zb, list),
    if za = zb then
        true
    else
        false
)$


tensorrank(a) := length(tensortype(a))$


tensorproduct(a, b) := block(
    mode_declare(b, list),
    if tensorp(a) and tensorp(b) then
        pseudoproduct(a, b)
    else
        []
)$


infix("oxo")$


a oxo b := tensorproduct(a, b)$
 

exchinext(a) := block([b, c, ni, nip],
    mode_declare(b, list, c, list, ni, fixnum, nip, fixnum),
    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
)$


exchnext(a, i) := block([n],
    mode_declare(n, fixnum),
    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],
    mode_declare(b, list),
    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
    )
)$


exchangetensortype(a, i, j) := tensortranspose(a, i, j)$


exchttype(a, i, j) := tensortranspose(a, i, j)$


squarematrixlikep(a) := block(
    if tensorrank(a) # 2 then
        false
    else if length(a) # length(a[1]) then
        false
    else
        true
)$


pseudotrace(a) := block([n, tra],
    mode_declare(n, fixnum, tra, any),
    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
    )
)$


listtomatrix(a) := block([tv, m, b],
    mode_declare(tv, any, m, fixnum, b, list),
    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
        []
)$


matrixtolist(a) := block([tl, m],
    mode_declare(tl, list, m, fixnum),
    if matrixp(a) then
    (
        tl : [],
        m : length(a),
        for i : m thru 1 step -1 do
            tl : cons(a[i], tl),
        tl
    )
    else
        []
)$


ltom(a) := listtomatrix(a)$

mtol(a) := matrixtolist(a)$


tomatrixlast2(a) := block(
    if squarematrixlikep(a) then
        listtomatrix(a)
    else
        maplist(lambda([x], tomatrixlast2(x)), a)
)$


matrixtoscalar(a) := block([n],
    mode_declare(n, fixnum),
    if matrixp(a) then
    (
        n : length(a),
        sum(a[i][i], i, 1, n)
    )
    else
        maplist(lambda([x], matrixtoscalar(x)), a)
)$


contractlast2(a) := matrixtoscalar(tomatrixlast2(a))$


tensorcontract(a, i, j) := block([n, b],
    mode_declare(n, fixnum, b, list),
    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 j # n then
            for p : j thru n - 1 do
                b : exchnext(b, p),
        if i # n - 1 then
            for p : i thru n - 2 do
                b : exchnext(b, p),
        b : contractlast2(b),
        b
    )
)$


cntrct(a, i, j) := tensorcontract(a, i, j)$


simplistp(f) := block(
    if not listp(first(f)) then
        true
    else
        false
)$


tensorvectorcontract(a, b) := block([c : copylist(a)],
    mode_declare(c, list),
    if (tensorrank(b) # 1) or (not tensorp(a)) then
    (
        print("Objects have invalid types"),
        []
    )
    else if first(reverse(tensortype(c))) = length(b) then
        contractlast2(tensorproduct(a, b))
    else
        []
)$


infix("tvc")$

a tvc b := tensorvectorcontract(a, b)$


tensortensorcontract(a, b, i, j) := block([n],
    mode_declare(n, fixnum),
    n : tensorrank(a),
    if (n = 0) or (not tensorp(b)) then
    (
        print("Objects may not be tensor"),
        []
    )
    else
        tensorcontract(tensorproduct(a, b), i, n + j)
)$


ttc(a, b, i, j) := tensortensorcontract(a, b, i, j)$


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))
)$


lflatten(x) := rec(x, [])$


flatten1(a) := maplist(lambda([x], lflatten(x)), a)$


demote(x) := block([tx, tt],
    mode_declare(tx, list, tt, list),
    tt : [],
    tx : tensortranspose(copylist(x), 2, 3),
    loop,
        if first(tx) # [] then
        (
            tt : append(tt, flatten1(first(tx))),
            tx : rest(tx)
        ),
    if tx # [] then go(loop),
    tt
)$


pauli(g1, g2) := demote(tensorproduct(g1, 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))$