simple-squiggle

A restricted subset of Squiggle
Log | Files | Refs | README

det.js (3787B)

      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', 'unaryMinus', 'lup'];
      7 export var createDet = /* #__PURE__ */factory(name, dependencies, _ref => {
      8   var {
      9     typed,
     10     matrix,
     11     subtract,
     12     multiply,
     13     unaryMinus,
     14     lup
     15   } = _ref;
     16 
     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 
     49       if (isMatrix(x)) {
     50         size = x.size();
     51       } else if (Array.isArray(x)) {
     52         x = matrix(x);
     53         size = x.size();
     54       } else {
     55         // a scalar
     56         size = [];
     57       }
     58 
     59       switch (size.length) {
     60         case 0:
     61           // scalar
     62           return clone(x);
     63 
     64         case 1:
     65           // vector
     66           if (size[0] === 1) {
     67             return clone(x.valueOf()[0]);
     68           } else {
     69             throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
     70           }
     71 
     72         case 2:
     73           {
     74             // two dimensional array
     75             var rows = size[0];
     76             var cols = size[1];
     77 
     78             if (rows === cols) {
     79               return _det(x.clone().valueOf(), rows, cols);
     80             } else {
     81               throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
     82             }
     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 
    100   function _det(matrix, rows, cols) {
    101     if (rows === 1) {
    102       // this is a 1 x 1 matrix
    103       return clone(matrix[0][0]);
    104     } else if (rows === 2) {
    105       // this is a 2 x 2 matrix
    106       // the determinant of [a11,a12;a21,a22] is det = a11*a22-a21*a12
    107       return subtract(multiply(matrix[0][0], matrix[1][1]), multiply(matrix[1][0], matrix[0][1]));
    108     } else {
    109       // Compute the LU decomposition
    110       var decomp = lup(matrix); // The determinant is the product of the diagonal entries of U (and those of L, but they are all 1)
    111 
    112       var det = decomp.U[0][0];
    113 
    114       for (var _i = 1; _i < rows; _i++) {
    115         det = multiply(det, decomp.U[_i][_i]);
    116       } // The determinant will be multiplied by 1 or -1 depending on the parity of the permutation matrix.
    117       // This can be determined by counting the cycles. This is roughly a linear time algorithm.
    118 
    119 
    120       var evenCycles = 0;
    121       var i = 0;
    122       var visited = [];
    123 
    124       while (true) {
    125         while (visited[i]) {
    126           i++;
    127         }
    128 
    129         if (i >= rows) break;
    130         var j = i;
    131         var cycleLen = 0;
    132 
    133         while (!visited[decomp.p[j]]) {
    134           visited[decomp.p[j]] = true;
    135           j = decomp.p[j];
    136           cycleLen++;
    137         }
    138 
    139         if (cycleLen % 2 === 0) {
    140           evenCycles++;
    141         }
    142       }
    143 
    144       return evenCycles % 2 === 0 ? det : unaryMinus(det);
    145     }
    146   }
    147 });