/** 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
   }

   /** Gets the specified (1-based) column as a row array. */
   getColumnAsArray[col] :=
   {
      return deepCopy[array.getColumn[col-1]]
   }

   /** Gets the specified (1-based) row as a row array. */
   getRowAsArray[row] :=
   {
      return deepCopy[array@(row-1)]
   }
   
   /** 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[] :=
   {
      return LUDecomposeCrout[]
   }

   
   /** 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 algorithm for Crout LU decomposition.  See:

       http://www.mosismath.com/Matrix_LU/Martix_LU.html

       This returns [L, U] as matrices which may be transposes of the matrices
       returned by LUDecomposeCrout[] .  See its notes for more about this
       decomposition and its properties.
   */
   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 the Cholesky-Banachiewicz algorithm to decompose a square
       Hermitian matrix into a lower triangular matrix.  If you transpose the
       lower triangular matrix L by L.transpose[], you get an upper triangular
       matrix which is symmetrical to the lower triangular matrix.
       Multiplying the lower triangular matrix by the upper triangular matrix,
       that is,

       L.multiply[L.transpose[]]
   
       gives you back the original matrix!

       See:  https://en.m.wikipedia.org/wiki/Cholesky_decomposition
   */
   CholeskyB[] :=
   {
      if ! isHermitian[]
      {
         println["Matrix.CholeskyB only works on square Hermitian matrices."]
         return undef
      }

      n = rows
      L = new array[[n,n], 0]

      for i = 0 to n-1
         for j = 0 to i
         {
            sum = 0
            for k = 0 to j-1
               sum = sum + L@i@k * L@j@k

            if i == j
               L@i@j = sqrt[array@i@i - sum]
            else
               L@i@j = (1 / L@j@j * (array@i@j -sum))
         }

      return new Matrix[L]
   }
   
   /** This uses the Cholesky-Crout algorithm to decompose a square Hermitian
       matrix into a lower triangular matrix.  If you transpose the lower
       triangular matrix L by L.transpose[], you get an upper triangular
       matrix which is symmetrical to the lower triangular matrix.
       Multiplying the lower triangular matrix by the upper triangular matrix,
       that is,

       L.multiply[L.transpose[]]
   
       gives you back the original matrix!

       See:  https://en.m.wikipedia.org/wiki/Cholesky_decomposition
   */
   CholeskyCrout[] :=
   {
      if ! isHermitian[]
      {
         println["Matrix.CholeskyCrout only works on square Hermitian matrices."]
         return undef
      }

      n = rows
      L = new array[[n,n], 0]

      for j = 0 to n-1
      {
         sum = 0
         for k = 0 to j-1
            sum = sum + (L@j@k)^2

         L@j@j = sqrt[array@j@j - sum]

         for i = j+1 to n-1
         {
            sum = 0
            for k = 0 to j-1
               sum = sum + L@i@k * L@j@k

            L@i@j = (1 / L@j@j * (array@i@j -sum))
         }
      }

      return new Matrix[L]
   }
   
   /** 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]

   
   /** Returns true if both matrices are equal (that is, they have the same
       dimensions and all array elements are equal.)
   */
   equals[other is Matrix] :=
   {
      return rows==other.rows and cols==other.cols and array==other.array
   }

   
   /** Returns a new Matrix where each element is the complex conjugate of
       the corresponding element in original Matrix. */
   conjugate[] :=
   {
      c = makeArray[[rows, cols], {|r,c,data| conjugate[data@r@c]}, array]
      return new Matrix[c]
   }


   /** Returns a new Matrix where the Matrix is first transposed and then
       each element is the complex conjugate of
       the corresponding element in the transposed Matrix. */
   conjugateTranspose[] :=
   {
      t = array.transpose[]
      c = makeArray[[rows, cols], {|r,c,data| conjugate[data@r@c]}, t]
      return new Matrix[c]
   }

   /** Returns true if this matrix is "Hermitian", or "self-adjoint", that is,
       the matrix must be square and equal to its own conjugate transpose.  In
       other words, for all elements, array@i@j == conjugate[array@j@i] where
       conjugate is the complex conjugate of a number.  */
   isHermitian[] :=
   {
      if rows != cols
         return false

      for i=0 to rows-1
         for j = i+1 to rows-1
            if array@i@j != conjugate[array@j@i]
               return false

      return true           
   }

   /** Performs the QR decomposition of a matrix.

       This is based on the (real-only) implementation at:
       https://github.com/fiji/Jama/blob/master/src/main/java/Jama/QRDecomposition.java

       This is an internal implementation, primarily for use by the leastSquares
       method, that returns raw arrays [QR, rdiag].  If you're solving for
       least squares, use that routine directly instead of this.
   */
   QRDecomposeInternal[] :=
   {
      QR = deepCopy[array]   // Copy the original matrix's array
      rdiag = new array[[cols], 0]       // Stores the diagonal

      for k = 0 to cols-1
      {
         // Compute 2-norm of k-th column
         nrm = 0 QR@0@k   // Make units work right
         for i = k to rows-1
            nrm = sqrt[nrm^2 + (QR@i@k)^2]

         if nrm != 0 QR@0@k  // Make units work right
         {
            // Form k-th Householder vector
            if isNegativeUnit[QR@k@k]
               nrm = -nrm

            for i = k to rows-1
               QR@i@k = QR@i@k / nrm

            QR@k@k = QR@k@k + 1

            // Apply transformation to remaining columns
            for j = k+1 to cols-1
            {
               s = undef
               for i=k to rows-1
               {
                  if s == undef
                     s = QR@i@k * QR@i@j
                  else
                     s = s + QR@i@k * QR@i@j
               }

               s = -s / QR@k@k

               for i = k to rows-1
                  QR@i@j = QR@i@j + s * QR@i@k
            }
         }

         rdiag@k = -nrm
      }

      return [QR, rdiag]
   }

   /** This performs a QR decomposition of this matrix and returns [Q, R]
       as new matrices.

       Q is the orthogonal factor.
       R is the upper triangular factor.

       If you are doing a least-squares solve, call leastSquares directly
       instead.

       THINK ABOUT:  Should this return the Householder vectors and rdiag?
   */
   QRDecompose[] :=
   {
      // The Q and R matrices are packed into the results of the following.
      [QR, rdiag] = QRDecomposeInternal[]

      // Extract the Q matrix
      Q = new array[[rows, cols], 0]
      for k = cols-1 to 0 step -1
      {
         for i = 0 to rows-1
            Q@i@k = 0

         Q@k@k = 1

         for j = k to cols-1
         {
            if QR@k@k != 0 QR@k@k  // Make units work right
            {
               s = 0 (QR@0@k) * Q@0@j // Make units work right
               for i = k to rows-1
                  s = s + QR@i@k * Q@i@j

               s = -s / QR@k@k
               
               for i = k to rows-1
                  Q@i@j = Q@i@j + s * QR@i@k
            }
         }
      }

      // Extract the R matrix
      R = new array[[cols, cols], 0]
      for i = 0 to cols-1
      {
         for j = 0 to cols-1
         {
            if i < j
               R@i@j = QR@i@j
            else
               if i == j
                  R@i@j = rdiag@i
               else
                  R@i@j = 0
         }
      }

      return [new Matrix[Q], new Matrix[R]]
   }

   /** Return the least-squares solution (called X) of the system
       A * X = B
       where this Matrix is A and parameter B is B.

       This can be performed on an overdetermined system, where there are
       more measurements than equations.

       TODO:  Try to solve exactly using x = A.inverse[] * b when the system
       is not overdetermined?

       See MatrixQRTest.frink for examples.

       This is the best discussion I've seen of least-squares fitting:
        https://www.aleksandrhovhannisyan.com/blog/the-method-of-least-squares/
        https://www.aleksandrhovhannisyan.com/blog/least-squares-fitting/
   */
   leastSquares[B is Matrix] :=
   {
      if rows != B.rows
      {
         println["Matrix.leastSquares:  Matrix row dimensions must agree."]
         return undef
      }
      
      [QR, rdiag] = QRDecomposeInternal[]

      // Check if the matrix is "full rank," that is, that each row is
      // independent (not a multiple of) other rows.
      for j = 0 to cols-1
      {
         if rdiag@j == 0 rdiag@j
         {
            println["Matrix.leastSquares: Row entries are not all independent."]
            return undef
         }
      }

      X = deepCopy[B.array]   // Copy B's array

      // Compute Y = Q.transpose * B
      for k = 0 to cols-1
      {
         for j = 0 to B.cols-1
         {
            s = undef
            for i = k to rows-1
               if s == undef
                  s = QR@i@k * X@i@j
               else
                  s = s + QR@i@k * X@i@j

            s = -s / QR@k@k

            for i = k to rows-1
               X@i@j = X@i@j + s * QR@i@k
         }
      }

      // Solve R * X = y
      for k = cols-1 to 0 step -1
      {
         for j = 0 to B.cols-1
            X@k@j = X@k@j / rdiag@k

         for i = 0 to k-1
            for j = 0 to B.cols-1
               X@i@j = X@i@j - X@k@j * QR@i@k
      }

      X1 = new Matrix[X]
      return X1.getSubMatrix[0, cols-1, 0, B.cols-1]
   }

   /** Gets a submatrix of this matrix. */
   getSubMatrix[fromRow, toRow, fromCol, toCol] :=
   {
      a = new array[[toRow-fromRow+1, toCol-fromCol+1], 0]
      for row = fromRow to toRow
         for col = fromCol to toCol
            a@(row-fromRow)@(col-fromCol) = array@row@col

//      println["a is $a"]
      return new Matrix[a]       
   }

   /** Rounds values to the nearest integer if they are less than the specified
       relative error away from an integer and returns a new array. */
   roundToInt[relerror = 1e-14] :=
   {
      ret = new array[[rows, cols]]
      for r = 0 to rows-1
         for c = 0 to cols-1
         {
            v = array@r@c
            rnd = round[v]
            if abs[(rnd-v)/rnd] <= relerror
               ret@r@c = rnd
            else
               ret@r@c = v
         }

      return new Matrix[ret]
   }
}



"class Matrix loaded OK."