%	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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[[_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].