| | | 1 | | using System.Diagnostics; |
| | | 2 | | using System.Text; |
| | | 3 | | |
| | | 4 | | namespace Towel.Mathematics; |
| | | 5 | | |
| | | 6 | | /// <summary>A matrix of arbitrary dimensions implemented as a flattened array.</summary> |
| | | 7 | | /// <typeparam name="T">The numeric type of this Matrix.</typeparam> |
| | | 8 | | [DebuggerDisplay("{" + nameof(DebuggerString) + "}")] |
| | | 9 | | public class Matrix<T> |
| | | 10 | | { |
| | | 11 | | internal readonly T[] _matrix; |
| | | 12 | | internal int _rows; |
| | | 13 | | internal int _columns; |
| | | 14 | | |
| | | 15 | | #region Properties |
| | | 16 | | |
| | 11 | 17 | | internal int Length => _matrix.Length; |
| | | 18 | | /// <summary>The number of rows in the matrix.</summary> |
| | 629 | 19 | | public int Rows => _rows; |
| | | 20 | | /// <summary>The number of columns in the matrix.</summary> |
| | 3642 | 21 | | public int Columns => _columns; |
| | | 22 | | /// <summary>Determines if the matrix is square.</summary> |
| | 38 | 23 | | public bool IsSquare => _rows == _columns; |
| | | 24 | | /// <summary>Determines if the matrix is a vector.</summary> |
| | 0 | 25 | | public bool IsVector { get { return _columns is 1; } } |
| | | 26 | | /// <summary>Determines if the matrix is a 2 component vector.</summary> |
| | 0 | 27 | | public bool Is2x1 => _rows is 2 && _columns is 1; |
| | | 28 | | /// <summary>Determines if the matrix is a 3 component vector.</summary> |
| | 0 | 29 | | public bool Is3x1 => _rows is 3 && _columns is 1; |
| | | 30 | | /// <summary>Determines if the matrix is a 4 component vector.</summary> |
| | 0 | 31 | | public bool Is4x1 => _rows is 4 && _columns is 1; |
| | | 32 | | /// <summary>Determines if the matrix is a 2 square matrix.</summary> |
| | 0 | 33 | | public bool Is2x2 => _rows is 2 && _columns is 2; |
| | | 34 | | /// <summary>Determines if the matrix is a 3 square matrix.</summary> |
| | 0 | 35 | | public bool Is3x3 => _rows is 3 && _columns is 3; |
| | | 36 | | /// <summary>Determines if the matrix is a 4 square matrix.</summary> |
| | 0 | 37 | | public bool Is4x4 => _rows is 4 && _columns is 4; |
| | | 38 | | |
| | | 39 | | /// <summary>Standard row-major matrix indexing.</summary> |
| | | 40 | | /// <param name="row">The row index.</param> |
| | | 41 | | /// <param name="column">The column index.</param> |
| | | 42 | | /// <returns>The value at the given indeces.</returns> |
| | | 43 | | public T this[int row, int column] |
| | | 44 | | { |
| | | 45 | | get |
| | 76 | 46 | | { |
| | 76 | 47 | | if (row < 0 || row >= Rows) |
| | 0 | 48 | | { |
| | 0 | 49 | | throw new ArgumentOutOfRangeException(nameof(row), row, "!(" + nameof(row) + " >= 0) || !(" + nameof(row) + " < |
| | | 50 | | } |
| | 76 | 51 | | if (column < 0 || column >= Columns) |
| | 0 | 52 | | { |
| | 0 | 53 | | throw new ArgumentOutOfRangeException(nameof(column), row, "!(" + nameof(column) + " >= 0) || !(" + nameof(colum |
| | | 54 | | } |
| | 76 | 55 | | return Get(row, column); |
| | 76 | 56 | | } |
| | | 57 | | set |
| | 45 | 58 | | { |
| | 45 | 59 | | if (row < 0 || row >= Rows) |
| | 0 | 60 | | { |
| | 0 | 61 | | throw new ArgumentOutOfRangeException(nameof(row), row, "!(" + nameof(row) + " >= 0) || !(" + nameof(row) + " < |
| | | 62 | | } |
| | 45 | 63 | | if (column < 0 || column >= Columns) |
| | 0 | 64 | | { |
| | 0 | 65 | | throw new ArgumentOutOfRangeException(nameof(column), row, "!(" + nameof(column) + " >= 0) || !(" + nameof(colum |
| | | 66 | | } |
| | 45 | 67 | | Set(row, column, value); |
| | 45 | 68 | | } |
| | | 69 | | } |
| | | 70 | | |
| | | 71 | | /// <summary>Indexing of the flattened array representing the matrix.</summary> |
| | | 72 | | /// <param name="flatIndex">The flattened index of the matrix.</param> |
| | | 73 | | /// <returns>The value at the given flattened index.</returns> |
| | | 74 | | public T this[int flatIndex] |
| | | 75 | | { |
| | | 76 | | get |
| | 0 | 77 | | { |
| | 0 | 78 | | if (flatIndex < 0 || _matrix.Length < flatIndex) |
| | 0 | 79 | | { |
| | 0 | 80 | | throw new ArgumentOutOfRangeException(nameof(flatIndex), flatIndex, "!(" + nameof(flatIndex) + " >= 0) || !(" + |
| | | 81 | | } |
| | 0 | 82 | | return _matrix[flatIndex]; |
| | 0 | 83 | | } |
| | | 84 | | set |
| | 9 | 85 | | { |
| | 9 | 86 | | if (flatIndex < 0 || _matrix.Length < flatIndex) |
| | 0 | 87 | | { |
| | 0 | 88 | | throw new ArgumentOutOfRangeException(nameof(flatIndex), flatIndex, "!(" + nameof(flatIndex) + " >= 0) || !(" + |
| | | 89 | | } |
| | 9 | 90 | | _matrix[flatIndex] = value; |
| | 9 | 91 | | } |
| | | 92 | | } |
| | | 93 | | |
| | | 94 | | #endregion |
| | | 95 | | |
| | | 96 | | #region Debugger Properties |
| | | 97 | | |
| | | 98 | | internal string? DebuggerString |
| | | 99 | | { |
| | | 100 | | get |
| | 0 | 101 | | { |
| | 0 | 102 | | if (Rows < 5 && Columns < 5) |
| | 0 | 103 | | { |
| | 0 | 104 | | StringBuilder stringBuilder = new(); |
| | 0 | 105 | | stringBuilder.Append('['); |
| | 0 | 106 | | for (int i = 0; i < Rows; i++) |
| | 0 | 107 | | { |
| | 0 | 108 | | stringBuilder.Append('['); |
| | 0 | 109 | | for (int j = 0; j < Columns; j++) |
| | 0 | 110 | | { |
| | 0 | 111 | | stringBuilder.Append(Get(i, j)); |
| | 0 | 112 | | if (j < Columns - 1) |
| | 0 | 113 | | { |
| | 0 | 114 | | stringBuilder.Append(','); |
| | 0 | 115 | | } |
| | 0 | 116 | | } |
| | 0 | 117 | | stringBuilder.Append(']'); |
| | 0 | 118 | | } |
| | 0 | 119 | | stringBuilder.Append(']'); |
| | 0 | 120 | | return stringBuilder.ToString(); |
| | | 121 | | } |
| | 0 | 122 | | return ToString(); |
| | 0 | 123 | | } |
| | | 124 | | } |
| | | 125 | | |
| | | 126 | | #endregion |
| | | 127 | | |
| | | 128 | | #region Constructors |
| | | 129 | | |
| | 57 | 130 | | internal Matrix(int rows, int columns, int length) |
| | 57 | 131 | | { |
| | 57 | 132 | | _matrix = new T[length]; |
| | 57 | 133 | | _rows = rows; |
| | 57 | 134 | | _columns = columns; |
| | 57 | 135 | | } |
| | | 136 | | |
| | | 137 | | /// <summary>Constructs a new default matrix of the given dimensions.</summary> |
| | | 138 | | /// <param name="rows">The number of row dimensions.</param> |
| | | 139 | | /// <param name="columns">The number of column dimensions.</param> |
| | 284 | 140 | | public Matrix(int rows, int columns) |
| | 284 | 141 | | { |
| | 284 | 142 | | if (rows < 1) |
| | 0 | 143 | | { |
| | 0 | 144 | | throw new ArgumentOutOfRangeException(nameof(rows), rows, "!(rows > 0)"); |
| | | 145 | | } |
| | 284 | 146 | | if (columns < 1) |
| | 0 | 147 | | { |
| | 0 | 148 | | throw new ArgumentOutOfRangeException(nameof(columns), columns, "!(columns > 0)"); |
| | | 149 | | } |
| | 284 | 150 | | _matrix = new T[rows * columns]; |
| | 284 | 151 | | _rows = rows; |
| | 284 | 152 | | _columns = columns; |
| | 284 | 153 | | } |
| | | 154 | | |
| | | 155 | | /// <summary>Constructs a new matrix and initializes it via function.</summary> |
| | | 156 | | /// <param name="rows">The number of rows to construct.</param> |
| | | 157 | | /// <param name="columns">The number of columns to construct.</param> |
| | | 158 | | /// <param name="function">The initialization function.</param> |
| | 237 | 159 | | public Matrix(int rows, int columns, Func<int, int, T> function) : this(rows, columns) |
| | 237 | 160 | | { |
| | 237 | 161 | | Format(this, function); |
| | 237 | 162 | | } |
| | | 163 | | |
| | | 164 | | /// <summary>Constructs a new matrix and initializes it via function.</summary> |
| | | 165 | | /// <param name="rows">The number of rows to construct.</param> |
| | | 166 | | /// <param name="columns">The number of columns to construct.</param> |
| | | 167 | | /// <param name="function">The initialization function.</param> |
| | 0 | 168 | | public Matrix(int rows, int columns, Func<int, T> function) : this(rows, columns) |
| | 0 | 169 | | { |
| | 0 | 170 | | Format(this, function); |
| | 0 | 171 | | } |
| | | 172 | | |
| | | 173 | | /// <summary> |
| | | 174 | | /// Creates a new matrix using ROW MAJOR ORDER. The data will be referenced to, so |
| | | 175 | | /// changes to it will modify the constructed matrix. |
| | | 176 | | /// </summary> |
| | | 177 | | /// <param name="rows">The number of rows to construct.</param> |
| | | 178 | | /// <param name="columns">The number of columns to construct.</param> |
| | | 179 | | /// <param name="data">The data of the matrix in ROW MAJOR ORDER.</param> |
| | 0 | 180 | | public Matrix(int rows, int columns, params T[] data) : this(rows, columns) |
| | 0 | 181 | | { |
| | 0 | 182 | | if (data is null) throw new ArgumentNullException(nameof(data)); |
| | 0 | 183 | | if (sourceof(data.Length != rows * columns, out string c1)) throw new ArgumentException(c1); |
| | 0 | 184 | | _rows = rows; |
| | 0 | 185 | | _columns = columns; |
| | 0 | 186 | | _matrix = data; |
| | 0 | 187 | | } |
| | | 188 | | |
| | 16 | 189 | | internal Matrix(Matrix<T> matrix) |
| | 16 | 190 | | { |
| | 16 | 191 | | _rows = matrix._rows; |
| | 16 | 192 | | _columns = matrix.Columns; |
| | 16 | 193 | | _matrix = (T[])matrix._matrix.Clone(); |
| | 16 | 194 | | } |
| | | 195 | | |
| | 0 | 196 | | internal Matrix(Vector<T> vector) |
| | 0 | 197 | | { |
| | 0 | 198 | | _rows = vector.Dimensions; |
| | 0 | 199 | | _columns = 1; |
| | 0 | 200 | | _matrix = (T[])vector._vector.Clone(); |
| | 0 | 201 | | } |
| | | 202 | | |
| | | 203 | | #endregion |
| | | 204 | | |
| | | 205 | | #region Factories |
| | | 206 | | |
| | | 207 | | /// <summary>Constructs a new zero-matrix of the given dimensions.</summary> |
| | | 208 | | /// <param name="rows">The number of rows of the matrix.</param> |
| | | 209 | | /// <param name="columns">The number of columns of the matrix.</param> |
| | | 210 | | /// <returns>The newly constructed zero-matrix.</returns> |
| | | 211 | | public static Matrix<T> FactoryZero(int rows, int columns) |
| | 6 | 212 | | { |
| | 6 | 213 | | if (rows < 1) |
| | 1 | 214 | | { |
| | 1 | 215 | | throw new ArgumentOutOfRangeException(nameof(rows), rows, "!(" + nameof(rows) + " > 0)"); |
| | | 216 | | } |
| | 5 | 217 | | if (columns < 1) |
| | 1 | 218 | | { |
| | 1 | 219 | | throw new ArgumentOutOfRangeException(nameof(columns), columns, "!(" + nameof(columns) + " > 0)"); |
| | | 220 | | } |
| | 4 | 221 | | return FactoryZeroImplementation(rows, columns); |
| | 4 | 222 | | } |
| | | 223 | | |
| | 1 | 224 | | internal static Func<int, int, Matrix<T>> FactoryZeroImplementation = (rows, columns) => |
| | 1 | 225 | | { |
| | 1 | 226 | | if (Equate(default, Constant<T>.Zero)) |
| | 1 | 227 | | { |
| | 5 | 228 | | FactoryZeroImplementation = (ROWS, COLUMNS) => new Matrix<T>(ROWS, COLUMNS); |
| | 1 | 229 | | } |
| | 1 | 230 | | else |
| | 0 | 231 | | { |
| | 0 | 232 | | FactoryZeroImplementation = (ROWS, COLUMNS) => |
| | 0 | 233 | | { |
| | 0 | 234 | | Matrix<T> matrix = new(ROWS, COLUMNS); |
| | 0 | 235 | | matrix._matrix.Format(Constant<T>.Zero); |
| | 0 | 236 | | return matrix; |
| | 0 | 237 | | }; |
| | 0 | 238 | | } |
| | 1 | 239 | | return FactoryZeroImplementation(rows, columns); |
| | 2 | 240 | | }; |
| | | 241 | | |
| | | 242 | | /// <summary>Constructs a new identity-matrix of the given dimensions.</summary> |
| | | 243 | | /// <param name="rows">The number of rows of the matrix.</param> |
| | | 244 | | /// <param name="columns">The number of columns of the matrix.</param> |
| | | 245 | | /// <returns>The newly constructed identity-matrix.</returns> |
| | | 246 | | public static Matrix<T> FactoryIdentity(int rows, int columns) |
| | 6 | 247 | | { |
| | 6 | 248 | | if (rows < 1) |
| | 1 | 249 | | { |
| | 1 | 250 | | throw new ArgumentOutOfRangeException(nameof(rows), rows, "!(" + nameof(rows) + " > 0)"); |
| | | 251 | | } |
| | 5 | 252 | | if (columns < 1) |
| | 1 | 253 | | { |
| | 1 | 254 | | throw new ArgumentOutOfRangeException(nameof(columns), columns, "!(" + nameof(columns) + " > 0)"); |
| | | 255 | | } |
| | 4 | 256 | | return FactoryIdentityImplementation(rows, columns); |
| | 4 | 257 | | } |
| | | 258 | | |
| | 1 | 259 | | internal static Func<int, int, Matrix<T>> FactoryIdentityImplementation = (rows, columns) => |
| | 1 | 260 | | { |
| | 1 | 261 | | if (Equate(default, Constant<T>.Zero)) |
| | 1 | 262 | | { |
| | 1 | 263 | | FactoryIdentityImplementation = (ROWS, COLUMNS) => |
| | 4 | 264 | | { |
| | 4 | 265 | | Matrix<T> matrix = new(ROWS, COLUMNS); |
| | 4 | 266 | | T[] MATRIX = matrix._matrix; |
| | 4 | 267 | | int minimum = Math.Min(ROWS, COLUMNS); |
| | 28 | 268 | | for (int i = 0; i < minimum; i++) |
| | 10 | 269 | | { |
| | 10 | 270 | | MATRIX[i * COLUMNS + i] = Constant<T>.One; |
| | 10 | 271 | | } |
| | 4 | 272 | | return matrix; |
| | 5 | 273 | | }; |
| | 1 | 274 | | } |
| | 1 | 275 | | else |
| | 0 | 276 | | { |
| | 0 | 277 | | FactoryIdentityImplementation = (ROWS, COLUMNS) => |
| | 0 | 278 | | { |
| | 0 | 279 | | Matrix<T> matrix = new(ROWS, COLUMNS); |
| | 0 | 280 | | Format(matrix, (x, y) => x == y ? Constant<T>.One : Constant<T>.Zero); |
| | 0 | 281 | | return matrix; |
| | 0 | 282 | | }; |
| | 0 | 283 | | } |
| | 1 | 284 | | return FactoryIdentityImplementation(rows, columns); |
| | 2 | 285 | | }; |
| | | 286 | | |
| | | 287 | | /// <summary>Constructs a new matrix where every entry is the same uniform value.</summary> |
| | | 288 | | /// <param name="rows">The number of rows of the matrix.</param> |
| | | 289 | | /// <param name="columns">The number of columns of the matrix.</param> |
| | | 290 | | /// <param name="value">The value to assign every spot in the matrix.</param> |
| | | 291 | | /// <returns>The newly constructed matrix filled with the uniform value.</returns> |
| | | 292 | | public static Matrix<T> FactoryUniform(int rows, int columns, T value) |
| | 0 | 293 | | { |
| | 0 | 294 | | if (rows < 1) |
| | 0 | 295 | | { |
| | 0 | 296 | | throw new ArgumentOutOfRangeException(nameof(rows), rows, "!(" + nameof(rows) + " > 0)"); |
| | | 297 | | } |
| | 0 | 298 | | if (columns < 1) |
| | 0 | 299 | | { |
| | 0 | 300 | | throw new ArgumentOutOfRangeException(nameof(columns), columns, "!(" + nameof(columns) + " > 0)"); |
| | | 301 | | } |
| | 0 | 302 | | Matrix<T> matrix = new(rows, columns); |
| | 0 | 303 | | matrix._matrix.Format(value); |
| | 0 | 304 | | return matrix; |
| | 0 | 305 | | } |
| | | 306 | | |
| | | 307 | | #endregion |
| | | 308 | | |
| | | 309 | | #region Mathematics |
| | | 310 | | |
| | | 311 | | #region RowMultiplication |
| | | 312 | | |
| | | 313 | | internal static void RowMultiplication(Matrix<T> matrix, int row, T scalar) |
| | 0 | 314 | | { |
| | 0 | 315 | | int columns = matrix.Columns; |
| | 0 | 316 | | for (int i = 0; i < columns; i++) |
| | 0 | 317 | | { |
| | 0 | 318 | | matrix.Set(row, i, Multiplication(matrix.Get(row, i), scalar)); |
| | 0 | 319 | | } |
| | 0 | 320 | | } |
| | | 321 | | |
| | | 322 | | #endregion |
| | | 323 | | |
| | | 324 | | #region RowAddition |
| | | 325 | | |
| | | 326 | | internal static void RowAddition(Matrix<T> matrix, int target, int second, T scalar) |
| | 0 | 327 | | { |
| | 0 | 328 | | int columns = matrix.Columns; |
| | 0 | 329 | | for (int i = 0; i < columns; i++) |
| | 0 | 330 | | { |
| | 0 | 331 | | matrix.Set(target, i, Addition(matrix.Get(target, i), Multiplication(matrix.Get(second, i), scalar))); |
| | 0 | 332 | | } |
| | 0 | 333 | | } |
| | | 334 | | |
| | | 335 | | #endregion |
| | | 336 | | |
| | | 337 | | #region SwapRows |
| | | 338 | | |
| | | 339 | | internal static void SwapRows(Matrix<T> matrix, int row1, int row2) => |
| | 0 | 340 | | SwapRows(matrix, row1, row2, 0, matrix.Columns - 1); |
| | | 341 | | |
| | | 342 | | internal static void SwapRows(Matrix<T> matrix, int row1, int row2, int columnStart, int columnEnd) |
| | 54 | 343 | | { |
| | 340 | 344 | | for (int i = columnStart; i <= columnEnd; i++) |
| | 116 | 345 | | { |
| | 116 | 346 | | T temp = matrix.Get(row1, i); |
| | 116 | 347 | | matrix.Set(row1, i, matrix.Get(row2, i)); |
| | 116 | 348 | | matrix.Set(row2, i, temp); |
| | 116 | 349 | | } |
| | 54 | 350 | | } |
| | | 351 | | |
| | | 352 | | #endregion |
| | | 353 | | |
| | | 354 | | #region GetCofactor |
| | | 355 | | |
| | | 356 | | internal static void GetCofactor(Matrix<T> a, Matrix<T> temp, int p, int q, int n) |
| | 98 | 357 | | { |
| | 196 | 358 | | int i = 0, j = 0; |
| | 720 | 359 | | for (int row = 0; row < n; row++) |
| | 262 | 360 | | { |
| | 1968 | 361 | | for (int col = 0; col < n; col++) |
| | 722 | 362 | | { |
| | 722 | 363 | | if (row != p && col != q) |
| | 296 | 364 | | { |
| | 296 | 365 | | temp.Set(i, j++, a.Get(row, col)); |
| | 296 | 366 | | if (j == n - 1) |
| | 164 | 367 | | { |
| | 164 | 368 | | j = 0; |
| | 164 | 369 | | i++; |
| | 164 | 370 | | } |
| | 296 | 371 | | } |
| | 722 | 372 | | } |
| | 262 | 373 | | } |
| | 98 | 374 | | } |
| | | 375 | | |
| | | 376 | | #endregion |
| | | 377 | | |
| | | 378 | | #region GetDeterminant Gaussian elimination |
| | | 379 | | |
| | | 380 | | #region MatrixElementFraction |
| | | 381 | | /// <summary> |
| | | 382 | | /// Used to avoid issues when 1/2 + 1/2 = 0 + 0 = 0 instead of 1 for types, where division results in precision loss |
| | | 383 | | /// </summary> |
| | | 384 | | #pragma warning disable CS0660 // Type defines operator == or operator != but does not override Object.Equals(object o) |
| | | 385 | | #pragma warning disable CS0661 // Type defines operator == or operator != but does not override Object.GetHashCode() |
| | | 386 | | private sealed class MatrixElementFraction |
| | | 387 | | #pragma warning restore CS0661 // Type defines operator == or operator != but does not override Object.GetHashCode() |
| | | 388 | | #pragma warning restore CS0660 // Type defines operator == or operator != but does not override Object.Equals(object o) |
| | | 389 | | { |
| | | 390 | | private readonly T Numerator; |
| | | 391 | | private readonly T Denominator; |
| | | 392 | | |
| | 620 | 393 | | internal MatrixElementFraction(T value) |
| | 620 | 394 | | { |
| | 620 | 395 | | Numerator = value; |
| | 620 | 396 | | Denominator = Constant<T>.One; |
| | 620 | 397 | | } |
| | | 398 | | |
| | 692 | 399 | | internal MatrixElementFraction(T num, T den) |
| | 692 | 400 | | { |
| | 692 | 401 | | Numerator = num; |
| | 692 | 402 | | Denominator = den; |
| | 692 | 403 | | } |
| | | 404 | | |
| | 64 | 405 | | public T Value => Division(Numerator, Denominator); |
| | | 406 | | |
| | | 407 | | // a / b + c / d = (a * d + b * c) / (b * d) |
| | | 408 | | public static MatrixElementFraction operator +(MatrixElementFraction a, MatrixElementFraction b) |
| | 88 | 409 | | => new(Addition( |
| | 88 | 410 | | Multiplication(a.Numerator, b.Denominator), |
| | 88 | 411 | | Multiplication(a.Denominator, b.Numerator) |
| | 88 | 412 | | ), Multiplication(a.Denominator, b.Denominator)); |
| | | 413 | | |
| | | 414 | | public static MatrixElementFraction operator *(MatrixElementFraction a, MatrixElementFraction b) |
| | 222 | 415 | | => new(Multiplication(a.Numerator, b.Numerator), Multiplication(a.Denominator, b.Denominator)); |
| | | 416 | | |
| | | 417 | | public static MatrixElementFraction operator -(MatrixElementFraction a) |
| | 142 | 418 | | => new(Multiplication(a.Numerator, Constant<T>.NegativeOne), a.Denominator); |
| | | 419 | | |
| | | 420 | | public static MatrixElementFraction operator -(MatrixElementFraction a, MatrixElementFraction b) |
| | 88 | 421 | | => a + (-b); |
| | | 422 | | |
| | | 423 | | // a / b / (d / c) = (a * c) / (b * d) |
| | | 424 | | public static MatrixElementFraction operator /(MatrixElementFraction a, MatrixElementFraction b) |
| | 88 | 425 | | => new(Multiplication(a.Numerator, b.Denominator), Multiplication(a.Denominator, b.Numerator)); |
| | | 426 | | |
| | | 427 | | public static bool operator <(MatrixElementFraction a, MatrixElementFraction b) |
| | 76 | 428 | | { |
| | 76 | 429 | | var c = Multiplication(a.Numerator, b.Denominator); |
| | 76 | 430 | | var d = Multiplication(b.Numerator, a.Denominator); |
| | 76 | 431 | | if (IsPositive(Multiplication(a.Denominator, b.Denominator))) |
| | 76 | 432 | | return LessThan(c, d); |
| | | 433 | | else |
| | 0 | 434 | | return !LessThan(c, d); |
| | 76 | 435 | | } |
| | | 436 | | |
| | | 437 | | public static bool operator ==(MatrixElementFraction a, MatrixElementFraction b) => |
| | 64 | 438 | | Equate(a.Numerator, b.Numerator) && |
| | 64 | 439 | | Equate(a.Denominator, b.Denominator); |
| | | 440 | | public static bool operator !=(MatrixElementFraction a, MatrixElementFraction b) |
| | 0 | 441 | | => !(a == b); |
| | | 442 | | |
| | | 443 | | public static bool operator >(MatrixElementFraction a, MatrixElementFraction b) |
| | 76 | 444 | | => a >= b && !(a == b); |
| | | 445 | | |
| | | 446 | | public static bool operator >=(MatrixElementFraction a, MatrixElementFraction b) |
| | 76 | 447 | | => !(a < b); |
| | | 448 | | |
| | | 449 | | public static bool operator <=(MatrixElementFraction a, MatrixElementFraction b) |
| | 0 | 450 | | => a < b || a == b; |
| | | 451 | | |
| | | 452 | | public static explicit operator MatrixElementFraction(int val) |
| | 0 | 453 | | => new(Convert<int, T>(val)); |
| | | 454 | | |
| | | 455 | | public MatrixElementFraction Abs() |
| | 152 | 456 | | => new(AbsoluteValue(Numerator), AbsoluteValue(Denominator)); |
| | | 457 | | |
| | 64 | 458 | | public bool IsDividedByZero => Equate(Denominator, Constant<T>.Zero); |
| | | 459 | | } |
| | | 460 | | #endregion |
| | | 461 | | |
| | | 462 | | // Reference: https://codereview.stackexchange.com/questions/204135/determinant-using-gauss-elimination |
| | | 463 | | internal static T GetDeterminantGaussian(Matrix<T> matrix, int n) |
| | 64 | 464 | | { |
| | | 465 | | // Note: the reasoning behind using MatrixElementFraction is to account for |
| | | 466 | | // rounding errors such as 2.00...01 instead of 2 |
| | | 467 | | |
| | 64 | 468 | | if (n is 1) |
| | 0 | 469 | | { |
| | 0 | 470 | | return matrix.Get(0, 0); |
| | | 471 | | } |
| | 64 | 472 | | MatrixElementFraction determinant = new(Constant<T>.One); |
| | 620 | 473 | | Matrix<MatrixElementFraction> fractioned = new(matrix.Rows, matrix.Columns, (r, c) => new(matrix.Get(r, c))); |
| | 396 | 474 | | for (int i = 0; i < n; i++) |
| | 134 | 475 | | { |
| | 134 | 476 | | var pivotElement = fractioned.Get(i, i); |
| | 134 | 477 | | var pivotRow = i; |
| | 420 | 478 | | for (int row = i + 1; row < n; ++row) |
| | 76 | 479 | | { |
| | 76 | 480 | | if (fractioned[row, i].Abs() > pivotElement.Abs()) |
| | 58 | 481 | | { |
| | 58 | 482 | | pivotElement = fractioned.Get(row, i); |
| | 58 | 483 | | pivotRow = row; |
| | 58 | 484 | | } |
| | 76 | 485 | | } |
| | 134 | 486 | | if (pivotElement.Equals(Constant<T>.Zero)) |
| | 0 | 487 | | { |
| | 0 | 488 | | return Constant<T>.Zero; |
| | | 489 | | } |
| | 134 | 490 | | if (pivotRow != i) |
| | 54 | 491 | | { |
| | 54 | 492 | | Matrix<MatrixElementFraction>.SwapRows(fractioned, i, pivotRow, 0, n - 1); |
| | 54 | 493 | | determinant = Negation(determinant); |
| | 54 | 494 | | } |
| | 134 | 495 | | determinant *= pivotElement; |
| | 420 | 496 | | for (int row = i + 1; row < n; row++) |
| | 76 | 497 | | { |
| | 328 | 498 | | for (int column = i + 1; column < n; column++) |
| | 88 | 499 | | { |
| | | 500 | | // reference: matrix[row][column] -= matrix[row][i] * matrix[i][column] / pivotElement; |
| | 88 | 501 | | fractioned.Set(row, column, |
| | 88 | 502 | | // D - A * B / C |
| | 88 | 503 | | D_subtract_A_multiply_B_divide_C<MatrixElementFraction>.Function( |
| | 88 | 504 | | /* A */ fractioned.Get(row, i), |
| | 88 | 505 | | /* B */ fractioned.Get(i, column), |
| | 88 | 506 | | /* C */ pivotElement, |
| | 88 | 507 | | /* D */ fractioned.Get(row, column))); |
| | 88 | 508 | | } |
| | 76 | 509 | | } |
| | 134 | 510 | | } |
| | | 511 | | #warning TODO: should we return zero if determinant's denominator is zero? |
| | 64 | 512 | | return determinant.IsDividedByZero ? Constant<T>.Zero : determinant.Value; |
| | 64 | 513 | | } |
| | | 514 | | |
| | | 515 | | #endregion |
| | | 516 | | |
| | | 517 | | #region GetDeterminant Laplace method |
| | | 518 | | |
| | | 519 | | internal static T GetDeterminantLaplace(Matrix<T> a, int n) |
| | 52 | 520 | | { |
| | 52 | 521 | | T determinant = Constant<T>.Zero; |
| | 52 | 522 | | if (n is 1) |
| | 32 | 523 | | { |
| | 32 | 524 | | return a.Get(0, 0); |
| | | 525 | | } |
| | 20 | 526 | | Matrix<T> temp = new(n, n); |
| | 20 | 527 | | T sign = Constant<T>.One; |
| | 128 | 528 | | for (int f = 0; f < n; f++) |
| | 44 | 529 | | { |
| | 44 | 530 | | GetCofactor(a, temp, 0, f, n); |
| | 44 | 531 | | determinant = |
| | 44 | 532 | | Addition(determinant, |
| | 44 | 533 | | Multiplication(sign, |
| | 44 | 534 | | Multiplication(a.Get(0, f), |
| | 44 | 535 | | GetDeterminantLaplace(temp, n - 1)))); |
| | 44 | 536 | | sign = Negation(sign); |
| | 44 | 537 | | } |
| | 20 | 538 | | return determinant; |
| | 52 | 539 | | } |
| | | 540 | | |
| | | 541 | | #endregion |
| | | 542 | | |
| | | 543 | | #region IsSymetric |
| | | 544 | | |
| | | 545 | | /// <summary>Determines if the matrix is symetric.</summary> |
| | | 546 | | /// <param name="a">The matrix to determine if symetric.</param> |
| | | 547 | | /// <returns>True if the matrix is symetric; false if not.</returns> |
| | | 548 | | public static bool GetIsSymetric(Matrix<T> a) |
| | 0 | 549 | | { |
| | 0 | 550 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 0 | 551 | | int rows = a._rows; |
| | 0 | 552 | | if (rows != a._columns) |
| | 0 | 553 | | { |
| | 0 | 554 | | return false; |
| | | 555 | | } |
| | 0 | 556 | | T[] A = a._matrix; |
| | 0 | 557 | | for (int row = 0; row < rows; row++) |
| | 0 | 558 | | { |
| | 0 | 559 | | for (int column = row + 1; column < rows; column++) |
| | 0 | 560 | | { |
| | 0 | 561 | | if (Inequate(A[row * rows + column], A[column * rows + row])) |
| | 0 | 562 | | { |
| | 0 | 563 | | return false; |
| | | 564 | | } |
| | 0 | 565 | | } |
| | 0 | 566 | | } |
| | 0 | 567 | | return true; |
| | 0 | 568 | | } |
| | | 569 | | |
| | | 570 | | /// <summary>Determines if the matrix is symetric.</summary> |
| | | 571 | | /// <returns>True if the matrix is symetric; false if not.</returns> |
| | | 572 | | public bool IsSymetric |
| | | 573 | | { |
| | | 574 | | get |
| | 0 | 575 | | { |
| | 0 | 576 | | return GetIsSymetric(this); |
| | 0 | 577 | | } |
| | | 578 | | } |
| | | 579 | | |
| | | 580 | | #endregion |
| | | 581 | | |
| | | 582 | | #region Negate |
| | | 583 | | |
| | | 584 | | /// <summary>Negates all the values in a matrix.</summary> |
| | | 585 | | /// <param name="a">The matrix to have its values negated.</param> |
| | | 586 | | /// <param name="b">The resulting matrix after the negation.</param> |
| | | 587 | | public static void Negate(Matrix<T> a, ref Matrix<T>? b) |
| | 4 | 588 | | { |
| | 4 | 589 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 4 | 590 | | int Length = a._matrix.Length; |
| | 4 | 591 | | T[] A = a._matrix; |
| | | 592 | | T[] B; |
| | 4 | 593 | | if (b is not null && b._matrix.Length == Length) |
| | 0 | 594 | | { |
| | 0 | 595 | | B = b._matrix; |
| | 0 | 596 | | b._rows = a._rows; |
| | 0 | 597 | | b._columns = a._columns; |
| | 0 | 598 | | } |
| | | 599 | | else |
| | 4 | 600 | | { |
| | 4 | 601 | | b = new Matrix<T>(a._rows, a._columns, Length); |
| | 4 | 602 | | B = b._matrix; |
| | 4 | 603 | | } |
| | 80 | 604 | | for (int i = 0; i < Length; i++) |
| | 36 | 605 | | { |
| | 36 | 606 | | B[i] = Negation(A[i]); |
| | 36 | 607 | | } |
| | 4 | 608 | | } |
| | | 609 | | |
| | | 610 | | /// <summary>Negates all the values in a matrix.</summary> |
| | | 611 | | /// <param name="a">The matrix to have its values negated.</param> |
| | | 612 | | /// <returns>The resulting matrix after the negation.</returns> |
| | | 613 | | public static Matrix<T> Negate(Matrix<T> a) |
| | 4 | 614 | | { |
| | 4 | 615 | | Matrix<T>? b = null; |
| | 4 | 616 | | Negate(a, ref b); |
| | 4 | 617 | | return b!; |
| | 4 | 618 | | } |
| | | 619 | | |
| | | 620 | | /// <summary>Negates all the values in a matrix.</summary> |
| | | 621 | | /// <param name="a">The matrix to have its values negated.</param> |
| | | 622 | | /// <returns>The resulting matrix after the negation.</returns> |
| | | 623 | | public static Matrix<T> operator -(Matrix<T> a) |
| | 4 | 624 | | { |
| | 4 | 625 | | return Negate(a); |
| | 4 | 626 | | } |
| | | 627 | | |
| | | 628 | | /// <summary>Negates all the values in a matrix.</summary> |
| | | 629 | | /// <param name="b">The resulting matrix after the negation.</param> |
| | | 630 | | public void Negate(ref Matrix<T>? b) |
| | 0 | 631 | | { |
| | 0 | 632 | | Negate(this, ref b); |
| | 0 | 633 | | } |
| | | 634 | | |
| | | 635 | | /// <summary>Negates all the values in this matrix.</summary> |
| | | 636 | | /// <returns>The resulting matrix after the negation.</returns> |
| | | 637 | | public Matrix<T> Negate() |
| | 0 | 638 | | { |
| | 0 | 639 | | return -this; |
| | 0 | 640 | | } |
| | | 641 | | |
| | | 642 | | #endregion |
| | | 643 | | |
| | | 644 | | #region Add |
| | | 645 | | |
| | | 646 | | /// <summary>Does standard addition of two matrices.</summary> |
| | | 647 | | /// <param name="a">The left matrix of the addition.</param> |
| | | 648 | | /// <param name="b">The right matrix of the addition.</param> |
| | | 649 | | /// <param name="c">The resulting matrix after the addition.</param> |
| | | 650 | | public static void Add(Matrix<T> a, Matrix<T> b, ref Matrix<T>? c) |
| | 5 | 651 | | { |
| | 5 | 652 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 5 | 653 | | if (b is null) throw new ArgumentNullException(nameof(b)); |
| | 6 | 654 | | if (sourceof(a._rows != b._rows || a._columns != b._columns, out string c1)) throw new ArgumentException(c1); |
| | 4 | 655 | | int Length = a._matrix.Length; |
| | 4 | 656 | | T[] A = a._matrix; |
| | 4 | 657 | | T[] B = b._matrix; |
| | | 658 | | T[] C; |
| | 4 | 659 | | if (c is not null && c._matrix.Length == Length) |
| | 0 | 660 | | { |
| | 0 | 661 | | C = c._matrix; |
| | 0 | 662 | | c._rows = a._rows; |
| | 0 | 663 | | c._columns = a._columns; |
| | 0 | 664 | | } |
| | | 665 | | else |
| | 4 | 666 | | { |
| | 4 | 667 | | c = new Matrix<T>(a._rows, a._columns, Length); |
| | 4 | 668 | | C = c._matrix; |
| | 4 | 669 | | } |
| | 80 | 670 | | for (int i = 0; i < Length; i++) |
| | 36 | 671 | | { |
| | 36 | 672 | | C[i] = Addition(A[i], B[i]); |
| | 36 | 673 | | } |
| | 4 | 674 | | } |
| | | 675 | | |
| | | 676 | | /// <summary>Does standard addition of two matrices.</summary> |
| | | 677 | | /// <param name="a">The left matrix of the addition.</param> |
| | | 678 | | /// <param name="b">The right matrix of the addition.</param> |
| | | 679 | | /// <returns>The resulting matrix after the addition.</returns> |
| | | 680 | | public static Matrix<T> Add(Matrix<T> a, Matrix<T> b) |
| | 5 | 681 | | { |
| | 5 | 682 | | Matrix<T>? c = null; |
| | 5 | 683 | | Add(a, b, ref c); |
| | 4 | 684 | | return c!; |
| | 4 | 685 | | } |
| | | 686 | | |
| | | 687 | | /// <summary>Does a standard matrix addition.</summary> |
| | | 688 | | /// <param name="a">The left matrix of the addition.</param> |
| | | 689 | | /// <param name="b">The right matrix of the addition.</param> |
| | | 690 | | /// <returns>The resulting matrix after the addition.</returns> |
| | | 691 | | public static Matrix<T> operator +(Matrix<T> a, Matrix<T> b) |
| | 5 | 692 | | { |
| | 5 | 693 | | return Add(a, b); |
| | 4 | 694 | | } |
| | | 695 | | |
| | | 696 | | /// <summary>Does standard addition of two matrices.</summary> |
| | | 697 | | /// <param name="b">The right matrix of the addition.</param> |
| | | 698 | | /// <param name="c">The resulting matrix after the addition.</param> |
| | | 699 | | public void Add(Matrix<T> b, ref Matrix<T>? c) |
| | 0 | 700 | | { |
| | 0 | 701 | | Add(this, b, ref c); |
| | 0 | 702 | | } |
| | | 703 | | |
| | | 704 | | /// <summary>Does a standard matrix addition.</summary> |
| | | 705 | | /// <param name="b">The matrix to add to this matrix.</param> |
| | | 706 | | /// <returns>The resulting matrix after the addition.</returns> |
| | | 707 | | public Matrix<T> Add(Matrix<T> b) |
| | 0 | 708 | | { |
| | 0 | 709 | | return this + b; |
| | 0 | 710 | | } |
| | | 711 | | |
| | | 712 | | #endregion |
| | | 713 | | |
| | | 714 | | #region Subtract |
| | | 715 | | |
| | | 716 | | /// <summary>Does a standard matrix subtraction.</summary> |
| | | 717 | | /// <param name="a">The left matrix of the subtraction.</param> |
| | | 718 | | /// <param name="b">The right matrix of the subtraction.</param> |
| | | 719 | | /// <param name="c">The resulting matrix after the subtraction.</param> |
| | | 720 | | public static void Subtract(Matrix<T> a, Matrix<T> b, ref Matrix<T>? c) |
| | 5 | 721 | | { |
| | 5 | 722 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 5 | 723 | | if (b is null) throw new ArgumentNullException(nameof(b)); |
| | 6 | 724 | | if (sourceof(a._rows != b._rows || a._columns != b._columns, out string c1)) throw new ArgumentException(c1); |
| | 4 | 725 | | T[] A = a._matrix; |
| | 4 | 726 | | T[] B = b._matrix; |
| | 4 | 727 | | int Length = A.Length; |
| | | 728 | | T[] C; |
| | 4 | 729 | | if (c is not null && c._matrix.Length == Length) |
| | 0 | 730 | | { |
| | 0 | 731 | | C = c._matrix; |
| | 0 | 732 | | c._rows = a._rows; |
| | 0 | 733 | | c._columns = a._columns; |
| | 0 | 734 | | } |
| | | 735 | | else |
| | 4 | 736 | | { |
| | 4 | 737 | | c = new Matrix<T>(a._rows, a._columns, Length); |
| | 4 | 738 | | C = c._matrix; |
| | 4 | 739 | | } |
| | 80 | 740 | | for (int i = 0; i < Length; i++) |
| | 36 | 741 | | { |
| | 36 | 742 | | C[i] = Subtraction(A[i], B[i]); |
| | 36 | 743 | | } |
| | 4 | 744 | | } |
| | | 745 | | |
| | | 746 | | /// <summary>Does a standard matrix subtraction.</summary> |
| | | 747 | | /// <param name="a">The left matrix of the subtraction.</param> |
| | | 748 | | /// <param name="b">The right matrix of the subtraction.</param> |
| | | 749 | | /// <returns>The resulting matrix after the subtraction.</returns> |
| | | 750 | | public static Matrix<T> Subtract(Matrix<T> a, Matrix<T> b) |
| | 5 | 751 | | { |
| | 5 | 752 | | Matrix<T>? c = null; |
| | 5 | 753 | | Subtract(a, b, ref c); |
| | 4 | 754 | | return c!; |
| | 4 | 755 | | } |
| | | 756 | | |
| | | 757 | | /// <summary>Does a standard matrix subtraction.</summary> |
| | | 758 | | /// <param name="a">The left matrix of the subtraction.</param> |
| | | 759 | | /// <param name="b">The right matrix of the subtraction.</param> |
| | | 760 | | /// <returns>The resulting matrix after the subtraction.</returns> |
| | | 761 | | public static Matrix<T> operator -(Matrix<T> a, Matrix<T> b) |
| | 5 | 762 | | { |
| | 5 | 763 | | return Subtract(a, b); |
| | 4 | 764 | | } |
| | | 765 | | |
| | | 766 | | /// <summary>Does a standard matrix subtraction.</summary> |
| | | 767 | | /// <param name="b">The right matrix of the subtraction.</param> |
| | | 768 | | /// <param name="c">The resulting matrix after the subtraction.</param> |
| | | 769 | | public void Subtract(Matrix<T> b, ref Matrix<T>? c) |
| | 0 | 770 | | { |
| | 0 | 771 | | Subtract(this, b, ref c); |
| | 0 | 772 | | } |
| | | 773 | | |
| | | 774 | | /// <summary>Does a standard matrix subtraction.</summary> |
| | | 775 | | /// <param name="b">The right matrix of the subtraction.</param> |
| | | 776 | | /// <returns>The resulting matrix after the subtraction.</returns> |
| | | 777 | | public Matrix<T> Subtract(Matrix<T> b) |
| | 0 | 778 | | { |
| | 0 | 779 | | return this - b; |
| | 0 | 780 | | } |
| | | 781 | | |
| | | 782 | | #endregion |
| | | 783 | | |
| | | 784 | | #region Multiply (Matrix * Matrix) |
| | | 785 | | |
| | | 786 | | /// <summary>Does a standard (triple for looped) multiplication between matrices.</summary> |
| | | 787 | | /// <param name="a">The left matrix of the multiplication.</param> |
| | | 788 | | /// <param name="b">The right matrix of the multiplication.</param> |
| | | 789 | | /// <param name="c">The resulting matrix of the multiplication.</param> |
| | | 790 | | public static void Multiply(Matrix<T> a, Matrix<T> b, ref Matrix<T>? c) |
| | 18 | 791 | | { |
| | 18 | 792 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 18 | 793 | | if (b is null) throw new ArgumentNullException(nameof(b)); |
| | 19 | 794 | | if (sourceof(a._columns != b._rows, out string c1)) throw new ArgumentException(c1); |
| | | 795 | | |
| | 17 | 796 | | if (ReferenceEquals(a, b) && ReferenceEquals(a, c)) |
| | 0 | 797 | | { |
| | 0 | 798 | | Matrix<T> clone = a.Clone(); |
| | 0 | 799 | | a = clone; |
| | 0 | 800 | | b = clone; |
| | 0 | 801 | | } |
| | 17 | 802 | | else if (ReferenceEquals(a, c)) |
| | 8 | 803 | | { |
| | 8 | 804 | | a = a.Clone(); |
| | 8 | 805 | | } |
| | 9 | 806 | | else if (ReferenceEquals(b, c)) |
| | 0 | 807 | | { |
| | 0 | 808 | | b = b.Clone(); |
| | 0 | 809 | | } |
| | 17 | 810 | | int c_Rows = a._rows; |
| | 17 | 811 | | int a_Columns = a._columns; |
| | 17 | 812 | | int c_Columns = b._columns; |
| | 17 | 813 | | T[] A = a._matrix; |
| | 17 | 814 | | T[] B = b._matrix; |
| | | 815 | | T[] C; |
| | 17 | 816 | | if (c is not null && c._matrix.Length == c_Rows * c_Columns) |
| | 12 | 817 | | { |
| | 12 | 818 | | C = c._matrix; |
| | 12 | 819 | | c._rows = c_Rows; |
| | 12 | 820 | | c._columns = c_Columns; |
| | 12 | 821 | | } |
| | | 822 | | else |
| | 5 | 823 | | { |
| | 5 | 824 | | c = new Matrix<T>(c_Rows, c_Columns); |
| | 5 | 825 | | C = c._matrix; |
| | 5 | 826 | | } |
| | 104 | 827 | | for (int i = 0; i < c_Rows; i++) |
| | 35 | 828 | | { |
| | 35 | 829 | | int i_times_a_Columns = i * a_Columns; |
| | 35 | 830 | | int i_times_c_Columns = i * c_Columns; |
| | 234 | 831 | | for (int j = 0; j < c_Columns; j++) |
| | 82 | 832 | | { |
| | 82 | 833 | | T sum = Constant<T>.Zero; |
| | 632 | 834 | | for (int k = 0; k < a_Columns; k++) |
| | 234 | 835 | | { |
| | 234 | 836 | | sum = MultiplyAddImplementation<T>.Function(A[i_times_a_Columns + k], B[k * c_Columns + j], sum); |
| | 234 | 837 | | } |
| | 82 | 838 | | C[i_times_c_Columns + j] = sum; |
| | 82 | 839 | | } |
| | 35 | 840 | | } |
| | 17 | 841 | | } |
| | | 842 | | |
| | | 843 | | /// <summary>Does a standard (triple for looped) multiplication between matrices.</summary> |
| | | 844 | | /// <param name="a">The left matrix of the multiplication.</param> |
| | | 845 | | /// <param name="b">The right matrix of the multiplication.</param> |
| | | 846 | | /// <returns>The resulting matrix of the multiplication.</returns> |
| | | 847 | | public static Matrix<T> Multiply(Matrix<T> a, Matrix<T> b) |
| | 6 | 848 | | { |
| | 6 | 849 | | Matrix<T>? c = null; |
| | 6 | 850 | | Multiply(a, b, ref c); |
| | 5 | 851 | | return c!; |
| | 5 | 852 | | } |
| | | 853 | | |
| | | 854 | | /// <summary>Does a standard (triple for looped) multiplication between matrices.</summary> |
| | | 855 | | /// <param name="a">The left matrix of the multiplication.</param> |
| | | 856 | | /// <param name="b">The right matrix of the multiplication.</param> |
| | | 857 | | /// <returns>The resulting matrix of the multiplication.</returns> |
| | | 858 | | public static Matrix<T> operator *(Matrix<T> a, Matrix<T> b) |
| | 6 | 859 | | { |
| | 6 | 860 | | return Multiply(a, b); |
| | 5 | 861 | | } |
| | | 862 | | |
| | | 863 | | /// <summary>Does a standard (triple for looped) multiplication between matrices.</summary> |
| | | 864 | | /// <param name="b">The right matrix of the multiplication.</param> |
| | | 865 | | /// <param name="c">The resulting matrix of the multiplication.</param> |
| | | 866 | | public void Multiply(Matrix<T> b, ref Matrix<T>? c) |
| | 0 | 867 | | { |
| | 0 | 868 | | Multiply(this, b, ref c); |
| | 0 | 869 | | } |
| | | 870 | | |
| | | 871 | | /// <summary>Does a standard (triple for looped) multiplication between matrices.</summary> |
| | | 872 | | /// <param name="b">The right matrix of the multiplication.</param> |
| | | 873 | | /// <returns>The resulting matrix of the multiplication.</returns> |
| | | 874 | | public Matrix<T> Multiply(Matrix<T> b) |
| | 0 | 875 | | { |
| | 0 | 876 | | return this * b; |
| | 0 | 877 | | } |
| | | 878 | | |
| | | 879 | | #endregion |
| | | 880 | | |
| | | 881 | | #region Multiply (Matrix * Vector) |
| | | 882 | | |
| | | 883 | | /// <summary>Does a matrix-vector multiplication.</summary> |
| | | 884 | | /// <param name="a">The left matrix of the multiplication.</param> |
| | | 885 | | /// <param name="b">The right vector of the multiplication.</param> |
| | | 886 | | /// <param name="c">The resulting vector of the multiplication.</param> |
| | | 887 | | public static void Multiply(Matrix<T> a, Vector<T> b, ref Vector<T>? c) |
| | 5 | 888 | | { |
| | 5 | 889 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 5 | 890 | | if (b is null) throw new ArgumentNullException(nameof(b)); |
| | 5 | 891 | | int rows = a._rows; |
| | 5 | 892 | | int columns = a._columns; |
| | 6 | 893 | | if (sourceof(columns != b.Dimensions, out string c1)) throw new ArgumentException(c1); |
| | 4 | 894 | | T[] A = a._matrix; |
| | 4 | 895 | | T[] B = b._vector; |
| | | 896 | | T[] C; |
| | 4 | 897 | | if (c is not null && c.Dimensions == columns) |
| | 0 | 898 | | { |
| | 0 | 899 | | C = c._vector; |
| | 0 | 900 | | } |
| | | 901 | | else |
| | 4 | 902 | | { |
| | 4 | 903 | | c = new Vector<T>(columns); |
| | 4 | 904 | | C = c._vector; |
| | 4 | 905 | | } |
| | 32 | 906 | | for (int i = 0; i < rows; i++) |
| | 12 | 907 | | { |
| | 12 | 908 | | int i_times_columns = i * columns; |
| | 12 | 909 | | T sum = Constant<T>.Zero; |
| | 96 | 910 | | for (int j = 0; j < columns; j++) |
| | 36 | 911 | | { |
| | 36 | 912 | | sum = Addition(sum, Multiplication(A[i_times_columns + j], B[j])); |
| | 36 | 913 | | } |
| | 12 | 914 | | C[i] = sum; |
| | 12 | 915 | | } |
| | 4 | 916 | | } |
| | | 917 | | |
| | | 918 | | /// <summary>Does a matrix-vector multiplication.</summary> |
| | | 919 | | /// <param name="a">The left matrix of the multiplication.</param> |
| | | 920 | | /// <param name="b">The right vector of the multiplication.</param> |
| | | 921 | | /// <returns>The resulting vector of the multiplication.</returns> |
| | | 922 | | public static Vector<T> Multiply(Matrix<T> a, Vector<T> b) |
| | 5 | 923 | | { |
| | 5 | 924 | | Vector<T>? c = null; |
| | 5 | 925 | | Multiply(a, b, ref c); |
| | 4 | 926 | | return c!; |
| | 4 | 927 | | } |
| | | 928 | | |
| | | 929 | | /// <summary>Does a matrix-vector multiplication.</summary> |
| | | 930 | | /// <param name="a">The left matrix of the multiplication.</param> |
| | | 931 | | /// <param name="b">The right vector of the multiplication.</param> |
| | | 932 | | /// <returns>The resulting vector of the multiplication.</returns> |
| | | 933 | | public static Vector<T> operator *(Matrix<T> a, Vector<T> b) |
| | 5 | 934 | | { |
| | 5 | 935 | | return Multiply(a, b); |
| | 4 | 936 | | } |
| | | 937 | | |
| | | 938 | | /// <summary>Does a matrix-vector multiplication.</summary> |
| | | 939 | | /// <param name="b">The right vector of the multiplication.</param> |
| | | 940 | | /// <param name="c">The resulting vector of the multiplication.</param> |
| | | 941 | | public void Multiply(Vector<T> b, ref Vector<T>? c) |
| | 0 | 942 | | { |
| | 0 | 943 | | Multiply(this, b, ref c); |
| | 0 | 944 | | } |
| | | 945 | | |
| | | 946 | | /// <summary>Does a matrix-vector multiplication.</summary> |
| | | 947 | | /// <param name="b">The right vector of the multiplication.</param> |
| | | 948 | | /// <returns>The resulting vector of the multiplication.</returns> |
| | | 949 | | public Vector<T> Multiply(Vector<T> b) |
| | 0 | 950 | | { |
| | 0 | 951 | | return this * b; |
| | 0 | 952 | | } |
| | | 953 | | |
| | | 954 | | #endregion |
| | | 955 | | |
| | | 956 | | #region Multiply (Matrix * Scalar) |
| | | 957 | | |
| | | 958 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 959 | | /// <param name="a">The matrix to have the values multiplied.</param> |
| | | 960 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 961 | | /// <param name="c">The resulting matrix after the multiplications.</param> |
| | | 962 | | public static void Multiply(Matrix<T> a, T b, ref Matrix<T>? c) |
| | 4 | 963 | | { |
| | 4 | 964 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 4 | 965 | | T[] A = a._matrix; |
| | 4 | 966 | | int Length = A.Length; |
| | | 967 | | T[] C; |
| | 4 | 968 | | if (c is not null && c._matrix.Length == Length) |
| | 0 | 969 | | { |
| | 0 | 970 | | C = c._matrix; |
| | 0 | 971 | | c._rows = a._rows; |
| | 0 | 972 | | c._columns = a._columns; |
| | 0 | 973 | | } |
| | | 974 | | else |
| | 4 | 975 | | { |
| | 4 | 976 | | c = new Matrix<T>(a._rows, a._columns, Length); |
| | 4 | 977 | | C = c._matrix; |
| | 4 | 978 | | } |
| | 80 | 979 | | for (int i = 0; i < Length; i++) |
| | 36 | 980 | | { |
| | 36 | 981 | | C[i] = Multiplication(A[i], b); |
| | 36 | 982 | | } |
| | 4 | 983 | | } |
| | | 984 | | |
| | | 985 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 986 | | /// <param name="a">The matrix to have the values multiplied.</param> |
| | | 987 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 988 | | /// <returns>The resulting matrix after the multiplications.</returns> |
| | | 989 | | public static Matrix<T> Multiply(Matrix<T> a, T b) |
| | 4 | 990 | | { |
| | 4 | 991 | | Matrix<T>? c = null; |
| | 4 | 992 | | Multiply(a, b, ref c); |
| | 4 | 993 | | return c!; |
| | 4 | 994 | | } |
| | | 995 | | |
| | | 996 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 997 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 998 | | /// <param name="a">The matrix to have the values multiplied.</param> |
| | | 999 | | /// <returns>The resulting matrix after the multiplications.</returns> |
| | | 1000 | | public static Matrix<T> Multiply(T b, Matrix<T> a) |
| | 0 | 1001 | | { |
| | 0 | 1002 | | return Multiply(a, b); |
| | 0 | 1003 | | } |
| | | 1004 | | |
| | | 1005 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 1006 | | /// <param name="a">The matrix to have the values multiplied.</param> |
| | | 1007 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 1008 | | /// <returns>The resulting matrix after the multiplications.</returns> |
| | | 1009 | | public static Matrix<T> operator *(Matrix<T> a, T b) |
| | 4 | 1010 | | { |
| | 4 | 1011 | | return Multiply(a, b); |
| | 4 | 1012 | | } |
| | | 1013 | | |
| | | 1014 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 1015 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 1016 | | /// <param name="a">The matrix to have the values multiplied.</param> |
| | | 1017 | | /// <returns>The resulting matrix after the multiplications.</returns> |
| | | 1018 | | public static Matrix<T> operator *(T b, Matrix<T> a) |
| | 0 | 1019 | | { |
| | 0 | 1020 | | return Multiply(b, a); |
| | 0 | 1021 | | } |
| | | 1022 | | |
| | | 1023 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 1024 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 1025 | | /// <param name="c">The resulting matrix after the multiplications.</param> |
| | | 1026 | | public void Multiply(T b, ref Matrix<T>? c) |
| | 0 | 1027 | | { |
| | 0 | 1028 | | Multiply(this, b, ref c); |
| | 0 | 1029 | | } |
| | | 1030 | | |
| | | 1031 | | /// <summary>Multiplies all the values in a matrix by a scalar.</summary> |
| | | 1032 | | /// <param name="b">The scalar to multiply the values by.</param> |
| | | 1033 | | /// <returns>The resulting matrix after the multiplications.</returns> |
| | | 1034 | | public Matrix<T> Multiply(T b) |
| | 0 | 1035 | | { |
| | 0 | 1036 | | return this * b; |
| | 0 | 1037 | | } |
| | | 1038 | | |
| | | 1039 | | #endregion |
| | | 1040 | | |
| | | 1041 | | #region Divide (Matrix / Scalar) |
| | | 1042 | | |
| | | 1043 | | /// <summary>Divides all the values in the matrix by a scalar.</summary> |
| | | 1044 | | /// <param name="a">The matrix to divide the values of.</param> |
| | | 1045 | | /// <param name="b">The scalar to divide all the matrix values by.</param> |
| | | 1046 | | /// <param name="c">The resulting matrix after the division.</param> |
| | | 1047 | | public static void Divide(Matrix<T> a, T b, ref Matrix<T>? c) |
| | 4 | 1048 | | { |
| | 4 | 1049 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 4 | 1050 | | T[] A = a._matrix; |
| | 4 | 1051 | | int Length = A.Length; |
| | | 1052 | | T[] C; |
| | 4 | 1053 | | if (c is not null && c._matrix.Length == Length) |
| | 0 | 1054 | | { |
| | 0 | 1055 | | C = c._matrix; |
| | 0 | 1056 | | c._rows = a._rows; |
| | 0 | 1057 | | c._columns = a._columns; |
| | 0 | 1058 | | } |
| | | 1059 | | else |
| | 4 | 1060 | | { |
| | 4 | 1061 | | c = new Matrix<T>(a._rows, a._columns, Length); |
| | 4 | 1062 | | C = c._matrix; |
| | 4 | 1063 | | } |
| | 80 | 1064 | | for (int i = 0; i < Length; i++) |
| | 36 | 1065 | | { |
| | 36 | 1066 | | C[i] = Division(A[i], b); |
| | 36 | 1067 | | } |
| | 4 | 1068 | | } |
| | | 1069 | | |
| | | 1070 | | /// <summary>Divides all the values in the matrix by a scalar.</summary> |
| | | 1071 | | /// <param name="a">The matrix to divide the values of.</param> |
| | | 1072 | | /// <param name="b">The scalar to divide all the matrix values by.</param> |
| | | 1073 | | /// <returns>The resulting matrix after the division.</returns> |
| | | 1074 | | public static Matrix<T> Divide(Matrix<T> a, T b) |
| | 4 | 1075 | | { |
| | 4 | 1076 | | Matrix<T>? c = null; |
| | 4 | 1077 | | Divide(a, b, ref c); |
| | 4 | 1078 | | return c!; |
| | 4 | 1079 | | } |
| | | 1080 | | |
| | | 1081 | | /// <summary>Divides all the values in the matrix by a scalar.</summary> |
| | | 1082 | | /// <param name="a">The matrix to divide the values of.</param> |
| | | 1083 | | /// <param name="b">The scalar to divide all the matrix values by.</param> |
| | | 1084 | | /// <returns>The resulting matrix after the division.</returns> |
| | | 1085 | | public static Matrix<T> operator /(Matrix<T> a, T b) |
| | 4 | 1086 | | { |
| | 4 | 1087 | | return Divide(a, b); |
| | 4 | 1088 | | } |
| | | 1089 | | |
| | | 1090 | | /// <summary>Divides all the values in the matrix by a scalar.</summary> |
| | | 1091 | | /// <param name="b">The scalar to divide all the matrix values by.</param> |
| | | 1092 | | /// <param name="c">The resulting matrix after the division.</param> |
| | | 1093 | | public void Divide(T b, ref Matrix<T>? c) |
| | 0 | 1094 | | { |
| | 0 | 1095 | | Divide(this, b, ref c); |
| | 0 | 1096 | | } |
| | | 1097 | | |
| | | 1098 | | /// <summary>Divides all the values in the matrix by a scalar.</summary> |
| | | 1099 | | /// <param name="b">The scalar to divide all the matrix values by.</param> |
| | | 1100 | | /// <returns>The resulting matrix after the division.</returns> |
| | | 1101 | | public Matrix<T> Divide(T b) |
| | 0 | 1102 | | { |
| | 0 | 1103 | | return this / b; |
| | 0 | 1104 | | } |
| | | 1105 | | |
| | | 1106 | | #endregion |
| | | 1107 | | |
| | | 1108 | | #region Power (Matrix ^ Scalar) |
| | | 1109 | | |
| | | 1110 | | /// <summary>Applies a power to a square matrix.</summary> |
| | | 1111 | | /// <param name="a">The matrix to be powered by.</param> |
| | | 1112 | | /// <param name="b">The power to apply to the matrix.</param> |
| | | 1113 | | /// <param name="c">The resulting matrix of the power operation.</param> |
| | | 1114 | | public static void Power(Matrix<T> a, int b, ref Matrix<T>? c) |
| | 5 | 1115 | | { |
| | 5 | 1116 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 6 | 1117 | | if (sourceof(!a.IsSquare, out string c1)) throw new ArgumentException(c1); |
| | 4 | 1118 | | if (sourceof(b < 0, out string c2)) throw new ArgumentOutOfRangeException(nameof(b), b, c2); |
| | 4 | 1119 | | if (b is 0) |
| | 0 | 1120 | | { |
| | 0 | 1121 | | if (c is not null && c._matrix.Length == a._matrix.Length) |
| | 0 | 1122 | | { |
| | 0 | 1123 | | c._rows = a._rows; |
| | 0 | 1124 | | c._columns = a._columns; |
| | 0 | 1125 | | Format(c, (x, y) => x == y ? Constant<T>.One : Constant<T>.Zero); |
| | 0 | 1126 | | } |
| | | 1127 | | else |
| | 0 | 1128 | | { |
| | 0 | 1129 | | c = Matrix<T>.FactoryIdentity(a._rows, a._columns); |
| | 0 | 1130 | | } |
| | 0 | 1131 | | return; |
| | | 1132 | | } |
| | 4 | 1133 | | if (c is not null && c._matrix.Length == a._matrix.Length) |
| | 0 | 1134 | | { |
| | 0 | 1135 | | c._rows = a._rows; |
| | 0 | 1136 | | c._columns = a._columns; |
| | 0 | 1137 | | T[] A = a._matrix; |
| | 0 | 1138 | | T[] C = c._matrix; |
| | 0 | 1139 | | for (int i = 0; i < a._matrix.Length; i++) |
| | 0 | 1140 | | { |
| | 0 | 1141 | | C[i] = A[i]; |
| | 0 | 1142 | | } |
| | 0 | 1143 | | } |
| | | 1144 | | else |
| | 4 | 1145 | | { |
| | 4 | 1146 | | c = a.Clone(); |
| | 4 | 1147 | | } |
| | 4 | 1148 | | Matrix<T>? d = new(a._rows, a._columns, a._matrix.Length); |
| | 32 | 1149 | | for (int i = 0; i < b; i++) |
| | 12 | 1150 | | { |
| | 12 | 1151 | | Multiply(c, a, ref d); |
| | 12 | 1152 | | Matrix<T>? temp = d; |
| | 12 | 1153 | | d = c; |
| | 12 | 1154 | | c = d; |
| | 12 | 1155 | | } |
| | 4 | 1156 | | } |
| | | 1157 | | |
| | | 1158 | | /// <summary>Applies a power to a square matrix.</summary> |
| | | 1159 | | /// <param name="a">The matrix to be powered by.</param> |
| | | 1160 | | /// <param name="b">The power to apply to the matrix.</param> |
| | | 1161 | | /// <returns>The resulting matrix of the power operation.</returns> |
| | | 1162 | | public static Matrix<T> Power(Matrix<T> a, int b) |
| | 5 | 1163 | | { |
| | 5 | 1164 | | Matrix<T>? c = null; |
| | 5 | 1165 | | Power(a, b, ref c); |
| | 4 | 1166 | | return c!; |
| | 4 | 1167 | | } |
| | | 1168 | | |
| | | 1169 | | /// <summary>Applies a power to a square matrix.</summary> |
| | | 1170 | | /// <param name="a">The matrix to be powered by.</param> |
| | | 1171 | | /// <param name="b">The power to apply to the matrix.</param> |
| | | 1172 | | /// <returns>The resulting matrix of the power operation.</returns> |
| | | 1173 | | [Obsolete("Please use the Power method. The ^ operator in C# does not follow the mathematical order of operations.", t |
| | | 1174 | | public static Matrix<T> operator ^(Matrix<T> a, int b) |
| | 0 | 1175 | | { |
| | 0 | 1176 | | return Power(a, b); |
| | 0 | 1177 | | } |
| | | 1178 | | |
| | | 1179 | | /// <summary>Applies a power to a square matrix.</summary> |
| | | 1180 | | /// <param name="b">The power to apply to the matrix.</param> |
| | | 1181 | | /// <param name="c">The resulting matrix of the power operation.</param> |
| | | 1182 | | public void Power(int b, ref Matrix<T>? c) |
| | 0 | 1183 | | { |
| | 0 | 1184 | | Power(this, b, ref c); |
| | 0 | 1185 | | } |
| | | 1186 | | |
| | | 1187 | | /// <summary>Applies a power to a square matrix.</summary> |
| | | 1188 | | /// <param name="b">The power to apply to the matrix.</param> |
| | | 1189 | | /// <returns>The resulting matrix of the power operation.</returns> |
| | | 1190 | | public Matrix<T> Power(int b) |
| | 5 | 1191 | | { |
| | 5 | 1192 | | return Power(this, b); |
| | 4 | 1193 | | } |
| | | 1194 | | |
| | | 1195 | | #endregion |
| | | 1196 | | |
| | | 1197 | | #region Determinant |
| | | 1198 | | |
| | | 1199 | | /// <summary>Computes the determinant of a square matrix.</summary> |
| | | 1200 | | /// <param name="a">The matrix to compute the determinant of.</param> |
| | | 1201 | | /// <returns>The computed determinant.</returns> |
| | | 1202 | | public static T Determinant(Matrix<T> a) |
| | 2 | 1203 | | { |
| | 2 | 1204 | | return DeterminantGaussian(a); |
| | | 1205 | | |
| | | 1206 | | #region Old Version |
| | | 1207 | | |
| | | 1208 | | // if (a is null) |
| | | 1209 | | // { |
| | | 1210 | | // throw new ArgumentNullException(nameof(a)); |
| | | 1211 | | // } |
| | | 1212 | | // if (!a.IsSquare) |
| | | 1213 | | // { |
| | | 1214 | | // throw new MathematicsException("Argument invalid !(" + nameof(a) + "." + nameof(a.IsSquare) + ")"); |
| | | 1215 | | // } |
| | | 1216 | | // T determinant = Constant<T>.One; |
| | | 1217 | | // Matrix<T> rref = a.Clone(); |
| | | 1218 | | // int a_rows = a._rows; |
| | | 1219 | | // for (int i = 0; i < a_rows; i++) |
| | | 1220 | | // { |
| | | 1221 | | // if (Compute.Equal(rref[i, i], Constant<T>.Zero)) |
| | | 1222 | | // { |
| | | 1223 | | // for (int j = i + 1; j < rref.Rows; j++) |
| | | 1224 | | // { |
| | | 1225 | | // if (Compute.NotEqual(rref.Get(j, i), Constant<T>.Zero)) |
| | | 1226 | | // { |
| | | 1227 | | // SwapRows(rref, i, j); |
| | | 1228 | | // determinant = Compute.Multiply(determinant, Constant<T>.NegativeOne); |
| | | 1229 | | // } |
| | | 1230 | | // } |
| | | 1231 | | // } |
| | | 1232 | | // determinant = Compute.Multiply(determinant, rref.Get(i, i)); |
| | | 1233 | | // T temp_rowMultiplication = Compute.Divide(Constant<T>.One, rref.Get(i, i)); |
| | | 1234 | | // RowMultiplication(rref, i, temp_rowMultiplication); |
| | | 1235 | | // for (int j = i + 1; j < rref.Rows; j++) |
| | | 1236 | | // { |
| | | 1237 | | // T scalar = Compute.Negate(rref.Get(j, i)); |
| | | 1238 | | // RowAddition(rref, j, i, scalar); |
| | | 1239 | | // } |
| | | 1240 | | // for (int j = i - 1; j >= 0; j--) |
| | | 1241 | | // { |
| | | 1242 | | // T scalar = Compute.Negate(rref.Get(j, i)); |
| | | 1243 | | // RowAddition(rref, j, i, scalar); |
| | | 1244 | | // } |
| | | 1245 | | // } |
| | | 1246 | | // return determinant; |
| | | 1247 | | |
| | | 1248 | | #endregion |
| | | 1249 | | |
| | 2 | 1250 | | } |
| | | 1251 | | |
| | | 1252 | | /// <summary> |
| | | 1253 | | /// Computes the determinant of a square matrix via Gaussian elimination.<br/> |
| | | 1254 | | /// Runtime: O((n^3 + 2n^−3) / 3) |
| | | 1255 | | /// </summary> |
| | | 1256 | | /// <param name="a">The matrix to compute the determinant of.</param> |
| | | 1257 | | /// <returns>The computed determinant.</returns> |
| | | 1258 | | public static T DeterminantGaussian(Matrix<T> a) |
| | 11 | 1259 | | { |
| | 11 | 1260 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 12 | 1261 | | if (sourceof(!a.IsSquare, out string c1)) throw new ArgumentException(c1); |
| | 10 | 1262 | | return GetDeterminantGaussian(a, a.Rows); |
| | 10 | 1263 | | } |
| | | 1264 | | |
| | | 1265 | | /// <summary> |
| | | 1266 | | /// Computes the determinant of a square matrix via Laplace's method.<br/> |
| | | 1267 | | /// Runtime: O(n(2^(n − 1) − 1)) |
| | | 1268 | | /// </summary> |
| | | 1269 | | /// <param name="a">The matrix to compute the determinant of.</param> |
| | | 1270 | | /// <returns>The computed determinant.</returns> |
| | | 1271 | | public static T DeterminantLaplace(Matrix<T> a) |
| | 9 | 1272 | | { |
| | 9 | 1273 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 10 | 1274 | | if (sourceof(!a.IsSquare, out string c1)) throw new ArgumentException(c1); |
| | 8 | 1275 | | return GetDeterminantLaplace(a, a.Rows); |
| | 8 | 1276 | | } |
| | | 1277 | | |
| | | 1278 | | /// <summary>Computes the determinant of a square matrix.</summary> |
| | | 1279 | | /// <returns>The computed determinant.</returns> |
| | | 1280 | | public T Determinant() |
| | 0 | 1281 | | { |
| | 0 | 1282 | | return DeterminantGaussian(this); |
| | 0 | 1283 | | } |
| | | 1284 | | |
| | | 1285 | | /// <summary>Computes the determinant of a square matrix via Laplace's method.</summary> |
| | | 1286 | | /// <returns>The computed determinant.</returns> |
| | | 1287 | | public T DeterminantLaplace() |
| | 9 | 1288 | | { |
| | 9 | 1289 | | return DeterminantLaplace(this); |
| | 8 | 1290 | | } |
| | | 1291 | | |
| | | 1292 | | /// <summary>Computes the determinant of a square matrix via Gaussian elimination.</summary> |
| | | 1293 | | /// <returns>The computed determinant.</returns> |
| | | 1294 | | public T DeterminantGaussian() |
| | 9 | 1295 | | { |
| | 9 | 1296 | | return DeterminantGaussian(this); |
| | 8 | 1297 | | } |
| | | 1298 | | |
| | | 1299 | | #endregion |
| | | 1300 | | |
| | | 1301 | | #region Trace |
| | | 1302 | | |
| | | 1303 | | /// <summary>Computes the trace of a square matrix.</summary> |
| | | 1304 | | /// <param name="a">The matrix to compute the trace of.</param> |
| | | 1305 | | /// <returns>The computed trace.</returns> |
| | | 1306 | | public static T Trace(Matrix<T> a) |
| | 5 | 1307 | | { |
| | 5 | 1308 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 6 | 1309 | | if (sourceof(!a.IsSquare, out string c1)) throw new ArgumentException(c1); |
| | 4 | 1310 | | T trace = Addition((Action<T> step) => |
| | 4 | 1311 | | { |
| | 4 | 1312 | | T[] A = a._matrix; |
| | 4 | 1313 | | int rows = a.Rows; |
| | 32 | 1314 | | for (int i = 0; i < rows; i++) |
| | 12 | 1315 | | { |
| | 12 | 1316 | | step(A[i * rows + i]); |
| | 12 | 1317 | | } |
| | 8 | 1318 | | }); |
| | 4 | 1319 | | return trace; |
| | 4 | 1320 | | } |
| | | 1321 | | |
| | | 1322 | | /// <summary>Computes the trace of a square matrix.</summary> |
| | | 1323 | | /// <returns>The computed trace.</returns> |
| | | 1324 | | public T Trace() |
| | 5 | 1325 | | { |
| | 5 | 1326 | | return Trace(this); |
| | 4 | 1327 | | } |
| | | 1328 | | |
| | | 1329 | | #endregion |
| | | 1330 | | |
| | | 1331 | | #region Minor |
| | | 1332 | | |
| | | 1333 | | #pragma warning disable CS1572, SA1625, SA1617 |
| | | 1334 | | |
| | | 1335 | | /// <summary>Gets the minor of a matrix.</summary> |
| | | 1336 | | /// <param name="a">The matrix to get the minor of.</param> |
| | | 1337 | | /// <param name="row">The restricted row to form the minor.</param> |
| | | 1338 | | /// <param name="column">The restricted column to form the minor.</param> |
| | | 1339 | | /// <param name="b">The minor of the matrix.</param> |
| | | 1340 | | /// <returns>The minor of the matrix.</returns> |
| | | 1341 | | [Obsolete(NotIntended, true)] |
| | 1 | 1342 | | public static void XML_Minor() => throw new DocumentationMethodException(); |
| | | 1343 | | |
| | | 1344 | | #pragma warning restore CS1572, SA1625, SA1617 |
| | | 1345 | | |
| | | 1346 | | /// <inheritdoc cref="XML_Minor"/> |
| | | 1347 | | public Matrix<T> Minor(int row, int column) => |
| | 26 | 1348 | | Minor(this, row, column); |
| | | 1349 | | |
| | | 1350 | | /// <inheritdoc cref="XML_Minor"/> |
| | | 1351 | | public static Matrix<T> Minor(Matrix<T> a, int row, int column) |
| | 27 | 1352 | | { |
| | 27 | 1353 | | Matrix<T>? b = null; |
| | 27 | 1354 | | Minor(a, row, column, ref b); |
| | 13 | 1355 | | return b!; |
| | 13 | 1356 | | } |
| | | 1357 | | |
| | | 1358 | | /// <inheritdoc cref="XML_Minor"/> |
| | | 1359 | | public void Minor(int row, int column, ref Matrix<T>? b) => |
| | 0 | 1360 | | Minor(this, row, column, ref b); |
| | | 1361 | | |
| | | 1362 | | /// <inheritdoc cref="XML_Minor"/> |
| | | 1363 | | public static void Minor(Matrix<T> a, int row, int column, ref Matrix<T>? b) |
| | 27 | 1364 | | { |
| | 28 | 1365 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 27 | 1366 | | if (sourceof(a._rows < 2 || a._columns < 2, out string c1)) throw new ArgumentException(c1); |
| | 33 | 1367 | | if (sourceof(!(0 <= row && row < a._rows), out string c2)) throw new ArgumentOutOfRangeException(nameof(row), row, c |
| | 21 | 1368 | | if (sourceof(!(0 <= column && column < a._columns), out string c3)) throw new ArgumentOutOfRangeException(nameof(col |
| | 13 | 1369 | | if (ReferenceEquals(a, b)) |
| | 0 | 1370 | | { |
| | 0 | 1371 | | a = a.Clone(); |
| | 0 | 1372 | | } |
| | 13 | 1373 | | int a_rows = a._rows; |
| | 13 | 1374 | | int a_columns = a._columns; |
| | 13 | 1375 | | int b_rows = a_rows - 1; |
| | 13 | 1376 | | int b_columns = a_columns - 1; |
| | 13 | 1377 | | int b_length = b_rows * b_columns; |
| | 13 | 1378 | | if (b is null || b._matrix.Length != b_length) |
| | 13 | 1379 | | { |
| | 13 | 1380 | | b = new Matrix<T>(b_rows, b_columns, b_length); |
| | 13 | 1381 | | } |
| | | 1382 | | else |
| | 0 | 1383 | | { |
| | 0 | 1384 | | b._rows = b_rows; |
| | 0 | 1385 | | b._columns = b_columns; |
| | 0 | 1386 | | } |
| | 13 | 1387 | | T[] B = b._matrix; |
| | 13 | 1388 | | T[] A = a._matrix; |
| | 109 | 1389 | | for (int i = 0, m = 0; i < a_rows; i++) |
| | 35 | 1390 | | { |
| | 35 | 1391 | | if (i == row) |
| | 13 | 1392 | | { |
| | 13 | 1393 | | continue; |
| | | 1394 | | } |
| | 22 | 1395 | | int i_times_a_columns = i * a_columns; |
| | 22 | 1396 | | int m_times_b_columns = m * b_columns; |
| | 190 | 1397 | | for (int j = 0, n = 0; j < a_columns; j++) |
| | 62 | 1398 | | { |
| | 62 | 1399 | | if (j == column) |
| | 22 | 1400 | | { |
| | 22 | 1401 | | continue; |
| | | 1402 | | } |
| | 40 | 1403 | | T temp = A[i_times_a_columns + j]; |
| | 40 | 1404 | | B[m_times_b_columns + n] = temp; |
| | 40 | 1405 | | n++; |
| | 40 | 1406 | | } |
| | 22 | 1407 | | m++; |
| | 22 | 1408 | | } |
| | 13 | 1409 | | } |
| | | 1410 | | |
| | | 1411 | | #endregion |
| | | 1412 | | |
| | | 1413 | | #region ConcatenateRowWise |
| | | 1414 | | |
| | | 1415 | | /// <summary>Combines two matrices from left to right |
| | | 1416 | | /// (result.Rows = left.Rows AND result.Columns = left.Columns + right.Columns).</summary> |
| | | 1417 | | /// <param name="a">The left matrix of the concatenation.</param> |
| | | 1418 | | /// <param name="b">The right matrix of the concatenation.</param> |
| | | 1419 | | /// <param name="c">The resulting matrix of the concatenation.</param> |
| | | 1420 | | public static void ConcatenateRowWise(Matrix<T> a, Matrix<T> b, ref Matrix<T>? c) |
| | 3 | 1421 | | { |
| | 3 | 1422 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 3 | 1423 | | if (b is null) throw new ArgumentNullException(nameof(b)); |
| | 3 | 1424 | | if (sourceof(a.Rows != b.Rows, out string c1)) throw new ArgumentException(c1); |
| | 3 | 1425 | | int c_rows = a._rows; |
| | 3 | 1426 | | int c_columns = a._columns + b._columns; |
| | 3 | 1427 | | int c_length = c_rows * c_columns; |
| | 3 | 1428 | | if (c is null || |
| | 3 | 1429 | | c._matrix.Length != c_length || |
| | 3 | 1430 | | object.ReferenceEquals(a, c) || |
| | 3 | 1431 | | object.ReferenceEquals(b, c)) |
| | 3 | 1432 | | { |
| | 3 | 1433 | | c = new Matrix<T>(c_rows, c_columns, c_length); |
| | 3 | 1434 | | } |
| | | 1435 | | else |
| | 0 | 1436 | | { |
| | 0 | 1437 | | c._rows = c_rows; |
| | 0 | 1438 | | c._columns = c_columns; |
| | 0 | 1439 | | } |
| | 3 | 1440 | | T[] A = a._matrix; |
| | 3 | 1441 | | T[] B = b._matrix; |
| | 3 | 1442 | | T[] C = c._matrix; |
| | 20 | 1443 | | for (int i = 0; i < c_rows; i++) |
| | 7 | 1444 | | { |
| | 7 | 1445 | | int i_times_a_columns = i * a.Columns; |
| | 7 | 1446 | | int i_times_b_columns = i * b.Columns; |
| | 7 | 1447 | | int i_times_c_columns = i * c_columns; |
| | 88 | 1448 | | for (int j = 0; j < c_columns; j++) |
| | 37 | 1449 | | { |
| | 37 | 1450 | | if (j < a.Columns) |
| | 19 | 1451 | | { |
| | 19 | 1452 | | C[i_times_c_columns + j] = A[i_times_a_columns + j]; |
| | 19 | 1453 | | } |
| | | 1454 | | else |
| | 18 | 1455 | | { |
| | 18 | 1456 | | C[i_times_c_columns + j] = B[i_times_b_columns + j - a.Columns]; |
| | 18 | 1457 | | } |
| | 37 | 1458 | | } |
| | 7 | 1459 | | } |
| | 3 | 1460 | | } |
| | | 1461 | | |
| | | 1462 | | /// <summary>Combines two matrices from left to right |
| | | 1463 | | /// (result.Rows = left.Rows AND result.Columns = left.Columns + right.Columns).</summary> |
| | | 1464 | | /// <param name="a">The left matrix of the concatenation.</param> |
| | | 1465 | | /// <param name="b">The right matrix of the concatenation.</param> |
| | | 1466 | | /// <returns>The resulting matrix of the concatenation.</returns> |
| | | 1467 | | public static Matrix<T> ConcatenateRowWise(Matrix<T> a, Matrix<T> b) |
| | 3 | 1468 | | { |
| | 3 | 1469 | | Matrix<T>? c = null; |
| | 3 | 1470 | | ConcatenateRowWise(a, b, ref c); |
| | 3 | 1471 | | return c!; |
| | 3 | 1472 | | } |
| | | 1473 | | |
| | | 1474 | | /// <summary>Combines two matrices from left to right |
| | | 1475 | | /// (result.Rows = left.Rows AND result.Columns = left.Columns + right.Columns).</summary> |
| | | 1476 | | /// <param name="b">The right matrix of the concatenation.</param> |
| | | 1477 | | /// <param name="c">The resulting matrix of the concatenation.</param> |
| | | 1478 | | public void ConcatenateRowWise(Matrix<T> b, ref Matrix<T>? c) |
| | 0 | 1479 | | { |
| | 0 | 1480 | | ConcatenateRowWise(this, b, ref c); |
| | 0 | 1481 | | } |
| | | 1482 | | |
| | | 1483 | | /// <summary>Combines two matrices from left to right |
| | | 1484 | | /// (result.Rows = left.Rows AND result.Columns = left.Columns + right.Columns).</summary> |
| | | 1485 | | /// <param name="b">The right matrix of the concatenation.</param> |
| | | 1486 | | /// <returns>The resulting matrix of the concatenation.</returns> |
| | | 1487 | | public Matrix<T> ConcatenateRowWise(Matrix<T> b) |
| | 3 | 1488 | | { |
| | 3 | 1489 | | return ConcatenateRowWise(this, b); |
| | 3 | 1490 | | } |
| | | 1491 | | |
| | | 1492 | | #endregion |
| | | 1493 | | |
| | | 1494 | | #region Echelon |
| | | 1495 | | |
| | | 1496 | | #if false |
| | | 1497 | | |
| | | 1498 | | /// <summary>Calculates the echelon of a matrix (aka REF).</summary> |
| | | 1499 | | /// <param name="a">The matrix to calculate the echelon of (aka REF).</param> |
| | | 1500 | | /// <param name="b">The echelon of the matrix (aka REF).</param> |
| | | 1501 | | /// <bug>Failing for non-floating point rational types due to zero how values are being compared.</bug> |
| | | 1502 | | public static void Echelon(Matrix<T> a, ref Matrix<T> b) |
| | | 1503 | | { |
| | | 1504 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | | 1505 | | if (ReferenceEquals(a, b)) |
| | | 1506 | | { |
| | | 1507 | | a = a.Clone(); |
| | | 1508 | | } |
| | | 1509 | | int Rows = a.Rows; |
| | | 1510 | | if (b is not null && b._matrix.Length == a._matrix.Length) |
| | | 1511 | | { |
| | | 1512 | | b._rows = Rows; |
| | | 1513 | | b._columns = a._columns; |
| | | 1514 | | CloneContents(a, b); |
| | | 1515 | | } |
| | | 1516 | | else |
| | | 1517 | | { |
| | | 1518 | | b = a.Clone(); |
| | | 1519 | | } |
| | | 1520 | | for (int i = 0; i < Rows; i++) |
| | | 1521 | | { |
| | | 1522 | | if (Equality(b.Get(i, i), Constant<T>.Zero)) |
| | | 1523 | | { |
| | | 1524 | | for (int j = i + 1; j < Rows; j++) |
| | | 1525 | | { |
| | | 1526 | | if (Inequality(b.Get(j, i), Constant<T>.Zero)) |
| | | 1527 | | { |
| | | 1528 | | SwapRows(b, i, j); |
| | | 1529 | | } |
| | | 1530 | | } |
| | | 1531 | | } |
| | | 1532 | | if (Equality(b.Get(i, i), Constant<T>.Zero)) |
| | | 1533 | | { |
| | | 1534 | | continue; |
| | | 1535 | | } |
| | | 1536 | | if (Inequality(b.Get(i, i), Constant<T>.One)) |
| | | 1537 | | { |
| | | 1538 | | for (int j = i + 1; j < Rows; j++) |
| | | 1539 | | { |
| | | 1540 | | if (Equality(b.Get(j, i), Constant<T>.One)) |
| | | 1541 | | { |
| | | 1542 | | SwapRows(b, i, j); |
| | | 1543 | | } |
| | | 1544 | | } |
| | | 1545 | | } |
| | | 1546 | | T rowMultipier = Division(Constant<T>.One, b.Get(i, i)); |
| | | 1547 | | RowMultiplication(b, i, rowMultipier); |
| | | 1548 | | for (int j = i + 1; j < Rows; j++) |
| | | 1549 | | { |
| | | 1550 | | T rowAddend = Negation(b.Get(j, i)); |
| | | 1551 | | RowAddition(b, j, i, rowAddend); |
| | | 1552 | | } |
| | | 1553 | | } |
| | | 1554 | | } |
| | | 1555 | | |
| | | 1556 | | /// <summary>Calculates the echelon of a matrix (aka REF).</summary> |
| | | 1557 | | /// <param name="a">The matrix to calculate the echelon of (aka REF).</param> |
| | | 1558 | | /// <returns>The echelon of the matrix (aka REF).</returns> |
| | | 1559 | | public static Matrix<T> Echelon(Matrix<T> a) |
| | | 1560 | | { |
| | | 1561 | | Matrix<T> b = null; |
| | | 1562 | | Echelon(a, ref b); |
| | | 1563 | | return b; |
| | | 1564 | | } |
| | | 1565 | | |
| | | 1566 | | /// <summary>Calculates the echelon of a matrix (aka REF).</summary> |
| | | 1567 | | /// <param name="b">The echelon of the matrix (aka REF).</param> |
| | | 1568 | | public void Echelon(ref Matrix<T> b) |
| | | 1569 | | { |
| | | 1570 | | Echelon(this, ref b); |
| | | 1571 | | } |
| | | 1572 | | |
| | | 1573 | | /// <summary>Calculates the echelon of a matrix (aka REF).</summary> |
| | | 1574 | | /// <returns>The echelon of the matrix (aka REF).</returns> |
| | | 1575 | | public Matrix<T> Echelon() |
| | | 1576 | | { |
| | | 1577 | | return Echelon(this); |
| | | 1578 | | } |
| | | 1579 | | |
| | | 1580 | | #endif |
| | | 1581 | | |
| | | 1582 | | #endregion |
| | | 1583 | | |
| | | 1584 | | #region ReducedEchelon |
| | | 1585 | | |
| | | 1586 | | /// <summary>Calculates the echelon of a matrix and reduces it (aka RREF).</summary> |
| | | 1587 | | /// <param name="a">The matrix matrix to calculate the reduced echelon of (aka RREF).</param> |
| | | 1588 | | /// <param name="b">The reduced echelon of the matrix (aka RREF).</param> |
| | | 1589 | | public static void ReducedEchelon(Matrix<T> a, ref Matrix<T>? b) |
| | 4 | 1590 | | { |
| | 4 | 1591 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 4 | 1592 | | int Rows = a.Rows; |
| | 4 | 1593 | | int Columns = a.Columns; |
| | 4 | 1594 | | if (object.ReferenceEquals(a, b)) |
| | 0 | 1595 | | { |
| | 0 | 1596 | | b = a.Clone(); |
| | 0 | 1597 | | } |
| | 4 | 1598 | | else if (b is not null && b._matrix.Length == a._matrix.Length) |
| | 0 | 1599 | | { |
| | 0 | 1600 | | b._rows = Rows; |
| | 0 | 1601 | | b._columns = a._columns; |
| | 0 | 1602 | | CloneContents(a, b); |
| | 0 | 1603 | | } |
| | | 1604 | | else |
| | 4 | 1605 | | { |
| | 4 | 1606 | | b = a.Clone(); |
| | 4 | 1607 | | } |
| | 4 | 1608 | | int lead = 0; |
| | 32 | 1609 | | for (int r = 0; r < Rows; r++) |
| | 12 | 1610 | | { |
| | 12 | 1611 | | if (Columns <= lead) break; |
| | 12 | 1612 | | int i = r; |
| | 12 | 1613 | | while (Equate(b.Get(i, lead), Constant<T>.Zero)) |
| | 4 | 1614 | | { |
| | 4 | 1615 | | i++; |
| | 4 | 1616 | | if (i == Rows) |
| | 4 | 1617 | | { |
| | 4 | 1618 | | i = r; |
| | 4 | 1619 | | lead++; |
| | 4 | 1620 | | if (Columns == lead) |
| | 4 | 1621 | | { |
| | 4 | 1622 | | lead--; |
| | 4 | 1623 | | break; |
| | | 1624 | | } |
| | 0 | 1625 | | } |
| | 0 | 1626 | | } |
| | 96 | 1627 | | for (int j = 0; j < Columns; j++) |
| | 36 | 1628 | | { |
| | 36 | 1629 | | T temp = b.Get(r, j); |
| | 36 | 1630 | | b.Set(r, j, b.Get(i, j)); |
| | 36 | 1631 | | b.Set(i, j, temp); |
| | 36 | 1632 | | } |
| | 12 | 1633 | | T div = b.Get(r, lead); |
| | 12 | 1634 | | if (Inequate(div, Constant<T>.Zero)) |
| | 8 | 1635 | | { |
| | 64 | 1636 | | for (int j = 0; j < Columns; j++) |
| | 24 | 1637 | | { |
| | 24 | 1638 | | b.Set(r, j, Division(b.Get(r, j), div)); |
| | 24 | 1639 | | } |
| | 8 | 1640 | | } |
| | 96 | 1641 | | for (int j = 0; j < Rows; j++) |
| | 36 | 1642 | | { |
| | 36 | 1643 | | if (j != r) |
| | 24 | 1644 | | { |
| | 24 | 1645 | | T sub = b.Get(j, lead); |
| | 192 | 1646 | | for (int k = 0; k < Columns; k++) |
| | 72 | 1647 | | { |
| | 72 | 1648 | | b.Set(j, k, Subtraction(b.Get(j, k), Multiplication(sub, b.Get(r, k)))); |
| | 72 | 1649 | | } |
| | 24 | 1650 | | } |
| | 36 | 1651 | | } |
| | 12 | 1652 | | lead++; |
| | 12 | 1653 | | } |
| | | 1654 | | |
| | | 1655 | | #region Old Version |
| | | 1656 | | |
| | | 1657 | | // if (a is null) |
| | | 1658 | | // { |
| | | 1659 | | // throw new ArgumentNullException(nameof(a)); |
| | | 1660 | | // } |
| | | 1661 | | // if (object.ReferenceEquals(a, b)) |
| | | 1662 | | // { |
| | | 1663 | | // a = a.Clone(); |
| | | 1664 | | // } |
| | | 1665 | | // int Rows = a.Rows; |
| | | 1666 | | // if (b is not null && b._matrix.Length == a._matrix.Length) |
| | | 1667 | | // { |
| | | 1668 | | // b._rows = Rows; |
| | | 1669 | | // b._columns = a._columns; |
| | | 1670 | | // CloneContents(a, b); |
| | | 1671 | | // } |
| | | 1672 | | // else |
| | | 1673 | | // { |
| | | 1674 | | // b = a.Clone(); |
| | | 1675 | | // } |
| | | 1676 | | // for (int i = 0; i < Rows; i++) |
| | | 1677 | | // { |
| | | 1678 | | // if (Compute.Equal(b.Get(i, i), Constant<T>.Zero)) |
| | | 1679 | | // { |
| | | 1680 | | // for (int j = i + 1; j < Rows; j++) |
| | | 1681 | | // { |
| | | 1682 | | // if (Compute.NotEqual(b.Get(j, i), Constant<T>.Zero)) |
| | | 1683 | | // { |
| | | 1684 | | // SwapRows(b, i, j); |
| | | 1685 | | // } |
| | | 1686 | | // } |
| | | 1687 | | // } |
| | | 1688 | | // if (Compute.Equal(b.Get(i, i), Constant<T>.Zero)) |
| | | 1689 | | // { |
| | | 1690 | | // continue; |
| | | 1691 | | // } |
| | | 1692 | | // if (Compute.NotEqual(b.Get(i, i), Constant<T>.One)) |
| | | 1693 | | // { |
| | | 1694 | | // for (int j = i + 1; j < Rows; j++) |
| | | 1695 | | // { |
| | | 1696 | | // if (Compute.Equal(b.Get(j, i), Constant<T>.One)) |
| | | 1697 | | // { |
| | | 1698 | | // SwapRows(b, i, j); |
| | | 1699 | | // } |
| | | 1700 | | // } |
| | | 1701 | | // } |
| | | 1702 | | // T rowMiltiplier = Compute.Divide(Constant<T>.One, b.Get(i, i)); |
| | | 1703 | | // RowMultiplication(b, i, rowMiltiplier); |
| | | 1704 | | // for (int j = i + 1; j < Rows; j++) |
| | | 1705 | | // { |
| | | 1706 | | // T rowAddend = Compute.Negate(b.Get(j, i)); |
| | | 1707 | | // RowAddition(b, j, i, rowAddend); |
| | | 1708 | | // } |
| | | 1709 | | // for (int j = i - 1; j >= 0; j--) |
| | | 1710 | | // { |
| | | 1711 | | // T rowAddend = Compute.Negate(b.Get(j, i)); |
| | | 1712 | | // RowAddition(b, j, i, rowAddend); |
| | | 1713 | | // } |
| | | 1714 | | // } |
| | | 1715 | | |
| | | 1716 | | #endregion |
| | 4 | 1717 | | } |
| | | 1718 | | |
| | | 1719 | | /// <summary>Calculates the echelon of a matrix and reduces it (aka RREF).</summary> |
| | | 1720 | | /// <param name="a">The matrix matrix to calculate the reduced echelon of (aka RREF).</param> |
| | | 1721 | | /// <returns>The reduced echelon of the matrix (aka RREF).</returns> |
| | | 1722 | | public static Matrix<T> ReducedEchelon(Matrix<T> a) |
| | 4 | 1723 | | { |
| | 4 | 1724 | | Matrix<T>? b = null; |
| | 4 | 1725 | | ReducedEchelon(a, ref b); |
| | 4 | 1726 | | return b!; |
| | 4 | 1727 | | } |
| | | 1728 | | |
| | | 1729 | | /// <summary>Calculates the echelon of a matrix and reduces it (aka RREF).</summary> |
| | | 1730 | | /// <param name="b">The reduced echelon of the matrix (aka RREF).</param> |
| | | 1731 | | public void ReducedEchelon(ref Matrix<T>? b) |
| | 0 | 1732 | | { |
| | 0 | 1733 | | ReducedEchelon(this, ref b); |
| | 0 | 1734 | | } |
| | | 1735 | | |
| | | 1736 | | /// <summary>Matrixs the reduced echelon form of this matrix (aka RREF).</summary> |
| | | 1737 | | /// <returns>The computed reduced echelon form of this matrix (aka RREF).</returns> |
| | | 1738 | | public Matrix<T> ReducedEchelon() |
| | 4 | 1739 | | { |
| | 4 | 1740 | | return ReducedEchelon(this); |
| | 4 | 1741 | | } |
| | | 1742 | | |
| | | 1743 | | #endregion |
| | | 1744 | | |
| | | 1745 | | #region Inverse |
| | | 1746 | | |
| | | 1747 | | /// <summary>Calculates the inverse of a matrix.</summary> |
| | | 1748 | | /// <param name="a">The matrix to calculate the inverse of.</param> |
| | | 1749 | | /// <param name="b">The inverse of the matrix.</param> |
| | | 1750 | | public static void Inverse(Matrix<T> a, ref Matrix<T>? b) |
| | 2 | 1751 | | { |
| | 2 | 1752 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 2 | 1753 | | if (sourceof(!a.IsSquare, out string c1)) throw new ArgumentException(c1); |
| | 2 | 1754 | | T determinant = Determinant(a); |
| | 2 | 1755 | | if (Equate(determinant, Constant<T>.Zero)) |
| | 0 | 1756 | | { |
| | 0 | 1757 | | throw new ArgumentException("Singular matrix encountered during inverse caluculation (cannot be inversed)."); |
| | | 1758 | | } |
| | 2 | 1759 | | Matrix<T> adjoint = a.Adjoint(); |
| | 2 | 1760 | | int dimension = a.Rows; |
| | 2 | 1761 | | int Length = a.Length; |
| | 2 | 1762 | | if (object.ReferenceEquals(a, b)) |
| | 0 | 1763 | | { |
| | 0 | 1764 | | b = a.Clone(); |
| | 0 | 1765 | | } |
| | 2 | 1766 | | else if (b is not null && b.Length == Length) |
| | 0 | 1767 | | { |
| | 0 | 1768 | | b._rows = dimension; |
| | 0 | 1769 | | b._columns = dimension; |
| | 0 | 1770 | | } |
| | | 1771 | | else |
| | 2 | 1772 | | { |
| | 2 | 1773 | | b = new Matrix<T>(dimension, dimension, Length); |
| | 2 | 1774 | | } |
| | 16 | 1775 | | for (int i = 0; i < dimension; i++) |
| | 6 | 1776 | | { |
| | 48 | 1777 | | for (int j = 0; j < dimension; j++) |
| | 18 | 1778 | | { |
| | 18 | 1779 | | b.Set(i, j, Division(adjoint.Get(i, j), determinant)); |
| | 18 | 1780 | | } |
| | 6 | 1781 | | } |
| | | 1782 | | |
| | | 1783 | | #region Old Version |
| | | 1784 | | |
| | | 1785 | | //// Note: this method can be optimized... |
| | | 1786 | | |
| | | 1787 | | // if (a is null) |
| | | 1788 | | // { |
| | | 1789 | | // throw new ArgumentNullException(nameof(a)); |
| | | 1790 | | // } |
| | | 1791 | | // if (Compute.Equal(Determinant(a), Constant<T>.Zero)) |
| | | 1792 | | // { |
| | | 1793 | | // throw new MathematicsException("inverse calculation failed."); |
| | | 1794 | | // } |
| | | 1795 | | // Matrix<T> identity = FactoryIdentity(a.Rows, a.Columns); |
| | | 1796 | | // Matrix<T> rref = a.Clone(); |
| | | 1797 | | // for (int i = 0; i < a.Rows; i++) |
| | | 1798 | | // { |
| | | 1799 | | // if (Compute.Equal(rref[i, i], Constant<T>.Zero)) |
| | | 1800 | | // { |
| | | 1801 | | // for (int j = i + 1; j < rref.Rows; j++) |
| | | 1802 | | // { |
| | | 1803 | | // if (Compute.NotEqual(rref[j, i], Constant<T>.Zero)) |
| | | 1804 | | // { |
| | | 1805 | | // SwapRows(rref, i, j); |
| | | 1806 | | // SwapRows(identity, i, j); |
| | | 1807 | | // } |
| | | 1808 | | // } |
| | | 1809 | | // } |
| | | 1810 | | // T identityRowMultiplier = Compute.Divide(Constant<T>.One, rref[i, i]); |
| | | 1811 | | // RowMultiplication(identity, i, identityRowMultiplier); |
| | | 1812 | | // RowMultiplication(rref, i, identityRowMultiplier); |
| | | 1813 | | // for (int j = i + 1; j < rref.Rows; j++) |
| | | 1814 | | // { |
| | | 1815 | | // T rowAdder = Compute.Negate(rref[j, i]); |
| | | 1816 | | // RowAddition(identity, j, i, rowAdder); |
| | | 1817 | | // RowAddition(rref, j, i, rowAdder); |
| | | 1818 | | // } |
| | | 1819 | | // for (int j = i - 1; j >= 0; j--) |
| | | 1820 | | // { |
| | | 1821 | | // T rowAdder = Compute.Negate(rref[j, i]); |
| | | 1822 | | // RowAddition(identity, j, i, rowAdder); |
| | | 1823 | | // RowAddition(rref, j, i, rowAdder); |
| | | 1824 | | // } |
| | | 1825 | | // } |
| | | 1826 | | // b = identity; |
| | | 1827 | | |
| | | 1828 | | #endregion |
| | | 1829 | | |
| | | 1830 | | #region Alternate Version |
| | | 1831 | | |
| | | 1832 | | // Matrix<T> identity = Matrix<T>.FactoryIdentity(matrix.Rows, matrix.Columns); |
| | | 1833 | | // Matrix<T> rref = matrix.Clone(); |
| | | 1834 | | // for (int i = 0; i < matrix.Rows; i++) |
| | | 1835 | | // { |
| | | 1836 | | // if (Compute.Equate(rref[i, i], Compute.FromInt32(0))) |
| | | 1837 | | // for (int j = i + 1; j < rref.Rows; j++) |
| | | 1838 | | // if (!Compute.Equate(rref[j, i], Compute.FromInt32(0))) |
| | | 1839 | | // { |
| | | 1840 | | // Matrix<T>.SwapRows(rref, i, j); |
| | | 1841 | | // Matrix<T>.SwapRows(identity, i, j); |
| | | 1842 | | // } |
| | | 1843 | | // Matrix<T>.RowMultiplication(identity, i, Compute.Divide(Compute.FromInt32(1), rref[i, i])); |
| | | 1844 | | // Matrix<T>.RowMultiplication(rref, i, Compute.Divide(Compute.FromInt32(1), rref[i, i])); |
| | | 1845 | | // for (int j = i + 1; j < rref.Rows; j++) |
| | | 1846 | | // { |
| | | 1847 | | // Matrix<T>.RowAddition(identity, j, i, Compute.Negate(rref[j, i])); |
| | | 1848 | | // Matrix<T>.RowAddition(rref, j, i, Compute.Negate(rref[j, i])); |
| | | 1849 | | // } |
| | | 1850 | | // for (int j = i - 1; j >= 0; j--) |
| | | 1851 | | // { |
| | | 1852 | | // Matrix<T>.RowAddition(identity, j, i, Compute.Negate(rref[j, i])); |
| | | 1853 | | // Matrix<T>.RowAddition(rref, j, i, Compute.Negate(rref[j, i])); |
| | | 1854 | | // } |
| | | 1855 | | // } |
| | | 1856 | | // return identity; |
| | | 1857 | | |
| | | 1858 | | #endregion |
| | 2 | 1859 | | } |
| | | 1860 | | |
| | | 1861 | | /// <summary>Calculates the inverse of a matrix.</summary> |
| | | 1862 | | /// <param name="a">The matrix to calculate the inverse of.</param> |
| | | 1863 | | /// <returns>The inverse of the matrix.</returns> |
| | | 1864 | | public static Matrix<T> Inverse(Matrix<T> a) |
| | 2 | 1865 | | { |
| | 2 | 1866 | | Matrix<T>? b = null; |
| | 2 | 1867 | | Inverse(a, ref b); |
| | 2 | 1868 | | return b!; |
| | 2 | 1869 | | } |
| | | 1870 | | |
| | | 1871 | | /// <summary>Matrixs the inverse of this matrix.</summary> |
| | | 1872 | | /// <returns>The inverse of this matrix.</returns> |
| | | 1873 | | public Matrix<T> Inverse() |
| | 2 | 1874 | | { |
| | 2 | 1875 | | return Inverse(this); |
| | 2 | 1876 | | } |
| | | 1877 | | |
| | | 1878 | | #endregion |
| | | 1879 | | |
| | | 1880 | | #region Adjoint |
| | | 1881 | | |
| | | 1882 | | /// <summary>Calculates the adjoint of a matrix.</summary> |
| | | 1883 | | /// <param name="a">The matrix to calculate the adjoint of.</param> |
| | | 1884 | | /// <param name="b">The adjoint of the matrix.</param> |
| | | 1885 | | public static void Adjoint(Matrix<T> a, ref Matrix<T>? b) |
| | 6 | 1886 | | { |
| | 6 | 1887 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 6 | 1888 | | if (sourceof(!a.IsSquare, out string c1)) throw new ArgumentException(c1); |
| | 6 | 1889 | | int dimension = a.Rows; |
| | 6 | 1890 | | int Length = a.Length; |
| | 6 | 1891 | | if (ReferenceEquals(a, b)) |
| | 0 | 1892 | | { |
| | 0 | 1893 | | b = a.Clone(); |
| | 0 | 1894 | | } |
| | 6 | 1895 | | else if (b is not null && b.Length == Length) |
| | 0 | 1896 | | { |
| | 0 | 1897 | | b._rows = dimension; |
| | 0 | 1898 | | b._columns = dimension; |
| | 0 | 1899 | | } |
| | | 1900 | | else |
| | 6 | 1901 | | { |
| | 6 | 1902 | | b = new Matrix<T>(dimension, dimension, Length); |
| | 6 | 1903 | | } |
| | 6 | 1904 | | if (dimension is 1) |
| | 0 | 1905 | | { |
| | 0 | 1906 | | b.Set(0, 0, Constant<T>.One); |
| | 0 | 1907 | | return; |
| | | 1908 | | } |
| | 6 | 1909 | | Matrix<T> temp = new(dimension, dimension, Length); |
| | 48 | 1910 | | for (int i = 0; i < dimension; i++) |
| | 18 | 1911 | | { |
| | 144 | 1912 | | for (int j = 0; j < dimension; j++) |
| | 54 | 1913 | | { |
| | 54 | 1914 | | GetCofactor(a, temp, i, j, dimension); |
| | 54 | 1915 | | T sign = (i + j) % 2 is 0 ? Constant<T>.One : Constant<T>.NegativeOne; |
| | 54 | 1916 | | b.Set(j, i, Multiplication(sign, GetDeterminantGaussian(temp, dimension - 1))); |
| | 54 | 1917 | | } |
| | 18 | 1918 | | } |
| | | 1919 | | |
| | | 1920 | | #region Old Version |
| | | 1921 | | |
| | | 1922 | | // if (a is null) |
| | | 1923 | | // { |
| | | 1924 | | // throw new ArgumentNullException(nameof(a)); |
| | | 1925 | | // } |
| | | 1926 | | // if (!a.IsSquare) |
| | | 1927 | | // { |
| | | 1928 | | // throw new MathematicsException("Argument invalid !(" + nameof(a) + "." + nameof(a.IsSquare) + ")"); |
| | | 1929 | | // } |
| | | 1930 | | // if (object.ReferenceEquals(a, b)) |
| | | 1931 | | // { |
| | | 1932 | | // a = a.Clone(); |
| | | 1933 | | // } |
| | | 1934 | | // int Length = a.Length; |
| | | 1935 | | // int Rows = a.Rows; |
| | | 1936 | | // int Columns = a.Columns; |
| | | 1937 | | // if (b is not null && b.Length == Length) |
| | | 1938 | | // { |
| | | 1939 | | // b._rows = Rows; |
| | | 1940 | | // b._columns = Columns; |
| | | 1941 | | // } |
| | | 1942 | | // else |
| | | 1943 | | // { |
| | | 1944 | | // b = new Matrix<T>(Rows, Columns, Length); |
| | | 1945 | | // } |
| | | 1946 | | // for (int i = 0; i < Rows; i++) |
| | | 1947 | | // { |
| | | 1948 | | // for (int j = 0; j < Columns; j++) |
| | | 1949 | | // { |
| | | 1950 | | // if (Compute.IsEven(a.Get(i, j))) |
| | | 1951 | | // { |
| | | 1952 | | // b[i, j] = Determinant(Minor(a, i, j)); |
| | | 1953 | | // } |
| | | 1954 | | // else |
| | | 1955 | | // { |
| | | 1956 | | // b[i, j] = Compute.Negate(Determinant(Minor(a, i, j))); |
| | | 1957 | | // } |
| | | 1958 | | // } |
| | | 1959 | | // } |
| | | 1960 | | |
| | | 1961 | | #endregion |
| | 6 | 1962 | | } |
| | | 1963 | | |
| | | 1964 | | /// <summary>Calculates the adjoint of a matrix.</summary> |
| | | 1965 | | /// <param name="a">The matrix to calculate the adjoint of.</param> |
| | | 1966 | | /// <returns>The adjoint of the matrix.</returns> |
| | | 1967 | | public static Matrix<T> Adjoint(Matrix<T> a) |
| | 6 | 1968 | | { |
| | 6 | 1969 | | Matrix<T>? b = null; |
| | 6 | 1970 | | Adjoint(a, ref b); |
| | 6 | 1971 | | return b!; |
| | 6 | 1972 | | } |
| | | 1973 | | |
| | | 1974 | | /// <summary>Calculates the adjoint of a matrix.</summary> |
| | | 1975 | | /// <param name="b">The adjoint of the matrix.</param> |
| | | 1976 | | public void Adjoint(ref Matrix<T>? b) |
| | 0 | 1977 | | { |
| | 0 | 1978 | | Adjoint(this, ref b); |
| | 0 | 1979 | | } |
| | | 1980 | | |
| | | 1981 | | /// <summary>Calculates the adjoint of a matrix.</summary> |
| | | 1982 | | /// <returns>The adjoint of the matrix.</returns> |
| | | 1983 | | public Matrix<T> Adjoint() |
| | 6 | 1984 | | { |
| | 6 | 1985 | | return Adjoint(this); |
| | 6 | 1986 | | } |
| | | 1987 | | |
| | | 1988 | | #endregion |
| | | 1989 | | |
| | | 1990 | | #region Transpose |
| | | 1991 | | |
| | | 1992 | | /// <summary>Returns the transpose of a matrix.</summary> |
| | | 1993 | | /// <param name="a">The matrix to transpose.</param> |
| | | 1994 | | /// <param name="b">The transpose of the matrix.</param> |
| | | 1995 | | public static void Transpose(Matrix<T> a, ref Matrix<T>? b) |
| | 3 | 1996 | | { |
| | 3 | 1997 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 3 | 1998 | | if (ReferenceEquals(a, b)) |
| | 0 | 1999 | | { |
| | 0 | 2000 | | if (b.IsSquare) |
| | 0 | 2001 | | { |
| | 0 | 2002 | | TransposeContents(b); |
| | 0 | 2003 | | return; |
| | | 2004 | | } |
| | 0 | 2005 | | } |
| | 3 | 2006 | | int Length = a.Length; |
| | 3 | 2007 | | int Rows = a.Columns; |
| | 3 | 2008 | | int Columns = a.Rows; |
| | 3 | 2009 | | if (b is not null && b.Length == a.Length && !ReferenceEquals(a, b)) |
| | 0 | 2010 | | { |
| | 0 | 2011 | | b._rows = Rows; |
| | 0 | 2012 | | b._columns = Columns; |
| | 0 | 2013 | | } |
| | | 2014 | | else |
| | 3 | 2015 | | { |
| | 3 | 2016 | | b = new Matrix<T>(Rows, Columns, Length); |
| | 3 | 2017 | | } |
| | 18 | 2018 | | for (int i = 0; i < Rows; i++) |
| | 6 | 2019 | | { |
| | 34 | 2020 | | for (int j = 0; j < Columns; j++) |
| | 11 | 2021 | | { |
| | 11 | 2022 | | b.Set(i, j, a.Get(j, i)); |
| | 11 | 2023 | | } |
| | 6 | 2024 | | } |
| | 3 | 2025 | | } |
| | | 2026 | | |
| | | 2027 | | /// <summary>Returns the transpose of a matrix.</summary> |
| | | 2028 | | /// <param name="a">The matrix to transpose.</param> |
| | | 2029 | | /// <returns>The transpose of the matrix.</returns> |
| | | 2030 | | public static Matrix<T> Transpose(Matrix<T> a) |
| | 3 | 2031 | | { |
| | 3 | 2032 | | Matrix<T>? b = null; |
| | 3 | 2033 | | Transpose(a, ref b); |
| | 3 | 2034 | | return b!; |
| | 3 | 2035 | | } |
| | | 2036 | | |
| | | 2037 | | /// <summary>Returns the transpose of a matrix.</summary> |
| | | 2038 | | /// <param name="b">The transpose of the matrix.</param> |
| | | 2039 | | public void Transpose(ref Matrix<T>? b) |
| | 0 | 2040 | | { |
| | 0 | 2041 | | Transpose(this, ref b); |
| | 0 | 2042 | | } |
| | | 2043 | | |
| | | 2044 | | /// <summary>Returns the transpose of a matrix.</summary> |
| | | 2045 | | /// <returns>The transpose of the matrix.</returns> |
| | | 2046 | | public Matrix<T> Transpose() |
| | 3 | 2047 | | { |
| | 3 | 2048 | | return Transpose(this); |
| | 3 | 2049 | | } |
| | | 2050 | | |
| | | 2051 | | #endregion |
| | | 2052 | | |
| | | 2053 | | #region DecomposeLowerUpper |
| | | 2054 | | |
| | | 2055 | | /// <summary>Decomposes a matrix into lower-upper reptresentation.</summary> |
| | | 2056 | | /// <param name="matrix">The matrix to decompose.</param> |
| | | 2057 | | /// <param name="lower">The computed lower triangular matrix.</param> |
| | | 2058 | | /// <param name="upper">The computed upper triangular matrix.</param> |
| | | 2059 | | public static void DecomposeLowerUpper(Matrix<T> matrix, ref Matrix<T> lower, ref Matrix<T> upper) |
| | 0 | 2060 | | { |
| | | 2061 | | // Note: this method can be optimized... |
| | | 2062 | | |
| | 0 | 2063 | | lower = Matrix<T>.FactoryIdentity(matrix.Rows, matrix.Columns); |
| | 0 | 2064 | | upper = matrix.Clone(); |
| | 0 | 2065 | | int[] permutation = new int[matrix.Rows]; |
| | 0 | 2066 | | for (int i = 0; i < matrix.Rows; i++) |
| | 0 | 2067 | | { |
| | 0 | 2068 | | permutation[i] = i; |
| | 0 | 2069 | | } |
| | 0 | 2070 | | T pom2, detOfP = Constant<T>.One; |
| | 0 | 2071 | | int k0 = 0; |
| | 0 | 2072 | | for (int k = 0; k < matrix.Columns - 1; k++) |
| | 0 | 2073 | | { |
| | 0 | 2074 | | T p = Constant<T>.Zero; |
| | 0 | 2075 | | for (int i = k; i < matrix.Rows; i++) |
| | 0 | 2076 | | { |
| | 0 | 2077 | | if (GreaterThan(GreaterThan(upper[i, k], Constant<T>.Zero) ? upper[i, k] : Negation(upper[i, k]), p)) |
| | 0 | 2078 | | { |
| | 0 | 2079 | | p = GreaterThan(upper[i, k], Constant<T>.Zero) ? upper[i, k] : Negation(upper[i, k]); |
| | 0 | 2080 | | k0 = i; |
| | 0 | 2081 | | } |
| | 0 | 2082 | | } |
| | 0 | 2083 | | if (Equate(p, Constant<T>.Zero)) |
| | 0 | 2084 | | { |
| | 0 | 2085 | | throw new ArgumentException("The matrix is singular!"); |
| | | 2086 | | } |
| | 0 | 2087 | | (permutation[k0], permutation[k]) = (permutation[k], permutation[k0]); |
| | 0 | 2088 | | for (int i = 0; i < k; i++) |
| | 0 | 2089 | | { |
| | 0 | 2090 | | pom2 = lower[k, i]; |
| | 0 | 2091 | | lower[k, i] = lower[k0, i]; |
| | 0 | 2092 | | lower[k0, i] = pom2; |
| | 0 | 2093 | | } |
| | 0 | 2094 | | if (k != k0) |
| | 0 | 2095 | | { |
| | 0 | 2096 | | detOfP = Multiplication(detOfP, Constant<T>.NegativeOne); |
| | 0 | 2097 | | } |
| | 0 | 2098 | | for (int i = 0; i < matrix.Columns; i++) |
| | 0 | 2099 | | { |
| | 0 | 2100 | | pom2 = upper[k, i]; |
| | 0 | 2101 | | upper[k, i] = upper[k0, i]; |
| | 0 | 2102 | | upper[k0, i] = pom2; |
| | 0 | 2103 | | } |
| | 0 | 2104 | | for (int i = k + 1; i < matrix.Rows; i++) |
| | 0 | 2105 | | { |
| | 0 | 2106 | | lower[i, k] = Division(upper[i, k], upper[k, k]); |
| | 0 | 2107 | | for (int j = k; j < matrix.Columns; j++) |
| | 0 | 2108 | | { |
| | 0 | 2109 | | upper[i, j] = Subtraction(upper[i, j], Multiplication(lower[i, k], upper[k, j])); |
| | 0 | 2110 | | } |
| | 0 | 2111 | | } |
| | 0 | 2112 | | } |
| | | 2113 | | |
| | | 2114 | | #region Alternate Version |
| | | 2115 | | |
| | | 2116 | | // lower = Matrix<T>.FactoryIdentity(matrix.Rows, matrix.Columns); |
| | | 2117 | | // upper = matrix.Clone(); |
| | | 2118 | | // int[] permutation = new int[matrix.Rows]; |
| | | 2119 | | // for (int i = 0; i < matrix.Rows; i++) permutation[i] = i; |
| | | 2120 | | // T p = 0, pom2, detOfP = 1; |
| | | 2121 | | // int k0 = 0, pom1 = 0; |
| | | 2122 | | // for (int k = 0; k < matrix.Columns - 1; k++) |
| | | 2123 | | // { |
| | | 2124 | | // p = 0; |
| | | 2125 | | // for (int i = k; i < matrix.Rows; i++) |
| | | 2126 | | // if ((upper[i, k] > 0 ? upper[i, k] : -upper[i, k]) > p) |
| | | 2127 | | // { |
| | | 2128 | | // p = upper[i, k] > 0 ? upper[i, k] : -upper[i, k]; |
| | | 2129 | | // k0 = i; |
| | | 2130 | | // } |
| | | 2131 | | // if (p is 0) |
| | | 2132 | | // throw new System.Exception("The matrix is singular!"); |
| | | 2133 | | // pom1 = permutation[k]; |
| | | 2134 | | // permutation[k] = permutation[k0]; |
| | | 2135 | | // permutation[k0] = pom1; |
| | | 2136 | | // for (int i = 0; i < k; i++) |
| | | 2137 | | // { |
| | | 2138 | | // pom2 = lower[k, i]; |
| | | 2139 | | // lower[k, i] = lower[k0, i]; |
| | | 2140 | | // lower[k0, i] = pom2; |
| | | 2141 | | // } |
| | | 2142 | | // if (k != k0) |
| | | 2143 | | // detOfP *= -1; |
| | | 2144 | | // for (int i = 0; i < matrix.Columns; i++) |
| | | 2145 | | // { |
| | | 2146 | | // pom2 = upper[k, i]; |
| | | 2147 | | // upper[k, i] = upper[k0, i]; |
| | | 2148 | | // upper[k0, i] = pom2; |
| | | 2149 | | // } |
| | | 2150 | | // for (int i = k + 1; i < matrix.Rows; i++) |
| | | 2151 | | // { |
| | | 2152 | | // lower[i, k] = upper[i, k] / upper[k, k]; |
| | | 2153 | | // for (int j = k; j < matrix.Columns; j++) |
| | | 2154 | | // upper[i, j] = upper[i, j] - lower[i, k] * upper[k, j]; |
| | | 2155 | | // } |
| | | 2156 | | // } |
| | | 2157 | | |
| | | 2158 | | #endregion |
| | 0 | 2159 | | } |
| | | 2160 | | |
| | | 2161 | | /// <summary>Decomposes a matrix into lower-upper reptresentation.</summary> |
| | | 2162 | | /// <param name="lower">The computed lower triangular matrix.</param> |
| | | 2163 | | /// <param name="upper">The computed upper triangular matrix.</param> |
| | | 2164 | | public void DecomposeLowerUpper(ref Matrix<T> lower, ref Matrix<T> upper) |
| | 0 | 2165 | | { |
| | 0 | 2166 | | DecomposeLowerUpper(this, ref lower, ref upper); |
| | 0 | 2167 | | } |
| | | 2168 | | |
| | | 2169 | | #endregion |
| | | 2170 | | |
| | | 2171 | | #region Rotate |
| | | 2172 | | |
| | | 2173 | | /// <summary>Rotates a 4x4 matrix around an 3D axis by a specified angle.</summary> |
| | | 2174 | | /// <param name="matrix">The 4x4 matrix to rotate.</param> |
| | | 2175 | | /// <param name="angle">The angle of rotation around the axis.</param> |
| | | 2176 | | /// <param name="axis">The 3D axis to rotate the matrix around.</param> |
| | | 2177 | | /// <returns>The rotated matrix.</returns> |
| | | 2178 | | public static Matrix<T> Rotate4x4(Matrix<T> matrix, Angle<T> angle, Vector<T> axis) |
| | 0 | 2179 | | { |
| | 0 | 2180 | | if (axis is null) throw new ArgumentNullException(nameof(axis)); |
| | 0 | 2181 | | if (matrix is null) throw new ArgumentNullException(nameof(matrix)); |
| | 0 | 2182 | | if (sourceof(axis._vector.Length is not 3, out string c1)) throw new ArgumentException(c1); |
| | 0 | 2183 | | if (sourceof(matrix._rows is not 4 || matrix._columns is not 4, out string c2)) throw new ArgumentException(c2); |
| | | 2184 | | |
| | | 2185 | | // if the angle is zero, no rotation is required |
| | 0 | 2186 | | if (Equate(angle._measurement, Constant<T>.Zero)) |
| | 0 | 2187 | | { |
| | 0 | 2188 | | return matrix.Clone(); |
| | | 2189 | | } |
| | | 2190 | | |
| | 0 | 2191 | | T cosine = CosineSystem(angle); |
| | 0 | 2192 | | T sine = SineSystem(angle); |
| | 0 | 2193 | | T oneMinusCosine = Subtraction(Constant<T>.One, cosine); |
| | 0 | 2194 | | T xy = Multiplication(axis.X, axis.Y); |
| | 0 | 2195 | | T yz = Multiplication(axis.Y, axis.Z); |
| | 0 | 2196 | | T xz = Multiplication(axis.X, axis.Z); |
| | 0 | 2197 | | T xs = Multiplication(axis.X, sine); |
| | 0 | 2198 | | T ys = Multiplication(axis.Y, sine); |
| | 0 | 2199 | | T zs = Multiplication(axis.Z, sine); |
| | | 2200 | | |
| | 0 | 2201 | | T f00 = Addition(Multiplication(Multiplication(axis.X, axis.X), oneMinusCosine), cosine); |
| | 0 | 2202 | | T f01 = Addition(Multiplication(xy, oneMinusCosine), zs); |
| | 0 | 2203 | | T f02 = Subtraction(Multiplication(xz, oneMinusCosine), ys); |
| | | 2204 | | // n[3] not used |
| | 0 | 2205 | | T f10 = Subtraction(Multiplication(xy, oneMinusCosine), zs); |
| | 0 | 2206 | | T f11 = Addition(Multiplication(Multiplication(axis.Y, axis.Y), oneMinusCosine), cosine); |
| | 0 | 2207 | | T f12 = Addition(Multiplication(yz, oneMinusCosine), xs); |
| | | 2208 | | // n[7] not used |
| | 0 | 2209 | | T f20 = Addition(Multiplication(xz, oneMinusCosine), ys); |
| | 0 | 2210 | | T f21 = Subtraction(Multiplication(yz, oneMinusCosine), xs); |
| | 0 | 2211 | | T f22 = Addition(Multiplication(Multiplication(axis.Z, axis.Z), oneMinusCosine), cosine); |
| | | 2212 | | |
| | | 2213 | | // Row 1 |
| | 0 | 2214 | | T _0_0 = Addition(Addition(Multiplication(matrix[0, 0], f00), Multiplication(matrix[1, 0], f01)), Multiplication(mat |
| | 0 | 2215 | | T _0_1 = Addition(Addition(Multiplication(matrix[0, 1], f00), Multiplication(matrix[1, 1], f01)), Multiplication(mat |
| | 0 | 2216 | | T _0_2 = Addition(Addition(Multiplication(matrix[0, 2], f00), Multiplication(matrix[1, 2], f01)), Multiplication(mat |
| | 0 | 2217 | | T _0_3 = Addition(Addition(Multiplication(matrix[0, 3], f00), Multiplication(matrix[1, 3], f01)), Multiplication(mat |
| | | 2218 | | // Row 2 |
| | 0 | 2219 | | T _1_0 = Addition(Addition(Multiplication(matrix[0, 0], f10), Multiplication(matrix[1, 0], f11)), Multiplication(mat |
| | 0 | 2220 | | T _1_1 = Addition(Addition(Multiplication(matrix[0, 1], f10), Multiplication(matrix[1, 1], f11)), Multiplication(mat |
| | 0 | 2221 | | T _1_2 = Addition(Addition(Multiplication(matrix[0, 2], f10), Multiplication(matrix[1, 2], f11)), Multiplication(mat |
| | 0 | 2222 | | T _1_3 = Addition(Addition(Multiplication(matrix[0, 3], f10), Multiplication(matrix[1, 3], f11)), Multiplication(mat |
| | | 2223 | | // Row 3 |
| | 0 | 2224 | | T _2_0 = Addition(Addition(Multiplication(matrix[0, 0], f20), Multiplication(matrix[1, 0], f21)), Multiplication(mat |
| | 0 | 2225 | | T _2_1 = Addition(Addition(Multiplication(matrix[0, 1], f20), Multiplication(matrix[1, 1], f21)), Multiplication(mat |
| | 0 | 2226 | | T _2_2 = Addition(Addition(Multiplication(matrix[0, 2], f20), Multiplication(matrix[1, 2], f21)), Multiplication(mat |
| | 0 | 2227 | | T _2_3 = Addition(Addition(Multiplication(matrix[0, 3], f20), Multiplication(matrix[1, 3], f21)), Multiplication(mat |
| | | 2228 | | // Row 4 |
| | 0 | 2229 | | T _3_0 = Constant<T>.Zero; |
| | 0 | 2230 | | T _3_1 = Constant<T>.Zero; |
| | 0 | 2231 | | T _3_2 = Constant<T>.Zero; |
| | 0 | 2232 | | T _3_3 = Constant<T>.One; |
| | | 2233 | | |
| | | 2234 | | #pragma warning disable CS8509 // The switch expression does not handle all possible values of its input type (it is not |
| | 0 | 2235 | | return new Matrix<T>(4, 4, (row, column) => |
| | 0 | 2236 | | (row, column) switch |
| | 0 | 2237 | | { |
| | 0 | 2238 | | /* 0 */ (0, 0) => _0_0, (0, 1) => _0_1, (0, 2) => _0_2, (0, 3) => _0_3, |
| | 0 | 2239 | | /* 1 */ (1, 0) => _1_0, (1, 1) => _1_1, (1, 2) => _1_2, (1, 3) => _1_3, |
| | 0 | 2240 | | /* 2 */ (2, 0) => _2_0, (2, 1) => _2_1, (2, 2) => _2_2, (2, 3) => _2_3, |
| | 0 | 2241 | | /* 3 */ (3, 0) => _3_0, (3, 1) => _3_1, (3, 2) => _3_2, (3, 3) => _3_3, |
| | 0 | 2242 | | }); |
| | | 2243 | | #pragma warning restore CS8509 // The switch expression does not handle all possible values of its input type (it is not |
| | 0 | 2244 | | } |
| | | 2245 | | |
| | | 2246 | | /// <summary>Converts an angle around an axis into a 4x4 rotational matrix.</summary> |
| | | 2247 | | /// <param name="angle">The angle of the axis rotation.</param> |
| | | 2248 | | /// <param name="axis">The axis of the axis rotation.</param> |
| | | 2249 | | /// <returns>The rotation expressed as a 4x4 matrix.</returns> |
| | 0 | 2250 | | public Matrix<T> Rotate4x4(Angle<T> angle, Vector<T> axis) => Rotate4x4(this, angle, axis); |
| | | 2251 | | |
| | | 2252 | | #endregion |
| | | 2253 | | |
| | | 2254 | | #region Equal |
| | | 2255 | | |
| | | 2256 | | /// <summary>Does a value equality check.</summary> |
| | | 2257 | | /// <param name="a">The first matrix to check for equality.</param> |
| | | 2258 | | /// <param name="b">The second matrix to check for equality.</param> |
| | | 2259 | | /// <returns>True if values are equal, false if not.</returns> |
| | | 2260 | | public static bool Equal(Matrix<T> a, Matrix<T> b) |
| | 67 | 2261 | | { |
| | 67 | 2262 | | if (a is null) |
| | 0 | 2263 | | { |
| | 0 | 2264 | | if (b is null) |
| | 0 | 2265 | | { |
| | 0 | 2266 | | return true; |
| | | 2267 | | } |
| | | 2268 | | else |
| | 0 | 2269 | | { |
| | 0 | 2270 | | return false; |
| | | 2271 | | } |
| | | 2272 | | } |
| | 67 | 2273 | | if (b is null) |
| | 0 | 2274 | | { |
| | 0 | 2275 | | return false; |
| | | 2276 | | } |
| | 67 | 2277 | | int Rows = a.Rows; |
| | 67 | 2278 | | int Columns = a.Columns; |
| | 67 | 2279 | | if (Rows != b.Rows || Columns != b.Columns) |
| | 0 | 2280 | | { |
| | 0 | 2281 | | return false; |
| | | 2282 | | } |
| | 67 | 2283 | | T[] A = a._matrix; |
| | 67 | 2284 | | T[] B = b._matrix; |
| | 67 | 2285 | | int Length = A.Length; |
| | 1074 | 2286 | | for (int i = 0; i < Length; i++) |
| | 470 | 2287 | | { |
| | 470 | 2288 | | if (!Equate(A[i], B[i])) |
| | 0 | 2289 | | { |
| | 0 | 2290 | | return false; |
| | | 2291 | | } |
| | 470 | 2292 | | } |
| | 67 | 2293 | | return true; |
| | 67 | 2294 | | } |
| | | 2295 | | |
| | | 2296 | | /// <summary>Does a value equality check.</summary> |
| | | 2297 | | /// <param name="a">The first matrix to check for equality.</param> |
| | | 2298 | | /// <param name="b">The second matrix to check for equality.</param> |
| | | 2299 | | /// <returns>True if values are equal, false if not.</returns> |
| | | 2300 | | public static bool operator ==(Matrix<T> a, Matrix<T> b) |
| | 67 | 2301 | | { |
| | 67 | 2302 | | return Equal(a, b); |
| | 67 | 2303 | | } |
| | | 2304 | | |
| | | 2305 | | /// <summary>Does a value non-equality check.</summary> |
| | | 2306 | | /// <param name="a">The first matrix to check for non-equality.</param> |
| | | 2307 | | /// <param name="b">The second matrix to check for non-equality.</param> |
| | | 2308 | | /// <returns>True if values are not equal, false if not.</returns> |
| | | 2309 | | public static bool operator !=(Matrix<T> a, Matrix<T> b) |
| | 0 | 2310 | | { |
| | 0 | 2311 | | return !Equal(a, b); |
| | 0 | 2312 | | } |
| | | 2313 | | |
| | | 2314 | | /// <summary>Does a value equality check.</summary> |
| | | 2315 | | /// <param name="b">The second matrix to check for equality.</param> |
| | | 2316 | | /// <returns>True if values are equal, false if not.</returns> |
| | | 2317 | | public bool Equal(Matrix<T> b) |
| | 0 | 2318 | | { |
| | 0 | 2319 | | return this == b; |
| | 0 | 2320 | | } |
| | | 2321 | | |
| | | 2322 | | #endregion |
| | | 2323 | | |
| | | 2324 | | #region Equal (+leniency) |
| | | 2325 | | |
| | | 2326 | | /// <summary>Does a value equality check with leniency.</summary> |
| | | 2327 | | /// <param name="a">The first matrix to check for equality.</param> |
| | | 2328 | | /// <param name="b">The second matrix to check for equality.</param> |
| | | 2329 | | /// <param name="leniency">How much the values can vary but still be considered equal.</param> |
| | | 2330 | | /// <returns>True if values are equal, false if not.</returns> |
| | | 2331 | | public static bool Equal(Matrix<T> a, Matrix<T> b, T leniency) |
| | 10 | 2332 | | { |
| | 10 | 2333 | | if (a is null) |
| | 0 | 2334 | | { |
| | 0 | 2335 | | if (b is null) |
| | 0 | 2336 | | { |
| | 0 | 2337 | | return true; |
| | | 2338 | | } |
| | | 2339 | | else |
| | 0 | 2340 | | { |
| | 0 | 2341 | | return false; |
| | | 2342 | | } |
| | | 2343 | | } |
| | 10 | 2344 | | if (b is null) |
| | 0 | 2345 | | { |
| | 0 | 2346 | | return false; |
| | | 2347 | | } |
| | 10 | 2348 | | int Rows = a.Rows; |
| | 10 | 2349 | | int Columns = a.Columns; |
| | 10 | 2350 | | if (Rows != b.Rows || Columns != b.Columns) |
| | 0 | 2351 | | { |
| | 0 | 2352 | | return false; |
| | | 2353 | | } |
| | 10 | 2354 | | T[] A = a._matrix; |
| | 10 | 2355 | | T[] B = b._matrix; |
| | 10 | 2356 | | int Length = A.Length; |
| | 128 | 2357 | | for (int i = 0; i < Length; i++) |
| | 58 | 2358 | | { |
| | 58 | 2359 | | if (!EqualToLeniency(A[i], B[i], leniency)) |
| | 4 | 2360 | | { |
| | 4 | 2361 | | return false; |
| | | 2362 | | } |
| | 54 | 2363 | | } |
| | 6 | 2364 | | return true; |
| | 10 | 2365 | | } |
| | | 2366 | | |
| | | 2367 | | /// <summary>Does a value equality check with leniency.</summary> |
| | | 2368 | | /// <param name="b">The second matrix to check for equality.</param> |
| | | 2369 | | /// <param name="leniency">How much the values can vary but still be considered equal.</param> |
| | | 2370 | | /// <returns>True if values are equal, false if not.</returns> |
| | | 2371 | | public bool Equal(Matrix<T> b, T leniency) |
| | 10 | 2372 | | { |
| | 10 | 2373 | | return Equal(this, b, leniency); |
| | 10 | 2374 | | } |
| | | 2375 | | |
| | | 2376 | | #endregion |
| | | 2377 | | |
| | | 2378 | | #endregion |
| | | 2379 | | |
| | | 2380 | | #region Other Methods |
| | | 2381 | | |
| | | 2382 | | #region Get/Set |
| | | 2383 | | |
| | | 2384 | | internal T Get(int row, int column) |
| | 2009 | 2385 | | { |
| | 2009 | 2386 | | return _matrix[row * Columns + column]; |
| | 2009 | 2387 | | } |
| | | 2388 | | |
| | | 2389 | | internal void Set(int row, int column, T value) |
| | 912 | 2390 | | { |
| | 912 | 2391 | | _matrix[row * Columns + column] = value; |
| | 912 | 2392 | | } |
| | | 2393 | | |
| | | 2394 | | #endregion |
| | | 2395 | | |
| | | 2396 | | #region Format |
| | | 2397 | | |
| | | 2398 | | /// <summary>Fills a matrix with values using a delegate.</summary> |
| | | 2399 | | /// <param name="matrix">The matrix to fill the values of.</param> |
| | | 2400 | | /// <param name="function">The function to set the values at the relative indeces.</param> |
| | | 2401 | | public static void Format(Matrix<T> matrix, Func<int, int, T> function) |
| | 237 | 2402 | | { |
| | 237 | 2403 | | int Rows = matrix.Rows; |
| | 237 | 2404 | | int Columns = matrix.Columns; |
| | 237 | 2405 | | T[] MATRIX = matrix._matrix; |
| | 237 | 2406 | | int i = 0; |
| | 1764 | 2407 | | for (int row = 0; row < Rows; row++) |
| | 645 | 2408 | | { |
| | 4994 | 2409 | | for (int column = 0; column < Columns; column++) |
| | 1852 | 2410 | | { |
| | 1852 | 2411 | | MATRIX[i++] = function(row, column); |
| | 1852 | 2412 | | } |
| | 645 | 2413 | | } |
| | 237 | 2414 | | } |
| | | 2415 | | |
| | | 2416 | | /// <summary>Fills a matrix with values using a delegate.</summary> |
| | | 2417 | | /// <param name="matrix">The matrix to fill the values of.</param> |
| | | 2418 | | /// <param name="func">The function to set the values at the relative indeces.</param> |
| | | 2419 | | public static void Format(Matrix<T> matrix, Func<int, T> func) |
| | 0 | 2420 | | { |
| | 0 | 2421 | | matrix._matrix.Format(func); |
| | 0 | 2422 | | } |
| | | 2423 | | |
| | | 2424 | | #endregion |
| | | 2425 | | |
| | | 2426 | | #region CloneContents |
| | | 2427 | | |
| | | 2428 | | internal static void CloneContents(Matrix<T> a, Matrix<T> b) |
| | 0 | 2429 | | { |
| | 0 | 2430 | | T[] a_flat = a._matrix; |
| | 0 | 2431 | | T[] b_flat = b._matrix; |
| | 0 | 2432 | | int length = a_flat.Length; |
| | 0 | 2433 | | for (int i = 0; i < length; i++) |
| | 0 | 2434 | | { |
| | 0 | 2435 | | b_flat[i] = a_flat[i]; |
| | 0 | 2436 | | } |
| | 0 | 2437 | | } |
| | | 2438 | | |
| | | 2439 | | #endregion |
| | | 2440 | | |
| | | 2441 | | #region TransposeContents |
| | | 2442 | | |
| | | 2443 | | internal static void TransposeContents(Matrix<T> a) |
| | 0 | 2444 | | { |
| | 0 | 2445 | | int Rows = a.Rows; |
| | 0 | 2446 | | for (int i = 0; i < Rows; i++) |
| | 0 | 2447 | | { |
| | 0 | 2448 | | int Rows_Minus_i = Rows - i; |
| | 0 | 2449 | | for (int j = 0; j < Rows_Minus_i; j++) |
| | 0 | 2450 | | { |
| | 0 | 2451 | | T temp = a.Get(i, j); |
| | 0 | 2452 | | a.Set(i, j, a.Get(j, i)); |
| | 0 | 2453 | | a.Set(j, i, temp); |
| | 0 | 2454 | | } |
| | 0 | 2455 | | } |
| | 0 | 2456 | | } |
| | | 2457 | | |
| | | 2458 | | #endregion |
| | | 2459 | | |
| | | 2460 | | #region Clone |
| | | 2461 | | |
| | | 2462 | | /// <summary>Creates a copy of a matrix.</summary> |
| | | 2463 | | /// <param name="a">The matrix to copy.</param> |
| | | 2464 | | /// <returns>The copy of this matrix.</returns> |
| | | 2465 | | public static Matrix<T> Clone(Matrix<T> a) |
| | 16 | 2466 | | { |
| | 16 | 2467 | | if (a is null) throw new ArgumentNullException(nameof(a)); |
| | 16 | 2468 | | return new Matrix<T>(a); |
| | 16 | 2469 | | } |
| | | 2470 | | |
| | | 2471 | | /// <summary>Copies this matrix.</summary> |
| | | 2472 | | /// <returns>The copy of this matrix.</returns> |
| | | 2473 | | public Matrix<T> Clone() |
| | 16 | 2474 | | { |
| | 16 | 2475 | | return Clone(this); |
| | 16 | 2476 | | } |
| | | 2477 | | |
| | | 2478 | | #endregion |
| | | 2479 | | |
| | | 2480 | | #endregion |
| | | 2481 | | |
| | | 2482 | | #region Casting Operators |
| | | 2483 | | |
| | | 2484 | | /// <summary>Converts a T[,] into a matrix.</summary> |
| | | 2485 | | /// <param name="array">The T[,] to convert to a matrix.</param> |
| | | 2486 | | /// <returns>The resulting matrix after conversion.</returns> |
| | | 2487 | | public static implicit operator Matrix<T>(T[,] array) => |
| | 1469 | 2488 | | new(array.GetLength(0), array.GetLength(1), (i, j) => array[i, j]); |
| | | 2489 | | |
| | | 2490 | | /// <summary>Converts a matrix into a T[,].</summary> |
| | | 2491 | | /// <param name="matrix">The matrix toconvert to a T[,].</param> |
| | | 2492 | | /// <returns>The resulting T[,] after conversion.</returns> |
| | | 2493 | | public static explicit operator T[,](Matrix<T> matrix) |
| | 0 | 2494 | | { |
| | 0 | 2495 | | int rows = matrix._rows; |
| | 0 | 2496 | | int columns = matrix._columns; |
| | 0 | 2497 | | T[,] array = new T[rows, columns]; |
| | 0 | 2498 | | T[] MATRIX = matrix._matrix; |
| | 0 | 2499 | | int k = 0; |
| | 0 | 2500 | | for (int i = 0; i < rows; i++) |
| | 0 | 2501 | | { |
| | 0 | 2502 | | for (int j = 0; j < columns; j++) |
| | 0 | 2503 | | { |
| | 0 | 2504 | | array[i, j] = MATRIX[k++]; |
| | 0 | 2505 | | } |
| | 0 | 2506 | | } |
| | 0 | 2507 | | return array; |
| | 0 | 2508 | | } |
| | | 2509 | | |
| | | 2510 | | #endregion |
| | | 2511 | | |
| | | 2512 | | #region Overrides |
| | | 2513 | | |
| | | 2514 | | /// <summary>Matrixs a hash code from the values of this matrix.</summary> |
| | | 2515 | | /// <returns>A hash code for the matrix.</returns> |
| | | 2516 | | public override int GetHashCode() |
| | 0 | 2517 | | { |
| | 0 | 2518 | | int hashCode = HashCode.Combine(_rows, _columns); |
| | 0 | 2519 | | foreach (T value in _matrix) |
| | 0 | 2520 | | { |
| | 0 | 2521 | | hashCode = HashCode.Combine(hashCode, Hash(value)); |
| | 0 | 2522 | | } |
| | 0 | 2523 | | return hashCode; |
| | 0 | 2524 | | } |
| | | 2525 | | |
| | | 2526 | | /// <summary>Does an equality check by value.</summary> |
| | | 2527 | | /// <param name="b">The object to compare to.</param> |
| | | 2528 | | /// <returns>True if the references are equal, false if not.</returns> |
| | 0 | 2529 | | public override bool Equals(object? b) => b is Matrix<T> B && Equal(this, B); |
| | | 2530 | | |
| | | 2531 | | #endregion |
| | | 2532 | | } |