extra/math/matrices/laplace/laplace.factor

   1 ! Copyright (C) 2013 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors arrays fry kernel locals math math.matrices
   4 math.vectors sequences sequences.private ;
   5 IN: math.matrices.laplace
   6
   7 <PRIVATE
   8
   9 : 2x2-determinant ( matrix -- x )
  10     first2 [ first2 ] bi@ -rot [ * ] 2bi@ - ;
  11
  12 ! using a virtual "missing element" sequence for performance
  13 TUPLE: missing seq i ;
  14 C: <missing> missing
  15 M: missing nth-unsafe
  16     [ i>> dupd >= [ 1 + ] when ] [ seq>> nth-unsafe ] bi ;
  17 M: missing length seq>> length 1 - ;
  18 INSTANCE: missing immutable-sequence
  19
  20 : first-sub-matrix ( matrix -- first-row seq )
  21     [ unclip-slice swap ] [ length <iota> ] bi
  22     [ '[ _ <missing> ] map ] with map ;
  23
  24 :: laplace-expansion ( row matrix -- x )
  25     matrix length 2 =
  26     [ matrix 2x2-determinant ] [
  27         matrix first-sub-matrix ! cheat, always expand on first row
  28         [ row swap laplace-expansion ] map
  29         v* [ odd? [ neg ] when ] map-index sum
  30     ] if ;
  31
  32 ERROR: not-a-square-matrix matrix ;
  33
  34 : check-square-matrix ( matrix -- matrix )
  35     dup square-matrix? [ not-a-square-matrix ] unless ; inline
  36
  37 PRIVATE>
  38
  39 : determinant ( matrix -- x )
  40     check-square-matrix 0 swap laplace-expansion ;