--- /dev/null
+! Copyright (C) 2013 Doug Coleman.
+! See http://factorcode.org/license.txt for BSD license.
+USING: tools.test kernel ;
+IN: math.matrices.laplace
+
+{ -2 } [ { { 1 2 } { 3 4 } } determinant ] unit-test
+
+{ 0 } [
+ { { 1 2 3 } { 4 5 6 } { 7 8 9 } } determinant
+] unit-test
+
+{ -47860032 } [
+ {
+ { 40 39 38 37 }
+ { 1 1 1 831 }
+ { 22 22 1110 299 }
+ { 13 14 15 17 }
+ } determinant
+] unit-test
--- /dev/null
+! Copyright (C) 2013 Doug Coleman.
+! See http://factorcode.org/license.txt for BSD license.
+USING: accessors arrays fry kernel locals math math.matrices
+math.vectors sequences sequences.private ;
+IN: math.matrices.laplace
+
+<PRIVATE
+
+: 2x2-determinant ( matrix -- x )
+ first2 [ first2 ] bi@ -rot [ * ] 2bi@ - ;
+
+! using a virtual "missing element" sequence for performance
+TUPLE: missing seq i ;
+C: <missing> missing
+M: missing nth-unsafe
+ [ i>> dupd >= [ 1 + ] when ] [ seq>> nth-unsafe ] bi ;
+M: missing length seq>> length 1 - ;
+INSTANCE: missing immutable-sequence
+
+: first-sub-matrix ( matrix -- first-row seq )
+ [ unclip-slice swap ] [ length iota ] bi
+ [ '[ _ <missing> ] map ] with map ;
+
+:: laplace-expansion ( row matrix -- x )
+ matrix length 2 =
+ [ matrix 2x2-determinant ] [
+ matrix first-sub-matrix ! cheat, always expand on first row
+ [ row swap laplace-expansion ] map
+ v* [ odd? [ neg ] when ] map-index sum
+ ] if ;
+
+ERROR: not-a-square-matrix matrix ;
+
+: check-square-matrix ( matrix -- matrix )
+ dup square-matrix? [ not-a-square-matrix ] unless ; inline
+
+PRIVATE>
+
+: determinant ( matrix -- x )
+ check-square-matrix 0 swap laplace-expansion ;