inv.js 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184
  1. import { isMatrix } from '../../utils/is.js';
  2. import { arraySize } from '../../utils/array.js';
  3. import { factory } from '../../utils/factory.js';
  4. import { format } from '../../utils/string.js';
  5. var name = 'inv';
  6. var dependencies = ['typed', 'matrix', 'divideScalar', 'addScalar', 'multiply', 'unaryMinus', 'det', 'identity', 'abs'];
  7. export var createInv = /* #__PURE__ */factory(name, dependencies, _ref => {
  8. var {
  9. typed,
  10. matrix,
  11. divideScalar,
  12. addScalar,
  13. multiply,
  14. unaryMinus,
  15. det,
  16. identity,
  17. abs
  18. } = _ref;
  19. /**
  20. * Calculate the inverse of a square matrix.
  21. *
  22. * Syntax:
  23. *
  24. * math.inv(x)
  25. *
  26. * Examples:
  27. *
  28. * math.inv([[1, 2], [3, 4]]) // returns [[-2, 1], [1.5, -0.5]]
  29. * math.inv(4) // returns 0.25
  30. * 1 / 4 // returns 0.25
  31. *
  32. * See also:
  33. *
  34. * det, transpose
  35. *
  36. * @param {number | Complex | Array | Matrix} x Matrix to be inversed
  37. * @return {number | Complex | Array | Matrix} The inverse of `x`.
  38. */
  39. return typed(name, {
  40. 'Array | Matrix': function ArrayMatrix(x) {
  41. var size = isMatrix(x) ? x.size() : arraySize(x);
  42. switch (size.length) {
  43. case 1:
  44. // vector
  45. if (size[0] === 1) {
  46. if (isMatrix(x)) {
  47. return matrix([divideScalar(1, x.valueOf()[0])]);
  48. } else {
  49. return [divideScalar(1, x[0])];
  50. }
  51. } else {
  52. throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
  53. }
  54. case 2:
  55. // two dimensional array
  56. {
  57. var rows = size[0];
  58. var cols = size[1];
  59. if (rows === cols) {
  60. if (isMatrix(x)) {
  61. return matrix(_inv(x.valueOf(), rows, cols), x.storage());
  62. } else {
  63. // return an Array
  64. return _inv(x, rows, cols);
  65. }
  66. } else {
  67. throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
  68. }
  69. }
  70. default:
  71. // multi dimensional array
  72. throw new RangeError('Matrix must be two dimensional ' + '(size: ' + format(size) + ')');
  73. }
  74. },
  75. any: function any(x) {
  76. // scalar
  77. return divideScalar(1, x); // FIXME: create a BigNumber one when configured for bignumbers
  78. }
  79. });
  80. /**
  81. * Calculate the inverse of a square matrix
  82. * @param {Array[]} mat A square matrix
  83. * @param {number} rows Number of rows
  84. * @param {number} cols Number of columns, must equal rows
  85. * @return {Array[]} inv Inverse matrix
  86. * @private
  87. */
  88. function _inv(mat, rows, cols) {
  89. var r, s, f, value, temp;
  90. if (rows === 1) {
  91. // this is a 1 x 1 matrix
  92. value = mat[0][0];
  93. if (value === 0) {
  94. throw Error('Cannot calculate inverse, determinant is zero');
  95. }
  96. return [[divideScalar(1, value)]];
  97. } else if (rows === 2) {
  98. // this is a 2 x 2 matrix
  99. var d = det(mat);
  100. if (d === 0) {
  101. throw Error('Cannot calculate inverse, determinant is zero');
  102. }
  103. return [[divideScalar(mat[1][1], d), divideScalar(unaryMinus(mat[0][1]), d)], [divideScalar(unaryMinus(mat[1][0]), d), divideScalar(mat[0][0], d)]];
  104. } else {
  105. // this is a matrix of 3 x 3 or larger
  106. // calculate inverse using gauss-jordan elimination
  107. // https://en.wikipedia.org/wiki/Gaussian_elimination
  108. // http://mathworld.wolfram.com/MatrixInverse.html
  109. // http://math.uww.edu/~mcfarlat/inverse.htm
  110. // make a copy of the matrix (only the arrays, not of the elements)
  111. var A = mat.concat();
  112. for (r = 0; r < rows; r++) {
  113. A[r] = A[r].concat();
  114. }
  115. // create an identity matrix which in the end will contain the
  116. // matrix inverse
  117. var B = identity(rows).valueOf();
  118. // loop over all columns, and perform row reductions
  119. for (var c = 0; c < cols; c++) {
  120. // Pivoting: Swap row c with row r, where row r contains the largest element A[r][c]
  121. var ABig = abs(A[c][c]);
  122. var rBig = c;
  123. r = c + 1;
  124. while (r < rows) {
  125. if (abs(A[r][c]) > ABig) {
  126. ABig = abs(A[r][c]);
  127. rBig = r;
  128. }
  129. r++;
  130. }
  131. if (ABig === 0) {
  132. throw Error('Cannot calculate inverse, determinant is zero');
  133. }
  134. r = rBig;
  135. if (r !== c) {
  136. temp = A[c];
  137. A[c] = A[r];
  138. A[r] = temp;
  139. temp = B[c];
  140. B[c] = B[r];
  141. B[r] = temp;
  142. }
  143. // eliminate non-zero values on the other rows at column c
  144. var Ac = A[c];
  145. var Bc = B[c];
  146. for (r = 0; r < rows; r++) {
  147. var Ar = A[r];
  148. var Br = B[r];
  149. if (r !== c) {
  150. // eliminate value at column c and row r
  151. if (Ar[c] !== 0) {
  152. f = divideScalar(unaryMinus(Ar[c]), Ac[c]);
  153. // add (f * row c) to row r to eliminate the value
  154. // at column c
  155. for (s = c; s < cols; s++) {
  156. Ar[s] = addScalar(Ar[s], multiply(f, Ac[s]));
  157. }
  158. for (s = 0; s < cols; s++) {
  159. Br[s] = addScalar(Br[s], multiply(f, Bc[s]));
  160. }
  161. }
  162. } else {
  163. // normalize value at Acc to 1,
  164. // divide each value on row r with the value at Acc
  165. f = Ac[c];
  166. for (s = c; s < cols; s++) {
  167. Ar[s] = divideScalar(Ar[s], f);
  168. }
  169. for (s = 0; s < cols; s++) {
  170. Br[s] = divideScalar(Br[s], f);
  171. }
  172. }
  173. }
  174. }
  175. return B;
  176. }
  177. }
  178. });