det.js 4.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143
  1. import { isMatrix } from '../../utils/is.js';
  2. import { clone } from '../../utils/object.js';
  3. import { format } from '../../utils/string.js';
  4. import { factory } from '../../utils/factory.js';
  5. var name = 'det';
  6. var dependencies = ['typed', 'matrix', 'subtract', 'multiply', 'divideScalar', 'isZero', 'unaryMinus'];
  7. export var createDet = /* #__PURE__ */factory(name, dependencies, _ref => {
  8. var {
  9. typed,
  10. matrix,
  11. subtract,
  12. multiply,
  13. divideScalar,
  14. isZero,
  15. unaryMinus
  16. } = _ref;
  17. /**
  18. * Calculate the determinant of a matrix.
  19. *
  20. * Syntax:
  21. *
  22. * math.det(x)
  23. *
  24. * Examples:
  25. *
  26. * math.det([[1, 2], [3, 4]]) // returns -2
  27. *
  28. * const A = [
  29. * [-2, 2, 3],
  30. * [-1, 1, 3],
  31. * [2, 0, -1]
  32. * ]
  33. * math.det(A) // returns 6
  34. *
  35. * See also:
  36. *
  37. * inv
  38. *
  39. * @param {Array | Matrix} x A matrix
  40. * @return {number} The determinant of `x`
  41. */
  42. return typed(name, {
  43. any: function any(x) {
  44. return clone(x);
  45. },
  46. 'Array | Matrix': function det(x) {
  47. var size;
  48. if (isMatrix(x)) {
  49. size = x.size();
  50. } else if (Array.isArray(x)) {
  51. x = matrix(x);
  52. size = x.size();
  53. } else {
  54. // a scalar
  55. size = [];
  56. }
  57. switch (size.length) {
  58. case 0:
  59. // scalar
  60. return clone(x);
  61. case 1:
  62. // vector
  63. if (size[0] === 1) {
  64. return clone(x.valueOf()[0]);
  65. }
  66. if (size[0] === 0) {
  67. return 1; // det of an empty matrix is per definition 1
  68. } else {
  69. throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
  70. }
  71. case 2:
  72. {
  73. // two-dimensional array
  74. var rows = size[0];
  75. var cols = size[1];
  76. if (rows === cols) {
  77. return _det(x.clone().valueOf(), rows, cols);
  78. }
  79. if (cols === 0) {
  80. return 1; // det of an empty matrix is per definition 1
  81. } else {
  82. throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
  83. }
  84. }
  85. default:
  86. // multi dimensional array
  87. throw new RangeError('Matrix must be two dimensional ' + '(size: ' + format(size) + ')');
  88. }
  89. }
  90. });
  91. /**
  92. * Calculate the determinant of a matrix
  93. * @param {Array[]} matrix A square, two dimensional matrix
  94. * @param {number} rows Number of rows of the matrix (zero-based)
  95. * @param {number} cols Number of columns of the matrix (zero-based)
  96. * @returns {number} det
  97. * @private
  98. */
  99. function _det(matrix, rows, cols) {
  100. if (rows === 1) {
  101. // this is a 1 x 1 matrix
  102. return clone(matrix[0][0]);
  103. } else if (rows === 2) {
  104. // this is a 2 x 2 matrix
  105. // the determinant of [a11,a12;a21,a22] is det = a11*a22-a21*a12
  106. return subtract(multiply(matrix[0][0], matrix[1][1]), multiply(matrix[1][0], matrix[0][1]));
  107. } else {
  108. // Bareiss algorithm
  109. // this algorithm have same complexity as LUP decomposition (O(n^3))
  110. // but it preserve precision of floating point more relative to the LUP decomposition
  111. var negated = false;
  112. var rowIndices = new Array(rows).fill(0).map((_, i) => i); // matrix index of row i
  113. for (var k = 0; k < rows; k++) {
  114. var k_ = rowIndices[k];
  115. if (isZero(matrix[k_][k])) {
  116. var _k = void 0;
  117. for (_k = k + 1; _k < rows; _k++) {
  118. if (!isZero(matrix[rowIndices[_k]][k])) {
  119. k_ = rowIndices[_k];
  120. rowIndices[_k] = rowIndices[k];
  121. rowIndices[k] = k_;
  122. negated = !negated;
  123. break;
  124. }
  125. }
  126. if (_k === rows) return matrix[k_][k]; // some zero of the type
  127. }
  128. var piv = matrix[k_][k];
  129. var piv_ = k === 0 ? 1 : matrix[rowIndices[k - 1]][k - 1];
  130. for (var i = k + 1; i < rows; i++) {
  131. var i_ = rowIndices[i];
  132. for (var j = k + 1; j < rows; j++) {
  133. matrix[i_][j] = divideScalar(subtract(multiply(matrix[i_][j], piv), multiply(matrix[i_][k], matrix[k_][j])), piv_);
  134. }
  135. }
  136. }
  137. var det = matrix[rowIndices[rows - 1]][rows - 1];
  138. return negated ? unaryMinus(det) : det;
  139. }
  140. }
  141. });