inv.js 6.0 KB

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