/******************************************* * Kenrou's Library for Maxima * * (C) Adachi, Kenrou * * 2004/7/26 version 1.0 * * 2005/10/30 version 1.1 * * 2005/10/31 version 1.1.1 * * 2005/11/1 version 1.1.2 * * 2007/1/12 version 1.2.0 * * 2007/1/31 version 1.2.1 * * 2007/2/3 version 1.2.2 * * 2007/3/8 version 1.3.0 * * 2007/6/30 version 1.3.1 * * 2007/7/1 version 1.3.2 * * 2008/6/5 version 1.3.3 * *******************************************/ load(eigen)$ /*** 基本行列パッケージ ***/ /* _FP_{i,j} : 第i行(列)と第j(列)を入れ換える。 */ _FP(i, j, n) := block([I, A], I : ident(n), A : zeromatrix(n, n), A[i][i] : -1, A[j][j] : -1, A[i][j] : 1, A[j][i] : 1, I + A )$ /* _FP_i(a) : 第i行(列)を a倍する。 */ _FQ(a, i, n) := block([A], A : ident(n), A[i][i] : a, A )$ /* _FP_{i,j}(a) : 第i行(j列)に第j行(i列)の a倍を加える。 */ _FR(a, i, j, n) := block([I,A], I : ident(n), A : zeromatrix(n,n), A[i][j] : a, I + A )$ /*** ブロック行列操作パッケージ ***/ /* ブロック行列の切り出し */ matblock(Mat, i1, i2, j1, j2) := block([A, m, n], m : i2 - i1, n : j2 - j1, A : zeromatrix(m + 1, n + 1), for i : 1 thru m + 1 do for j : 1 thru n + 1 do A[i][j] : Mat[i1 + i - 1][j1 + j - 1], A )$ /* ブロック行列の貼り付け */ matreplace(Mat1, Mat2, i1, i2, j1, j2) := block([A, m, n], m : i2 - i1, n : j2 - j1, A : copymatrix(Mat1), for i : 1 thru m + 1 do for j : 1 thru n + 1 do A[i1 + i - 1][j1 + j - 1] : Mat2[i][j], A )$ /*** 計量パッケージ ***/ /* 計量行列から内積をつくる */ /* G はエルミート行列 */ geom(a, b, G) := block([ta], ta : conjugate(transpose(a)), ratsimp(ta . G . b) )$ /* 絶対値 */ absolute(a, G) := ratsimp(sqrt(geom(a, a, G)))$ /* 角の余弦 */ anglcos(a, b, G) := block([denom], denom : absolute(a, G) * absolute(b, G), ratsimp(geom(a, b, G) / denom) )$ /* Minkowski 行列 */ MinkMatrix : ematrix(4, 4, 2, 1, 1) - ident(4)$ /* Minkowski 内積 */ mink(a, b) := geom(a, b, MinkMatrix)$ /* 時間的ノルム */ tminknorm(a) := sqrt(mink(a,a))$ /* 空間的ノルム */ sminknorm(a) := sqrt(-mink(a, a))$ /* Minkowski ノルム */ minknorm(a) := block( if mink(a, a) < 0 then sqrt(-mink(a, a)) else sqrt(mink(a, a)) )$ /* 双曲的余弦 */ hypangcos(a, b) := ratsimp(mink(a, b) / (tminknorm(a) * tminknorm(b)))$ /*** 使っては駄目ここから ***/ /*** 古い計量パッケージ ***/ /* 一応直した 2007/2/7 */ /* エルミート内積 */ hermitian(x, y) := block([ans, lngth], lngth : length(x), ans : sum(conjugate(transpose(x[i,1])) * y[i,1], i, 1, lngth), ans )$ /* ユークリッド内積 */ riemannian(x,y) := block([ans, lngth, i], lngth : length(x), ans : sum(transpose(x[i,1]) * y[i,1], i, 1, lngth), ans )$ eucliddian(x, y) := riemannian(x, y)$ /* ローレンツ空間のベクトルの時間成分 */ ltfctr(x) := x[1][1]$ /* 空間成分 */ spfctr(x) := block([ans, slngth], slength : length(x) - 1, ans : zeromatrix(slngth, 1), for i : 1 thru slength do ans[i][1] : x[i + 1][1], ans )$ /* ローレンツ内積 */ lorentzian(x, y) := block([ans, tx, ty, sx, sy], tx : ltfctr(x), ty : ltfctr(y), sx : spfctr(x), sy : spfctr(y), ans : tx * ty - inprod(sx, sy), ans )$ /* ミンコフスキノルム */ minkowski(x) := block([tau], tau : lorentzian(x, y), if tau >= 0 then sqrt(tau) else sqrt(-tau) )$ /* 3元速度 */ newton_velocity(x, t) := diff(x, t)$ /* 固有時 */ /* 一般には計算できないかも? */ tvctr(x, t) := diff(x, t)$ proper_time(x, t, alpha, beta) := block([tau], /* minkowski(diff(x, t)) > 0 for all t */ tau : integrate(sqrt(lorentzian(tvctr(x,t), tvctr(x,t))), t, alpha, beta), tau )$ /* 使っては駄目ここまで */ /*** Gram-Schmidt 直交化 パッケージ ***/ /* Gram-Schmidt's Method A : (m,n)-matrix (m >= n), gso(A) : orthogonalizes A, gson(A) : orthonormalizes A */ /* Gram-Schmidt 直交化 */ gso(A) := block([m, n, B, C, D], m : length(A), B : transpose(A), n : length(B), C : gramschmidt(B), D : zeromatrix(n, m), for i:1 thru n do D[i] : C[i], transpose(D) )$ /* Gram-Schmidt 正規直交化 */ gson(A) := block([m, n, B, C, D], m : length(A), B : transpose(A), n : length(B), C : gramschmidt(B), D : zeromatrix(n, m), for i : 1 thru n do D[i] : unitvector(C[i]), transpose(D) )$ /*** 固有ベクトルパッケージ ***/ /* 固有ベクトル --> 行列 */ eigmat(A) := block([n, B, C], n : length(A), B : eigenvectors(A), C : zeromatrix(n, n), for i : 1 thru n do C[i] : B[i + 1], transpose(C) )$ /* (対角化)ユニタリ行列 */ dzmat(A) := ratsimp(gson(eigmat(A)))$ /* 対角行列にする */ diagonalize(A) := block([B, C], B : dzmat(A), C : B^^(-1) . A . B, ratsimp(C) )$ /*** ベクトル解析パッケージ ***/ /* 前置作用素宣言 */ prefix("Grad", 142, expr, expr)$ prefix("Div", 142, expr, expr)$ prefix("Rot", 142 , expr, expr)$ /* Gradient */ "Grad"(f) := block([gv], gv : zeromatrix(3,1), gv[1][1] : diff(f, x), gv[2][1] : diff(f, y), gv[3][1] : diff(f, z), return(gv) )$ /* Divergence */ "Div"(v) := block([f], f : diff(v[1][1], x) + diff(v[2][1], y) + diff(v[3][1], z), return(f) )$ /* Rotation */ "Rot"(v) := block([w], w : zeromatrix(3,1), w[1][1] : diff(v[3][1], y) - diff(v[2][1], z), w[2][1] : diff(v[1][1], z) - diff(v[3][1], x), w[3][1] : diff(v[2][1], x) - diff(v[1][1], y), return(w) )$ /* Laplacian */ Laplacian(f) := div (grad (f))$ /** (3次元)内積と外積 **/ /* 中置演算子宣言 */ infix("@*")$ infix("@x")$ /* "@*" : 内積 */ v @* w := block([val], val : v[1][1] * w[1][1] + v[2][1] * w[2][1] + v[3][1] * w[3][1] , return(val) )$ /* "@x" : 外積 */ v @x w := block([valvec], valvec : zeromatrix(3,1), valvec[1][1] : v[2][1] * w[3][1] - v[3][1] * w[2][1] , valvec[2][1] : v[3][1] * w[1][1] - v[1][1] * w[3][1] , valvec[3][1] : v[1][1] * w[2][1] - v[2][1] * w[1][1] , return(valvec) )$ /*** ユーティリティ ***/ /* 列ベクトルの生成 */ colvector(A) := transpose(matrix(A))$ /* 行ベクトルの生成 */ rowvector(A) := matrix(A)$ /*** 行列の基本変形パッケージ ***/ /** コピーをつくって変形 **/ /* Aの第i行をk倍する */ cp_col_sp(A, k, i) := block([B], B : copymatrix(A), B[i] : k * B[i], B )$ /* Aの第i行に第j行のk倍を加える */ cp_col_ad(A, k, i, j) := block([B], B : copymatrix(A), B[i] : B[i] + k * B[j], B )$ /* Aの第i行と第j行を入れ換える */ cp_col_ex(A, i, j) := block([B, dummy], B : copymatrix(A), dummy : B[i], B[i] : B[j], B[j] : dummy, B )$ /* Aの第i列をk倍する */ cp_row_sp(A, k, i) := block([B], B : cp_col_sp(transpose(A), k, i), transpose(B) )$ /* Aの第i列に第j列のk倍を加える */ cp_row_ad(A, k, i, j) := block([B], B : cp_col_ad(transpose(A), k, i, j), transpose(B) )$ /* Aの第i列と第j列を入れ換える */ cp_row_ex(A, i, j) := block([B], B : cp_col_ex(transpose(A), i, j), transpose(B) )$ /** 破壊的変形 **/ /* Aの第i行をk倍する */ col_sp(A, k, i) := block( A[i] : k * A[i], A )$ /* Aの第i行に第j行のk倍を加える */ col_ad(A, k, i, j) := block( A[i] : A[i] + k * A[j], A )$ /* Aの第i行と第j行を入れ換える */ col_ex(A, i, j) := block([dummy], dummy : A[i], A[i] : A[j], A[j] : dummy, A )$ /* Aの第i列をk倍する */ row_sp(A, k, i) := block([len], len : length(transpose(A)[1]), for p : 1 thru len do A[p][i] : k * A[p][i], A )$ /* Aの第i列に第j列のk倍を加える */ row_ad(A, k, i, j) := block([len], len : length(transpose(A)[1]), for p : 1 thru len do A[p][i] : A[p][i] + k * A[p][j], A )$ /* Aの第i列と第j列を入れ換える */ row_ex(A, i, j) := block([len, c], len : length(transpose(A)[1]), for p : 1 thru len do ( c : A[p][i], A[p][i] : A[p][j], A[p][j] : c ), A )$ /* 途中の変形を表示して階段行列に変形 */ /* あまり賢くないが教育的なアルゴリズムによる */ /* 破壊的変形 */ kaidan(A) := block([m, n, i, j, ir], if matrixp(A) = FALSE then ( print(A," is not a matrix"), return(FALSE) ), m : length(transpose(A)[1]), n : length(A[1]), findnonzero(A, pi, pj, mm) := block([pix], pix : pi, jump, if not A[pix, pj] = 0 then pix else if pix = mm then (-1) else ( pix : pix + 1, go(jump) ) ), i : 1, j : 1, loop, ir : findnonzero(A, i, j, m), if ir = -1 then ( if j = n then go(finish) else ( j : j + 1, go(loop) ) ) else if ir > i then ( col_ex(A, i, ir), print(A) ), col_sp(A, 1 / A[i, j], i), print(A), if i > 1 then for k : 1 thru i - 1 do if not A[k, j] = 0 then ( col_ad(A, -A[k, j], k, i), print(A) ), if i < m then for k : i + 1 thru m do if not A[k, j] = 0 then ( col_ad(A, -A[k, j], k, i), print(A) ), if i < m and j < n then ( i : i + 1, j : j + 1, go(loop) ), finish, A )$ /* コピーを作って変形 */ cp_kaidan(A) := block([CPA], CPA : copymatrix(A), kaidan(CPA) )$ /* 階段行列を標準形に */ /* A が階段行列でなければ意味がない */ /* 階段行列を標準形に */ /* A が階段行列でなければ意味がない */ kai2hyo(A) := block([m, n, minima], if matrixp(A) = FALSE then ( print(A, " is not a matrix!"), return(FALSE) ), n : length(A[1]), m : length(transpose(A)[1]), if m > n then minima : n else minima : m, findfirst(A, ii, nn) := block([ij, l], for l : ii thru nn do if (A[ii, l] = 1) then ij : l else ij : l + 1, if ij <= nn then return(ij) else return(-1) ), block([i, j, k], for i : 1 thru minima do ( j : findfirst(A, i, n), if (not j = -1) and (not i = j) then ( row_ex(A, i, j), print(A) ), if i < n then for k : i + 1 thru n do if not A[i, k] = 0 then ( row_ad(A, -A[i, k], k, i), print(A) ) ) ), A )$ /* 階段行列をコピーをつくってから標準形に */ cp_kai2hyo(A) := block([CPA], CPA : copymatrix(A), kai2hyo(CPA) )$ /* 行列を標準形に */ hyoujyun(A) := kai2hyo(kaidan(A))$ /* コピーをつくって標準形に */ cp_hyoujyun(A) := kai2hyo(cp_kaidan(A))$ /* 正則行列 A の逆行列を計算過程を表示して求める。 */ myinv(A) := block([I, AI, P, m, n], m : length(transpose(A)[1]), n : length(A[1]), if not m = n then ( print(A, " is not a square matrix"), return(FALSE) ), if not m = rank(A) then ( print(A, " is not a non-singular matrix"), return(FALSE) ), I : copymatrix(ident(m)), AI : addcol(A, I), kaidan(AI), P : matblock(AI, 1, m, n + 1, n + m), P )$ /* 行列 A を階段行列に変形する正則行列 P を求める */ queryp(A) := block([P, CPA, m, n, i, j, ir, dummy, c], if matrixp(A) = FALSE then ( print(A," is not a matrix!"), return(FALSE) ), m : length(transpose(A)[1]), n : length(A[1]), findnonzero(A, pi, pj, mm) := block([pix], pix : pi, jump, if not A[pix, pj] = 0 then pix else if pix = mm then (-1) else ( pix : pix + 1, go(jump) ) ), P : copymatrix(ident(m)), CPA : copymatrix(A), i : 1, j : 1, loop, ir : findnonzero(CPA, i, j, m), if ir = -1 then ( if j = n then go(finish) else ( j : j + 1, go(loop) ) ) else if ir > i then ( dummy : CPA[_i], CPA[i] : CPA[ir], CPA[ir] : dummy, dummy : P[i], P[i] : P[ir], P[ir] : dummy, print(CPA, " ", P) ), c : 1 / CPA[i, j], CPA[i] : c * CPA[i], P[i] : c * P[i], print(CPA, " ", P), if i > 1 then for k : 1 thru i - 1 do if not CPA[k, j] = 0 then ( c : CPA[k, j], CPA[k] : CPA[k] - c * CPA[i], P[k] : P[k] - c * P[i], print(CPA, " ", P) ), if i < m then for k : i + 1 thru m do if not CPA[k, j] = 0 then ( c : CPA[k, j], CPA[k] : CPA[k] - c * CPA[i], P[k] : P[k] - c * P[i], print(CPA, " ", P) ), if i < m and j < n then ( i : i + 1, j : j + 1, go(loop) ), finish, P )$ /* 途中の結果を印字せずに階段行列を求める */ quiet_kaidan(A) := block([m, n, i, j, ir], if matrixp(A) = FALSE then ( print(A," is not a matrix"), return(FALSE) ), m : length(transpose(A)[1]), n : length(A[1]), findnonzero(A, pi, pj, mm) := block([pix], pix : pi, jump, if not A[pix, pj] = 0 then pix else if pix = mm then (-1) else ( pix : pix + 1, go(jump) ) ), i : 1, j : 1, loop, ir : findnonzero(A, i, j, m), if ir = -1 then ( if j = n then go(finish) else ( j : j + 1, go(loop) ) ) else if ir > i then ( col_ex(A, i, ir) ), col_sp(A, 1 / A[i, j], i), if i > 1 then for k : 1 thru i - 1 do if not A[k, j] = 0 then ( col_ad(A, -A[k, j], k, i) ), if i < m then for k : i + 1 thru m do if not A[k, j] = 0 then ( col_ad(A, -A[k, j], k, i) ), if i < m and j < n then ( i : i + 1, j : j + 1, go(loop) ), finish, A )$ /* 階段行列を標準形に変形する正則行列 Q を求める。 */ queryq(A) := block([m, n, minima, Q, CPA], if matrixp(A) = FALSE then ( print(A, " is not a matrix!"), return(FALSE) ), n : length(A[1]), m : length(transpose(A)[1]), Q : copymatrix(ident(n)), CPA : copymatrix(A), if m > n then minima : n else minima : m, findfirst(A, ii, nn) := block([ij, l], for l : ii thru nn do if (A[ii, l] = 1) then ij : l else ij : l + 1, if ij <= nn then return(ij) else return(-1) ), block([i, j, k], for i : 1 thru minima do ( j : findfirst(CPA, i, n), if (not j = -1) and (not i = j) then ( row_ex(CPA, i, j), row_ex(Q, i, j), print(Q) ), if i < n then for k : i + 1 thru n do if not CPA[i, k] = 0 then ( row_ad(Q, -CPA[i, k], k, i), row_ad(CPA, -CPA[i, k], k, i), print(Q) ) ) ), Q )$ getpq(A) := block([P, Q, B, C], B : copymatrix(A), C : quiet_kaidan(B), P : queryp(A), Q : queryq(C), [P, Q] )$ /* 2008/6/5 追加 */ /* 正方行列を上三角行列に基本変形 */ /* ある行に他のある行の定数倍を加えるという操作だけを使う。*/ kamisankaku(A) := block([m, n, i, ir], if matrixp(A) = FALSE then ( print(A," is not a matrix"), return(FALSE) ), m : length(transpose(A)[1]), n : length(A[1]), findnonzero(A, pi, pj, mm) := block([pix], pix : pi, jump, if not A[pix, pj] = 0 then pix else if pix = mm then (-1) else ( pix : pix + 1, go(jump) ) ), if not m = n then ( print(A, " is not a square-matrix"), return(FALSE) ), i : 1, loop, if i = m then go(finish), if not A[i, i] = 0 then ( for k : i + 1 thru m do if not A[k, i] = 0 then ( col_ad(A, -(A[k, i]/A[i, i]), k, i), print(A) ) ) else ( ir : findnonzero(A, i, i, m), if ir = -1 then ( i : i + 1, go(loop) ) else ( col_ad(A, 1, i, ir), for k : i + 1 thru m do if not A[k, i] = 0 then ( col_ad(A, -(A[k, i]/A[i, i]), k, i), print(A) ) ) ), if i < m then ( i : i + 1, go(loop) ), finish, A )$ /* 正方行列の行列式を求める。 */ /* 上三角行列に変形して、対角成分の積を取る。*/ gyouretsushiki(A) := block([B, detm, m, n], B : copymatrix(A), m : length(transpose(B)[1]), n : length(B[1]), if not m = n then ( print(A, " is not a square-matrix"), return(FALSE) ), kamisankaku(B), detm : prod(B[k, k], k, 1, m), detm )$ /*** 組み立て除法パッケージ ***/ /* ;; 組み立て除法 ;; 多項式はリストで表現する、 例えば ;; x^3 + 3*x^2 - 5 なら [1, 3, 0, -5] のように表現する。 ;; 商と剰余をリストで返す、例えば ;; [1, 0, 1, -2] なら、商が x^2 + 1 で剰余が -2 であることを示す。 ;; kumitate(p, a) : 多項式 p(x) を x - a で割る。 ;; (akumitate(p, a) も同じ) ;; ekumitate(p d) : d は a*x + b の表現 [a, b] の形式 ;; (aekumitate(p, d), yaekumitate(p, d) も同じ) */ /* 組み立て除法 */ kumitate(p, a) := block( kurikaesi(ss, pp, aa) := block( if pp = [] then ss else kurikaesi(cons(aa * first(ss) + first(pp), ss), rest(pp), aa) ), if rest(p) = [] then cons(0, p) else reverse(kurikaesi(cons(first(p), []), rest(p), a)) )$ /* もう一つの組み立て除法 */ akumitate(p, a) := block( kurikaesi(ss, pp, aa) := block( if pp = [] then ss else kurikaesi(cons(aa * first(ss) + first(pp), ss), rest(pp), aa) ), if rest(p) = [] then cons(0, p) else block([s], s : [0], rest(reverse(kurikaesi(s, p, a))) ) )$ /* 非モニック多項式の組み立て除法 */ ekumitate(p, d) := block([a, b, tp, r], a : first(d), b : first(rest(d)), tp : reverse(kumitate(p, -b/a)), r : first(tp), reverse(cons(r, map(lambda([x], x/a), rest(tp)))) )$ aekumitate(p, d) := block([a, b, tp, r], a : first(d), b : second(d), tp : reverse(kumitate(p, -b/a)), r : first(tp), reverse(append(cons(r, []), map(lambda([x], x/a), rest(tp)))) )$ yaekumitate(p, d) := block([a, b, tp, r], a : first(d), b : second(d), tp : map(lambda([x], x/a), kumitate(p, -b/a)), r : last(tp), r : r * a, tp )$ /*** シンプソン法による数値積分 ***/ /************************************************* * simpson{2|3}(f, a, b, n) = \int_a^b f(x) dx * * n : odd * * * *************************************************/ /***************** example *********************** (%i2) fun(x) := 4 / (1 + x^2); 4 (%o2) fun(x) := ------ 2 1 + x (%i3) simpson3(fun, 0.0, 1.0, 1001); (%o3) 3.141592653589795 **************************************************/ /* 繰り返しによる */ simpson2(f, a, b, n) := block([sum, h, m, x], h : (b - a) / (n - 1), m : (n - 1) / 2, sum : 4 * f(h), for i:2 thru m do ( x : h * (2 * i - 1), sum : sum + 4 * f(x) + 2 * f(x - h) ), sum : sum + f(a) + f(b), sum : sum * h / 3, return(sum) )$ /* 再帰による */ simpson3(f, a, b, n) := block([h], h : (b - a) / (n - 1), dlt(i) := 4 * f(a + h * (2 * i - 1)) + 2 * f(a + h * (2 * i)), rsum(_dlt, _a, _inc, _b) := if _a > _b then 0 else _dlt(_a) + rsum(_dlt, _inc(_a), _inc, _b), (f(a) - f(b) + rsum(dlt, 1, lambda([_i], _i + 1), (n - 1) / 2) ) * h / 3 )$ /* 汎用テンソル演算パッケージ */ /* 2007/3/7 */ /* 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 kenrou.mac */