%	linalg.s
%
%	Purpose:	Matrix/vector operations script
%	Copyright:	R.A. van Engelen, Leiden University, 1997

banner("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].