/* 汎用テンソル演算 の 実験的パッケージ */
/* (C) Kenrou Adachi */
/* 2007/01/12 */
/* 2007/01/13 */
/* 2007/01/27 */
/* 2007/02/09 */
/* 2007/03/07 */
/* version 1.3 */

load(eigen)$

/* テンソルの型 */
/*
 \mathbb{R}^{r1} \otimes \mathbb{R}^{r2} \otimes \cdots \otimes \mathbb{R}^{rk}
 に属するテンソルの型をリスト [r1,r2,...,rk] で表現する。
*/

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


/* テンソルの階数 */

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


/* 与えられた型のゼロテンソルを生成する。 */

zerotensor(ttype) := block([rr, rk, rst, rtyp, tt1, tt2],
    rr : length(ttype),
    tt1 : [],
    rtyp : reverse(ttype),
    rk : first(rtyp),
    rst : rest(rtyp),

    for i : 1 thru rk do
        tt1 : cons(0, tt1),

    if not rst = [] then     /* if rr > 1 then */
    (
        rk : first(rst),
        for j : 2 thru rr do
        (
            tt2 : [],
            for i : 1 thru rk do
                tt2 : cons(tt1, tt2),
            tt1 : tt2,
            rst : rest(rst),         /* if not j = rr then   */
            if not rst = [] then     /* (  rst : rest(rst),  */
                rk : first(rst)      /*    rk  : first(rst)) */
        )
    ),
    tt1
)$


/* もう一つのゼロテンソル */
/* 2007/2/9 */

zerotensor2(ttype) := block([n, rtyp, rlst, zerot],
    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
        (
            zerot : zerolist(rlst),
            rtyp : rest(rtyp)
        )
        else
        (
            rlst : first(rtyp),
            zerot : rec(rlst, zerot),
            rtyp : rest(rtyp)
        ),

    zerot
)$


/* テンソル積 */

infix("@x@")$


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


a @x@ b := rp(a, b)$


/* 縮約 */

/* テンソルとベクトルの縮約 */

simplistp(f) := block(
    if not listp(first(f)) then
        TRUE
    else
        FALSE
)$


tvcontract(a, b) := block(
    if simplistp(a) then
        inprod(a, b)
    else if simplistp(first(a)) then
        maplist(lambda([x], inprod(x, b)), a)
    else
        maplist(lambda([x], tvcontract(x, b)), a)
)$


infix("@tvc")$

a @tvc b := tvcontract(a, b)$


/* 2007/1/29 */
/* テンソルの型を変換 */

/* [r1,r2,...] ---> [r2,r1,...] */

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


/* [..., ri, r(i+1), ...] ---> [..., r(i+1), ri, ...] */

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

/* [...,ri,...,rj,...] ---> [...,rj,...,ri,...] */


exchttype(a, i, j) := block([b],
    if i = j then
    (
        print("Invalid indices"),
        FALSE
    )
    else if i > j then
        exchttype(a, j, i)
    else if j = i + 1 then
        exchnext(a, i)
    else
    (
        b : copylist(a),
        for k : i thru j - 1 do         /* ---> [...,r(i+1),...,rj,ri,...] */
            b : exchnext(b, k),
        for k : j - 2 thru i step -1 do /* ---> [...,rj,r(i+1),...,ri,...] */
            b : exchnext(b, k),
        b
    )
)$


/* 縮約 */

/* 添え字が正方行列状か? */

simpsqmatp(a) := block(
    if listp(a[1][1]) then
        FALSE
    else if not length(a) = length(a[1]) then
        FALSE
    else
        TRUE
)$

/* 正方行列状 2階テンソルのトレース */

cntrmat(a) := block([n, tra],
    if not simpsqmatp(a) then
    (
        print("The indices are not square-matrix-like"),
        FALSE
    )
    else
    (
        n : length(a),
        tra : sum(a[i][i], i, 1, n),
        tra
    )
)$

/* 最後の2つの添字について縮約 */

cntrlst2(a) := block(
    if simpsqmatp(a) then
        cntrmat(a)
    else
        maplist(lambda([x], cntrlst2(x)), a)
)$

/* 任意の２つの添字について縮約 */

cntrct(a, i, j) := block([n, b],
    n : tensorrank(a),
    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
        cntrct(a, j, i)
    else if (i = n - 1) then
        cntrlst2(a)
    else
    (
        b : copylist(a),

/* bug fixed 2007/1/31 */
    /* [...,ri,...,rj,...,r(n-1),rn] ---> [...,ri,...,#,...,r(n-1),rn,rj] */

        if not j = n then
            for p : j thru n - 1 do
                b : exchnext(b, p),

    /* [...,ri,...,#,...r(n-1),rn,rj] ---> [...,#,...,#,...,r(n-1),rn,ri,rj] */

        if not i = n - 1 then
            for p : i thru n - 2 do
                b : exchnext(b, p),

    /* ---> [...,#,...,#,...,r(n-1),rn] */

        b : cntrlst2(b),
        b
    )
)$


/* テンソル積を取って縮約 */
/* 2007/2/3 */

ttcntrct(a, b, i, j) := block([n],
    n : tensorrank(a),
    cntrct(a @x@ b, i, n + j)
)$

/* テンソルかどうか判定 */
/* 2007/3/7 */


sametensortypep(a, b) := block([ta : 0 * a, tb : 0 * b],
    if ta = tb then
        TRUE
    else
        FALSE
)$


tensorp(a) := block([za : zerotensor(tensortype(a)), ta : 0 * a],
    if za = ta then
        TRUE
    else
        FALSE
)$