Download or view Matrix.frink in plain text format
/** This is an incomplete class that implements matrix operations. Please
feel free to implement missing functions! */
class Matrix
{
/** The array that will actually store the matrix data. Note that this
is a raw 2-dimensional Frink array, and the indices are zero-based in
contrast with most matrix notation which is one-based. The class should
try to enforce that this is a rectangular 2-dimensional array, but the
code doesn't currently enforce that with some constructors. */
var array
/** The number of rows in the matrix. */
var rows
/** The number of columns in the matrix. */
var cols
/** Construct a new Matrix with the specified number of rows and columns
and the specified default value for each element. */
new[rowCount, colCount, default=0] :=
{
rows = rowCount
cols = colCount
array = makeArray[[rows, cols], default]
}
/** Construct a new Matrix from a rectangular 2-dimensional array. This
currently does not verify that the array is rectangular. */
new[a is array] :=
{
rows = length[a]
cols = length[a@0]
array = a
}
/** Sets the value of the specified element, using 1-based coordinates, as
is common in matrix notation. Stupid mathematicians. */
set[row, col, val] := array@(row-1)@(col-1) = val
/** Gets the value of the specified element, using 1-based coordinates, as
is common in matrix notation. */
get[row, col] := array@(row-1)@(col-1)
/** Multiply two matrices. The number of columns in column a should equal
the number of rows in column b. The resulting matrix will have the
same number of rows as matrix a and the same number of columns as
matrix b.
TODO: Use a faster algorithm for large arrays? This is O(n^3).
TODO: Optimize this for small arrays by precalculating.
*/
class multiply[a, b] :=
{
if a.cols != b.rows
return undef
else
items = a.cols
resRows = a.rows
resCols = b.cols
resultArray = makeArray[[resRows, resCols], 0]
aa = a.array
bb = b.array
for row = 0 to resRows-1
for col = 0 to resCols-1
{
sum = 0
for k = 0 to items-1
sum = sum + aa@row@k * bb@k@col
resultArray@row@col = sum
}
return new Matrix[resultArray]
}
/** Multiply this matrix by the specified matrix and return the result. */
multiply[b] := multiply[this, b]
/** Multiplies all elements of a matrix by a scalar. */
multiplyByScalar[s] :=
{
a = makeArray[[rows, cols]]
for rowNum = 0 to rows-1
{
row = array@rowNum
for colNum = 0 to cols-1
a@rowNum@colNum = row@colNum * s
}
return new Matrix[a]
}
/** Transposes the Matrix and returns the new transposed Matrix. This
uses the built-in array.transpose[] method.
*/
transpose[] := new Matrix[array.transpose[]]
/** Returns true if this Matrix is square, false otherwise. */
isSquare[] := rows == cols
/** Exchange the specified-numbered rows in the matrix. The rows are
1-indexed to match standard matrix notation. */
exchangeRows[r1, r2] :=
{
temp = array@(r2-1)
array@(r2-1) = array@(r1-1)
array@(r1-1) = temp
}
/** Exchange the specified-numbered columns in the matrix. The columns are
1-indexed to match standard matrix notation. */
exchangeCols[c1, c2] :=
{
for r = 0 to rows-1
{
temp = array@r@(c2-1)
array@r@(c2-1) = array@r@(c1-1)
array@r@(c1-1) = temp
}
}
/** Returns a matrix as a string with rows separated with newlines. */
toString[] :=
{
result = ""
for r = 0 to rows-1
result = result + " " + array@r + "\n"
return result
}
/** Returns the matrix formatted by formatTable with rows separated with
newlines. */
toTable[] :=
{
return formatTable[array]
}
/** Formats a matrix with Unicode box drawing characters. The result is a
single string. */
formatMatrix[] :=
{
return joinln[formatMatrixLines[]]
}
/** Returns the matrix formatted as a matrix with Unicode box drawing
characters. The result is an array of lines that can be further
formatted.
*/
formatMatrixLines[] :=
{
m = formatTableLines[array]
rows = length[m]
width = length[m@0]
result = new array
result.push["\u250c" + repeat[" ", width] + "\u2510"]
for n=0 to rows-1
result.push["\u2502" + m@n + "\u2502"]
result.push["\u2514" + repeat[" ", width] + "\u2518"]
return result
}
/** Creates the LU (lower-upper) decomposition of the matrix.
This creates two matrices, L and U, where L has the value 1 on the
diagonal from top left to bottom right, like:
1 0 0
L = L21 1 0
L31 L32 1
U11 U12 U13
U = 0 U22 U23
0 0 U33
This is basically like solving equations by Gaussian elimination. */
LUDecompose[] :=
{
}
/** This uses Crout's method to decompose a square matrix into an lower and
upper triangular matrix.
See:
https://en.wikipedia.org/wiki/Crout_matrix_decomposition
This creates two matrices, L and U, where U has the value 1 on the
diagonal from top left to bottom right, like:
L11 0 0
L = L21 L22 0
L31 L32 L33
1 U12 U13
U = 0 1 U23
0 0 1
returns [L, U]
The original matrix can be obtained as L.multiply[U]
Note that the Crout algorithm fails on certain matrices, for example
m = new Matrix[[[1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[-1,1,1,1]]]
See:
https://semath.info/src/inverse-cofactor-ex4.html
for more on that matrix.
https://en.wikipedia.org/wiki/Crout_matrix_decomposition
*/
LUDecomposeCrout[] :=
{
n = rows
L = new array[[rows, cols], 0]
U = new array[[rows, cols], 0]
sum = 0
for i=0 to n-1
U@i@i = 1
for j=0 to n-1
{
for i=j to n-1
{
sum = 0
for k=0 to j-1
sum = sum + L@i@k * U@k@j
L@i@j = array@i@j - sum
}
for i=j to n-1
{
sum = 0
for k=0 to j-1
sum = sum + L@j@k * U@k@i
if L@j@j == 0
{
println["Matrix.LUDecomposeCrout: det(L) is zero. Can't divide by zero."]
return undef
}
U@j@i = (array@j@i - sum) / L@j@j
}
}
return [new Matrix[L], new Matrix[U]]
}
/** Another attempt at LU decomposition. See:
http://www.mosismath.com/Matrix_LU/Martix_LU.html
*/
LUDecomposeCrout2[] :=
{
if ! isSquare[]
{
println["Matrix.LUDecomposeCrout2 only works on square arrays."]
return undef
}
n = rows
L = new array[[n,n], 0]
U = new array[[n,n], 0]
a = 0
for i = 0 to n-1
L@i@i = 1
for j = 0 to n-1
{
for i = 0 to n-1
{
if i >= j
{
U@j@i = array@j@i
for k = 0 to j-1
{
U@j@i = U@j@i - U@k@i * L@j@k
}
}
if i > j
{
L@i@j = array@i@j
for k = 0 to j-1
L@i@j = L@i@j - U@k@j * L@i@k
diag = U@j@j
if diag == 0
{
println["Matrix.LUDecomposeCrout2: det(L) is zero. Can't divide by zero."]
println[formatMatrix[L]]
println[formatMatrix[U]]
return undef
}
L@i@j = L@i@j / diag
}
}
}
return [new Matrix[L], new Matrix[U]]
}
/** This uses Doolittle's method to decompose a square matrix into an lower
and upper triangular matrix.
This creates two matrices, L and U, where L has the value 1 on the
diagonal from top left to bottom right, like:
1 0 0
L = L21 1 0
L31 L32 1
U11 U12 U13
U = 0 U22 U23
0 0 U33
The original matrix can be obtained as L.multiply[U]
Note that the Doolittle algorithm fails on certain matrices, for example
m = new Matrix[[[1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[-1,1,1,1]]]
See:
https://semath.info/src/inverse-cofactor-ex4.html
for more on that matrix.
See:
https://www.geeksforgeeks.org/doolittle-algorithm-lu-decomposition/
*/
LUDecomposeDoolittle[] :=
{
if ! isSquare[]
{
println["Matrix.LUDecomposeDoolittle only works on square arrays."]
return undef
}
n = rows
L = new array[[n,n], 0]
U = new array[[n,n], 0]
// Decomposing matrix into Upper and Lower
// triangular matrix
for i = 0 to n-1
{
// Upper Triangular
for k = i to n-1
{
// Summation of L(i, j) * U(j, k)
sum = 0
for j = 0 to i-1
sum = sum + (L@i@j * U@j@k)
// Evaluating U(i, k)
U@i@k = array@i@k - sum
}
// Lower Triangular
for k = i to n-1
{
if i == k
L@i@i = 1 // Diagonal as 1
else
{
// Summation of L(k, j) * U(j, i)
sum = 0
for j = 0 to i-1
sum = sum + L@k@j * U@j@i
diag = U@i@i
if (diag == 0)
{
println["Matrix.LUDecomposeDoolittle: U@$i@$i is zero when solving for L@$k@$i. Can't divide by zero."]
println[formatMatrix[L]]
println[formatMatrix[U]]
return undef
}
// Evaluating L(k, i)
L@k@i = (array@k@i - sum) / diag
}
}
}
return [new Matrix[L], new Matrix[U]]
}
/** Returns the determinant of a square matrix. */
det[] :=
{
if ! isSquare[]
{
println["Determinants are only defined for square matrices! Size was $rows x $cols"]
return undef
}
// It's much faster to use precalculated determinant formulas for small
// matrices than, to, say LUDecompose the things.
if rows == 1
return array@0@0
if rows == 2
return array@0@0 array@1@1 - array@0@1 array@1@0 // ad - bc
if rows == 3
return -array@0@2 array@1@1 array@2@0 +
array@0@1 array@1@2 array@2@0 +
array@0@2 array@1@0 array@2@1 -
array@0@0 array@1@2 array@2@1 -
array@0@1 array@1@0 array@2@2 +
array@0@0 array@1@1 array@2@2
if rows == 4
return array@0@1 array@1@3 array@2@2 array@3@0 -
array@0@1 array@1@2 array@2@3 array@3@0 -
array@0@0 array@1@3 array@2@2 array@3@1 +
array@0@0 array@1@2 array@2@3 array@3@1 -
array@0@1 array@1@3 array@2@0 array@3@2 +
array@0@0 array@1@3 array@2@1 array@3@2 +
array@0@1 array@1@0 array@2@3 array@3@2 -
array@0@0 array@1@1 array@2@3 array@3@2 +
array@0@3 (array@1@2 array@2@1 array@3@0 -
array@1@1 array@2@2 array@3@0 -
array@1@2 array@2@0 array@3@1 +
array@1@0 array@2@2 array@3@1 +
array@1@1 array@2@0 array@3@2 -
array@1@0 array@2@1 array@3@2) +
(array@0@1 array@1@2 array@2@0 -
array@0@0 array@1@2 array@2@1 -
array@0@1 array@1@0 array@2@2 +
array@0@0 array@1@1 array@2@2) array@3@3 +
array@0@2 (-array@1@3 array@2@1 array@3@0 +
array@1@1 array@2@3 array@3@0 +
array@1@3 array@2@0 array@3@1 -
array@1@0 array@2@3 array@3@1 -
array@1@1 array@2@0 array@3@3 +
array@1@0 array@2@1 array@3@3)
if rows == 5
return array@0@2 array@1@4 array@2@3 array@3@1 array@4@0 -
array@0@2 array@1@3 array@2@4 array@3@1 array@4@0 -
array@0@1 array@1@4 array@2@3 array@3@2 array@4@0 +
array@0@1 array@1@3 array@2@4 array@3@2 array@4@0 -
array@0@2 array@1@4 array@2@1 array@3@3 array@4@0 +
array@0@1 array@1@4 array@2@2 array@3@3 array@4@0 +
array@0@2 array@1@1 array@2@4 array@3@3 array@4@0 -
array@0@1 array@1@2 array@2@4 array@3@3 array@4@0 +
array@0@2 array@1@3 array@2@1 array@3@4 array@4@0 -
array@0@1 array@1@3 array@2@2 array@3@4 array@4@0 -
array@0@2 array@1@1 array@2@3 array@3@4 array@4@0 +
array@0@1 array@1@2 array@2@3 array@3@4 array@4@0 -
array@0@2 array@1@4 array@2@3 array@3@0 array@4@1 +
array@0@2 array@1@3 array@2@4 array@3@0 array@4@1 +
array@0@0 array@1@4 array@2@3 array@3@2 array@4@1 -
array@0@0 array@1@3 array@2@4 array@3@2 array@4@1 +
array@0@2 array@1@4 array@2@0 array@3@3 array@4@1 -
array@0@0 array@1@4 array@2@2 array@3@3 array@4@1 -
array@0@2 array@1@0 array@2@4 array@3@3 array@4@1 +
array@0@0 array@1@2 array@2@4 array@3@3 array@4@1 -
array@0@2 array@1@3 array@2@0 array@3@4 array@4@1 +
array@0@0 array@1@3 array@2@2 array@3@4 array@4@1 +
array@0@2 array@1@0 array@2@3 array@3@4 array@4@1 -
array@0@0 array@1@2 array@2@3 array@3@4 array@4@1 +
array@0@1 array@1@4 array@2@3 array@3@0 array@4@2 -
array@0@1 array@1@3 array@2@4 array@3@0 array@4@2 -
array@0@0 array@1@4 array@2@3 array@3@1 array@4@2 +
array@0@0 array@1@3 array@2@4 array@3@1 array@4@2 -
array@0@1 array@1@4 array@2@0 array@3@3 array@4@2 +
array@0@0 array@1@4 array@2@1 array@3@3 array@4@2 +
array@0@1 array@1@0 array@2@4 array@3@3 array@4@2 -
array@0@0 array@1@1 array@2@4 array@3@3 array@4@2 +
array@0@1 array@1@3 array@2@0 array@3@4 array@4@2 -
array@0@0 array@1@3 array@2@1 array@3@4 array@4@2 -
array@0@1 array@1@0 array@2@3 array@3@4 array@4@2 +
array@0@0 array@1@1 array@2@3 array@3@4 array@4@2 +
array@0@2 array@1@4 array@2@1 array@3@0 array@4@3 -
array@0@1 array@1@4 array@2@2 array@3@0 array@4@3 -
array@0@2 array@1@1 array@2@4 array@3@0 array@4@3 +
array@0@1 array@1@2 array@2@4 array@3@0 array@4@3 -
array@0@2 array@1@4 array@2@0 array@3@1 array@4@3 +
array@0@0 array@1@4 array@2@2 array@3@1 array@4@3 +
array@0@2 array@1@0 array@2@4 array@3@1 array@4@3 -
array@0@0 array@1@2 array@2@4 array@3@1 array@4@3 +
array@0@1 array@1@4 array@2@0 array@3@2 array@4@3 -
array@0@0 array@1@4 array@2@1 array@3@2 array@4@3 -
array@0@1 array@1@0 array@2@4 array@3@2 array@4@3 +
array@0@0 array@1@1 array@2@4 array@3@2 array@4@3 +
array@0@2 array@1@1 array@2@0 array@3@4 array@4@3 -
array@0@1 array@1@2 array@2@0 array@3@4 array@4@3 -
array@0@2 array@1@0 array@2@1 array@3@4 array@4@3 +
array@0@0 array@1@2 array@2@1 array@3@4 array@4@3 +
array@0@1 array@1@0 array@2@2 array@3@4 array@4@3 -
array@0@0 array@1@1 array@2@2 array@3@4 array@4@3 +
array@0@4 (array@1@1 array@2@3 array@3@2 array@4@0 -
array@1@1 array@2@2 array@3@3 array@4@0 -
array@1@0 array@2@3 array@3@2 array@4@1 +
array@1@0 array@2@2 array@3@3 array@4@1 -
array@1@1 array@2@3 array@3@0 array@4@2 +
array@1@0 array@2@3 array@3@1 array@4@2 +
array@1@1 array@2@0 array@3@3 array@4@2 -
array@1@0 array@2@1 array@3@3 array@4@2 +
array@1@3 (array@2@2 array@3@1 array@4@0 -
array@2@1 array@3@2 array@4@0 -
array@2@2 array@3@0 array@4@1 +
array@2@0 array@3@2 array@4@1 +
array@2@1 array@3@0 array@4@2 -
array@2@0 array@3@1 array@4@2) +
(array@1@1 array@2@2 array@3@0 -
array@1@0 array@2@2 array@3@1 -
array@1@1 array@2@0 array@3@2 +
array@1@0 array@2@1 array@3@2) array@4@3 +
array@1@2 (-array@2@3 array@3@1 array@4@0 +
array@2@1 array@3@3 array@4@0 +
array@2@3 array@3@0 array@4@1 -
array@2@0 array@3@3 array@4@1 -
array@2@1 array@3@0 array@4@3 +
array@2@0 array@3@1 array@4@3)) +
(array@0@2 (-array@1@3 array@2@1 array@3@0 +
array@1@1 array@2@3 array@3@0 +
array@1@3 array@2@0 array@3@1 -
array@1@0 array@2@3 array@3@1 -
array@1@1 array@2@0 array@3@3 +
array@1@0 array@2@1 array@3@3) +
array@0@1 (array@1@3 array@2@2 array@3@0 -
array@1@2 array@2@3 array@3@0 -
array@1@3 array@2@0 array@3@2 +
array@1@0 array@2@3 array@3@2 +
array@1@2 array@2@0 array@3@3 -
array@1@0 array@2@2 array@3@3) +
array@0@0 (-array@1@3 array@2@2 array@3@1 +
array@1@2 array@2@3 array@3@1 +
array@1@3 array@2@1 array@3@2 -
array@1@1 array@2@3 array@3@2 -
array@1@2 array@2@1 array@3@3 +
array@1@1 array@2@2 array@3@3)) array@4@4 +
array@0@3 (-array@1@1 array@2@4 array@3@2 array@4@0 +
array@1@1 array@2@2 array@3@4 array@4@0 +
array@1@0 array@2@4 array@3@2 array@4@1 -
array@1@0 array@2@2 array@3@4 array@4@1 +
array@1@1 array@2@4 array@3@0 array@4@2 -
array@1@0 array@2@4 array@3@1 array@4@2 -
array@1@1 array@2@0 array@3@4 array@4@2 +
array@1@0 array@2@1 array@3@4 array@4@2 +
array@1@4 (-array@2@2 array@3@1 array@4@0 +
array@2@1 array@3@2 array@4@0 +
array@2@2 array@3@0 array@4@1 -
array@2@0 array@3@2 array@4@1 -
array@2@1 array@3@0 array@4@2 +
array@2@0 array@3@1 array@4@2) -
array@1@1 array@2@2 array@3@0 array@4@4 +
array@1@0 array@2@2 array@3@1 array@4@4 +
array@1@1 array@2@0 array@3@2 array@4@4 -
array@1@0 array@2@1 array@3@2 array@4@4 +
array@1@2 (array@2@4 array@3@1 array@4@0 -
array@2@1 array@3@4 array@4@0 -
array@2@4 array@3@0 array@4@1 +
array@2@0 array@3@4 array@4@1 +
array@2@1 array@3@0 array@4@4 -
array@2@0 array@3@1 array@4@4))
// If we fell through to here, we have a larger matrix and have to try
// to find its determinant through more inefficent methods.
// TODO: Calculate determinant in terms of other determinants instead of
// using LUDecompose when LUDecompose doesn't work.
product = 1
[L, U] = LUDecomposeDoolittle[]
// This will happen if the matrix is singular.
if U == undef
return undef
// Multiply diagonals of the lower triangular matrix. The
// upper triangular matrix has all ones on the diagonal.
// Should there be a permutation matrix here to get the signs
// right?
for i=0 to rows-1
product = product * (U.array)@i@i
return product
}
/** Create a square identity matrix with the specified number of rows and
columns. The elements on the diagonal will be set to 1, the rest to
zero. This requires Frink 2020-04-19 or later. */
class makeIdentity[dimension] :=
{
return new Matrix[makeArray[[dimension, dimension], {|a,b| a==b ? 1 : 0}]]
}
/** Makes a diagonal matrix. This is passed in an array of elements (e.g.
[1,2,3] that will make up the diagonal elements. This requires Frink
2020-04-22 or later.
*/
class makeDiagonal[array] :=
{
d = length[array]
return new Matrix[makeArray[[d,d], {|a,b,data| a==b ? data@a : 0}, array]]
}
/** Returns a matrix with the specified row and column removed. Row and
column numbers are 1-indexed. This is used in many matrix operations
including calculation of determinant or of the adjugate/adjoint matrix.
*/
removeRowColumn[rowToRemove, colToRemove] :=
{
a = makeArray[[rows-1, cols-1]]
newRow = 0
ROW:
for origRow = 0 to rows-1
{
if origRow == rowToRemove-1
next ROW
newCol = 0
COL:
for origCol = 0 to rows-1
{
if origCol == colToRemove-1
next COL
a@newRow@newCol = array@origRow@origCol
newCol = newCol + 1
}
newRow = newRow + 1
}
return new Matrix[a]
}
/** Returns the adjugate or adjoint matrix. See:
https://semath.info/src/inverse-cofactor-ex4.html
TODO: Calculate hardcoded adjugate matrices for small matrices because
this is expensive
*/
adjugate[] :=
{
a = makeArray[[rows,cols]]
for i = 0 to rows-1
for j = 0 to cols-1
a@i@j = (-1)^((i+1)+(j+1)) * removeRowColumn[j+1, i+1].det[]
return new Matrix[a]
}
/** Returns the inverse of a matrix. See:
https://semath.info/src/inverse-cofactor-ex4.html
TODO: Calculate hardcoded inverse matrices for small matrices because
this is expensive
*/
inverse[] :=
{
if ! isSquare[]
{
println["Matrix.inverse only works on square arrays."]
return undef
}
// It's much faster to calculate inverse matrices with hard-coded
// equations for small matrices.
if rows == 1
return 1 / array@0@0
if rows == 2
{
invdet = 1/det[]
// [ d -b ]
// [ -c a ] / det
return (new Matrix[[[array@1@1, -array@0@1],
[-array@1@0, array@0@0]]]).multiplyByScalar[invdet]
}
if rows == 3
{
invdet = 1/det[]
return (new Matrix[[[-array@1@2 array@2@1 + array@1@1 array@2@2,
array@0@2 array@2@1 - array@0@1 array@2@2,
-array@0@2 array@1@1 + array@0@1 array@1@2],
[ array@1@2 array@2@0 - array@1@0 array@2@2,
-array@0@2 array@2@0 + array@0@0 array@2@2,
array@0@2 array@1@0 - array@0@0 array@1@2],
[-array@1@1 array@2@0 + array@1@0 array@2@1,
array@0@1 array@2@0 - array@0@0 array@2@1,
-array@0@1 array@1@0 + array@0@0 array@1@1]]]).multiplyByScalar[invdet]
}
if rows == 4
{
invdet = 1/det[]
return (new Matrix[[[-array@1@3 array@2@2 array@3@1 +
array@1@2 array@2@3 array@3@1 +
array@1@3 array@2@1 array@3@2 -
array@1@1 array@2@3 array@3@2 -
array@1@2 array@2@1 array@3@3 +
array@1@1 array@2@2 array@3@3,
array@0@3 array@2@2 array@3@1 -
array@0@2 array@2@3 array@3@1 -
array@0@3 array@2@1 array@3@2 +
array@0@1 array@2@3 array@3@2 +
array@0@2 array@2@1 array@3@3 -
array@0@1 array@2@2 array@3@3,
-array@0@3 array@1@2 array@3@1 +
array@0@2 array@1@3 array@3@1 +
array@0@3 array@1@1 array@3@2 -
array@0@1 array@1@3 array@3@2 -
array@0@2 array@1@1 array@3@3 +
array@0@1 array@1@2 array@3@3,
array@0@3 array@1@2 array@2@1 -
array@0@2 array@1@3 array@2@1 -
array@0@3 array@1@1 array@2@2 +
array@0@1 array@1@3 array@2@2 +
array@0@2 array@1@1 array@2@3 -
array@0@1 array@1@2 array@2@3],
[array@1@3 array@2@2 array@3@0 -
array@1@2 array@2@3 array@3@0 -
array@1@3 array@2@0 array@3@2 +
array@1@0 array@2@3 array@3@2 +
array@1@2 array@2@0 array@3@3 -
array@1@0 array@2@2 array@3@3,
-array@0@3 array@2@2 array@3@0 +
array@0@2 array@2@3 array@3@0 +
array@0@3 array@2@0 array@3@2 -
array@0@0 array@2@3 array@3@2 -
array@0@2 array@2@0 array@3@3 +
array@0@0 array@2@2 array@3@3,
array@0@3 array@1@2 array@3@0 -
array@0@2 array@1@3 array@3@0 -
array@0@3 array@1@0 array@3@2 +
array@0@0 array@1@3 array@3@2 +
array@0@2 array@1@0 array@3@3 -
array@0@0 array@1@2 array@3@3,
-array@0@3 array@1@2 array@2@0 +
array@0@2 array@1@3 array@2@0 +
array@0@3 array@1@0 array@2@2 -
array@0@0 array@1@3 array@2@2 -
array@0@2 array@1@0 array@2@3 +
array@0@0 array@1@2 array@2@3],
[-array@1@3 array@2@1 array@3@0 +
array@1@1 array@2@3 array@3@0 +
array@1@3 array@2@0 array@3@1 -
array@1@0 array@2@3 array@3@1 -
array@1@1 array@2@0 array@3@3 +
array@1@0 array@2@1 array@3@3,
array@0@3 array@2@1 array@3@0 -
array@0@1 array@2@3 array@3@0 -
array@0@3 array@2@0 array@3@1 +
array@0@0 array@2@3 array@3@1 +
array@0@1 array@2@0 array@3@3 -
array@0@0 array@2@1 array@3@3,
-array@0@3 array@1@1 array@3@0 +
array@0@1 array@1@3 array@3@0 +
array@0@3 array@1@0 array@3@1 -
array@0@0 array@1@3 array@3@1 -
array@0@1 array@1@0 array@3@3 +
array@0@0 array@1@1 array@3@3,
array@0@3 array@1@1 array@2@0 -
array@0@1 array@1@3 array@2@0 -
array@0@3 array@1@0 array@2@1 +
array@0@0 array@1@3 array@2@1 +
array@0@1 array@1@0 array@2@3 -
array@0@0 array@1@1 array@2@3],
[array@1@2 array@2@1 array@3@0 -
array@1@1 array@2@2 array@3@0 -
array@1@2 array@2@0 array@3@1 +
array@1@0 array@2@2 array@3@1 +
array@1@1 array@2@0 array@3@2 -
array@1@0 array@2@1 array@3@2,
-array@0@2 array@2@1 array@3@0 +
array@0@1 array@2@2 array@3@0 +
array@0@2 array@2@0 array@3@1 -
array@0@0 array@2@2 array@3@1 -
array@0@1 array@2@0 array@3@2 +
array@0@0 array@2@1 array@3@2,
array@0@2 array@1@1 array@3@0 -
array@0@1 array@1@2 array@3@0 -
array@0@2 array@1@0 array@3@1 +
array@0@0 array@1@2 array@3@1 +
array@0@1 array@1@0 array@3@2 -
array@0@0 array@1@1 array@3@2,
-array@0@2 array@1@1 array@2@0 +
array@0@1 array@1@2 array@2@0 +
array@0@2 array@1@0 array@2@1 -
array@0@0 array@1@2 array@2@1 -
array@0@1 array@1@0 array@2@2 +
array@0@0 array@1@1 array@2@2]]]).multiplyByScalar[invdet]
}
// Otherwise use the adjugate function.
return adjugate[].multiplyByScalar[1/det[]]
}
/** Returns the Kronecker product of a and b.
https://en.wikipedia.org/wiki/Kronecker_product
*/
class KroneckerProduct[a is Matrix, b is Matrix] :=
{
m = a.rows
n = a.cols
p = b.rows
q = b.cols
newRows = m p
newCols = n q
n = new array[[newRows, newCols], 0]
for i=0 to newRows-1
for j=0 to newCols-1
n@i@j = (a.array)@(i div p)@(j div q) * (b.array)@(i mod p)@(j mod q)
return new Matrix[n]
}
/** Return the Kronecker product of this matrix and b. */
KroneckerProduct[b is Matrix] := KroneckerProduct[this, b]
}
"class Matrix loaded OK."
Download or view Matrix.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19302 days, 19 hours, 24 minutes ago.