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

[[_G207698 | _G207699] | _G207702] + [_G207704 | _G207705] := error(
"matrix + vector").
 [_G207698 | _G207699] + [[_G207701 | _G207702] | _G207705] := error(
"vector + matrix").
[[_G207698 | _G207699] | _G207702] - [_G207704 | _G207705] := error(
"matrix - vector").
 [_G207698 | _G207699] - [[_G207701 | _G207702] | _G207705] := 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.

_G207698 .* _G207699 := error("dot product of scalars").
[ _G207698
| _G207699 ] .* _G207702 := error("dot product of vector with scalar").
_G207701 .* [ _G207698
            | _G207699 ] := 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])
.

_G207698 #* _G207699 := error("cross product of scalars").
[ _G207698
| _G207699 ] #* _G207702 := error("cross product of vector with scalar").
_G207701 #* [ _G207698
            | _G207699 ] := 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]))
.

 [_G207698 | _G207699] &* _G207702 := error("matrix product with scalar").
[ [ _G207698 | _G207699 ]
| _G207702 ] &* _G207705 := error("matrix product with scalar").
_G207701 &* [_G207698 | _G207699] := error("matrix product with scalar").
_G207704 &* [ [ _G207698 | _G207699 ]
            | _G207702 ] := 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_([], [[_G207698 | _G207699] | _G207702]) := [].
(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_([], [_G207698 | _G207699]) := [].
matmult_(_G207698, _G207699) := error("incompatible matrix dimensions").

(dotprodrows_([X | Row1], [[Y | Row2] | Rows2]) :=
 [[X | Row1] .* [Y | Row2] | dotprodrows_([X | Row1], Rows2)])
.
dotprodrows_([_G207698 | _G207699], []) := [].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	Gaussian elimination
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gauss(_G207698) 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] | _G207702] ^ 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].