% linalg.s % % Purpose: Matrix/vector operations script % Copyright: R.A. van Engelen, Leiden University, 1997banner("Linear Algebra Package").infix ++ is associative(grouping = left).[X] ++ [Y] := [[X, Y]].[X | ColX] ++ [Y | ColY] := [[X, Y] | ColX ++ ColY].[X] ++ [[Y | RowY]] := [[X, Y | RowY]].[X | ColX] ++ [[Y | RowY] | MatY] := [[X, Y | RowY] | ColX ++ MatY].[[X | RowX]] ++ [[Y | RowY]] := [[X | RowX] // [Y | RowY]].([[X | RowX] | MatX] ++ [[Y | RowY] | MatY] :=

[[X | RowX] // [Y | RowY] | MatX ++ MatY])

.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Create zero matrix % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%zero(Rank) := zero(Rank, Rank).zero(RowRank, ColRank) := fill(RowRank, fill(ColRank, 0)).%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Create identity matrix % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ident(1) := [[1]].ident(N) := (iff N <= 0 then error("Cannot create identity matrix of rank " // sprint(N)) else

[[1 | fill(N - 1, 0)] | zero(N - 1, 1) ++ ident(N - 1)]

).%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Matrix transpose % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%prefix transpose is intrinsic.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Matrix-matrix and vector-vector addition % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[[_G2509|_G2510] |_G2513] + [_G2515|_G2516] := error("matrix + vector").[_G2509|_G2510] + [[_G2512|_G2513] |_G2516] := error("vector + matrix").[[_G2509|_G2510] |_G2513] - [_G2515|_G2516] := error("matrix - vector").[_G2509|_G2510] - [[_G2512|_G2513] |_G2516] := error("vector - matrix").%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Dot/inner product % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%infix .* has precedence(*) has latex_name("\cdot").[X] .* [Y] := X * Y.[X | Xs] .* [Y | Ys] := X * Y + Xs .* Ys._G2509.*_G2510:= error("dot product of scalars").[_G2509

|_G2510] .*_G2513:= error("dot product of vector with scalar")._G2512.* [_G2509

|_G2510] := error("dot product of scalar with vector").%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Cross/outer product % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%infix #* has precedence(*) has latex_name("\times").[X] #* [Y] := [X * Y].[X1, X2] #* [Y1, Y2] := [X1 * Y2 - X2 * Y1, X2 * Y1 - X1 * Y2].([X1, X2, X3] #* [Y1, Y2, Y3] :=

[X2 * Y3 - X3 * Y2, X3 * Y1 - X1 * Y3, X1 * Y2 - X2 * Y1])

._G2509#*_G2510:= error("cross product of scalars").[_G2509

|_G2510] #*_G2513:= error("cross product of vector with scalar")._G2512#* [_G2509

|_G2510] := error("cross product of scalar with vector").%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Matrix-matrix and matrix-vector product % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%infix &* is associative(grouping = left) has latex_name("\,").[X | Col] &* [Y] := [X | Col] * Y.[X | Col] &* [[Y | Row]] := matmult_([X | Col], [Y | Row]).[[X | Row]] &* [Y | Col] := [X | Row] .* [Y | Col].([[X | Row1] | Rows1] &* [Y | Row2] :=

matmult_([[X | Row1] | Rows1], [[Y | Row2]]))

.([[X | Row1] | Rows1] &* [[Y | Row2] | Rows2] :=

matmult_([[X | Row1] | Rows1], transpose [[Y | Row2] | Rows2]))

.[_G2509|_G2510] &*_G2513:= error("matrix product with scalar").[ [_G2509|_G2510]

|_G2513] &*_G2516:= error("matrix product with scalar")._G2512&* [_G2509|_G2510] := error("matrix product with scalar")._G2515&* [ [_G2509|_G2510]

|_G2513] := error("matrix product with scalar").% matmult_: matrix product with transpose of matrix(matmult_([[X | Row] | Rows], [[Y | Col]]) :=

[[X | Row] .* [Y | Col] | matmult_(Rows, [[Y | Col]])])

.(matmult_([[X | Row] | Rows], [[Y | Col] | Cols]) :=

[ dotprodrows_([X | Row], [[Y | Col] | Cols])

| matmult_(Rows, [[Y | Col] | Cols]) ])

.matmult_([], [[_G2509|_G2510] |_G2513]) := [].(matmult_([[X] | Rows], [Y | Col]) :=

[X * [Y | Col] | matmult_(Rows, [Y | Col])])

.(matmult_([X | Rows], [Y | Col]) :=

[X * [Y | Col] | matmult_(Rows, [Y | Col])])

.matmult_([], [_G2509|_G2510]) := [].matmult_(_G2509,_G2510) := error("incompatible matrix dimensions").(dotprodrows_([X | Row1], [[Y | Row2] | Rows2]) :=

[[X | Row1] .* [Y | Row2] | dotprodrows_([X | Row1], Rows2)])

.dotprodrows_([_G2509|_G2510], []) := [].%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Gaussian elimination % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%gauss(_G2509) is intrinsic(module = gauss).% Prolog predicate gauss/2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Matrix power % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%([[X | Row] | Mat] ^ -1 :=

(map(drop(Rank, gauss([[X | Row] | Mat] ++ ident(Rank)))) where

Rank = length([[X | Row] | Mat])

))

.[[X | Row] |_G2513] ^ 0 := ident(length([X | Row])).[[X | Row] | Mat] ^ 1 := [[X | Row] | Mat].[[X | Row] | Mat] ^ 2 := [[X | Row] | Mat] &* [[X | Row] | Mat].([[X | Row] | Mat] ^ N :=

(iff op(0, N) == "T" then transpose [[X | Row] | Mat] else

iff N mod 2 == 0 then ( [[X | Row] | Mat]

^ (N div 2)

) ^ 2 else

[[X | Row] | Mat] &* [[X | Row] | Mat] ^ (N - 1)

))

.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Solve (MATLAB left division) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[[X | Row] | Mat] \ [Y | Vec] := [[X | Row] | Mat] ^ -1 &* [Y | Vec].