2005-04-13 20:44:06 -04:00
|
|
|
! Copyright (C) 2005 Slava Pestov.
|
|
|
|
! See http://factor.sf.net/license.txt for BSD license.
|
|
|
|
IN: matrices
|
2005-05-03 04:40:13 -04:00
|
|
|
USING: errors generic kernel lists math namespaces sequences
|
|
|
|
vectors ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-23 00:25:52 -04:00
|
|
|
: n*v ( n vec -- vec ) [ * ] map-with ;
|
2005-06-23 22:35:41 -04:00
|
|
|
: v*n ( vec n -- vec ) swap n*v ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
|
|
|
! Vector operations
|
2005-05-14 17:18:45 -04:00
|
|
|
: v+ ( v v -- v ) [ + ] 2map ;
|
|
|
|
: v- ( v v -- v ) [ - ] 2map ;
|
|
|
|
: v* ( v v -- v ) [ * ] 2map ;
|
2005-06-22 02:32:17 -04:00
|
|
|
: v/ ( v v -- v ) [ / ] 2map ;
|
2005-05-19 15:16:25 -04:00
|
|
|
: v** ( v v -- v ) [ conjugate * ] 2map ;
|
2005-06-22 02:32:17 -04:00
|
|
|
: vmax ( v v -- v ) [ max ] 2map ;
|
|
|
|
: vmin ( v v -- v ) [ min ] 2map ;
|
|
|
|
: vneg ( v -- v ) [ neg ] map ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-06-25 20:39:53 -04:00
|
|
|
: sum ( v -- n ) 0 [ + ] reduce ;
|
|
|
|
: product 1 [ * ] reduce ;
|
|
|
|
|
|
|
|
: set-axis ( x y axis -- v )
|
|
|
|
2dup v* >r >r drop dup r> v* v- r> v+ ;
|
2005-05-28 20:52:23 -04:00
|
|
|
|
2005-05-19 15:16:25 -04:00
|
|
|
! Later, this will fixed when 2each works properly
|
2005-05-23 00:25:52 -04:00
|
|
|
! : v. ( v v -- x ) 0 swap [ conjugate * + ] 2each ;
|
2005-05-28 20:52:23 -04:00
|
|
|
: v. ( v v -- x ) v** sum ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-22 22:08:46 -04:00
|
|
|
: cross-trace ( v1 v2 i1 i2 -- v1 v2 n )
|
|
|
|
pick nth >r pick nth r> * ;
|
|
|
|
|
|
|
|
: cross-minor ( v1 v2 i1 i2 -- n )
|
|
|
|
[ cross-trace -rot ] 2keep swap cross-trace 2nip - ;
|
2005-05-21 02:28:23 -04:00
|
|
|
|
|
|
|
: cross ( { x1 y1 z1 } { x2 y2 z2 } -- { z1 z2 z3 } )
|
|
|
|
#! Cross product of two 3-dimensional vectors.
|
2005-05-22 22:08:46 -04:00
|
|
|
3 <vector>
|
|
|
|
[ >r 2dup 1 2 cross-minor 0 r> set-nth ] keep
|
|
|
|
[ >r 2dup 2 0 cross-minor 1 r> set-nth ] keep
|
|
|
|
[ >r 2dup 0 1 cross-minor 2 r> set-nth ] keep
|
|
|
|
2nip ;
|
2005-05-20 23:52:31 -04:00
|
|
|
|
2005-04-30 02:01:04 -04:00
|
|
|
! Matrices
|
2005-05-03 04:40:13 -04:00
|
|
|
! The major dimension is the number of elements per row.
|
|
|
|
TUPLE: matrix rows cols sequence ;
|
|
|
|
: >matrix<
|
|
|
|
[ matrix-rows ] keep
|
|
|
|
[ matrix-cols ] keep
|
|
|
|
matrix-sequence ;
|
|
|
|
|
2005-04-13 20:44:06 -04:00
|
|
|
M: matrix clone ( matrix -- matrix )
|
|
|
|
clone-tuple
|
|
|
|
dup matrix-sequence clone over set-matrix-sequence ;
|
|
|
|
|
2005-05-23 00:25:52 -04:00
|
|
|
: matrix@ ( row col matrix -- n ) matrix-cols rot * + ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
|
|
|
: matrix-get ( row col matrix -- elt )
|
|
|
|
[ matrix@ ] keep matrix-sequence nth ;
|
|
|
|
|
|
|
|
: matrix-set ( elt row col matrix -- )
|
|
|
|
[ matrix@ ] keep matrix-sequence set-nth ;
|
|
|
|
|
|
|
|
: <zero-matrix> ( rows cols -- matrix )
|
|
|
|
2dup * zero-vector <matrix> ;
|
|
|
|
|
2005-05-28 20:52:23 -04:00
|
|
|
: <row-matrix> ( vector -- matrix )
|
2005-04-30 02:01:04 -04:00
|
|
|
#! Turn a vector into a matrix of one row.
|
|
|
|
[ 1 swap length ] keep <matrix> ;
|
|
|
|
|
2005-05-28 20:52:23 -04:00
|
|
|
: <col-matrix> ( vector -- matrix )
|
2005-04-30 02:01:04 -04:00
|
|
|
#! Turn a vector into a matrix of one column.
|
|
|
|
[ length 1 ] keep <matrix> ;
|
|
|
|
|
2005-05-05 15:31:57 -04:00
|
|
|
: make-matrix ( rows cols quot -- matrix | quot: i j -- elt )
|
|
|
|
-rot [
|
|
|
|
[ [ [ rot call , ] 3keep ] 2repeat ] make-vector nip
|
|
|
|
] 2keep rot <matrix> ; inline
|
2005-04-13 20:44:06 -04:00
|
|
|
|
|
|
|
: <identity-matrix> ( n -- matrix )
|
|
|
|
#! Make a nxn identity matrix.
|
|
|
|
dup [ = 1 0 ? ] make-matrix ;
|
|
|
|
|
|
|
|
: transpose ( matrix -- matrix )
|
|
|
|
dup matrix-cols over matrix-rows [
|
2005-05-22 22:08:46 -04:00
|
|
|
swap pick matrix-get
|
2005-04-13 20:44:06 -04:00
|
|
|
] make-matrix nip ;
|
|
|
|
|
|
|
|
! Sequence of elements in a row of a matrix.
|
|
|
|
TUPLE: row index matrix ;
|
|
|
|
: >row< dup row-index swap row-matrix ;
|
|
|
|
M: row length row-matrix matrix-cols ;
|
2005-05-28 20:52:23 -04:00
|
|
|
M: row nth ( n row -- n ) >row< swapd matrix-get ;
|
2005-04-30 02:01:04 -04:00
|
|
|
M: row thaw >vector ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
|
|
|
! Sequence of elements in a column of a matrix.
|
|
|
|
TUPLE: col index matrix ;
|
|
|
|
: >col< dup col-index swap col-matrix ;
|
|
|
|
M: col length col-matrix matrix-rows ;
|
2005-05-28 20:52:23 -04:00
|
|
|
M: col nth ( n column -- n ) >col< matrix-get ;
|
2005-04-30 02:01:04 -04:00
|
|
|
M: col thaw >vector ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-28 20:52:23 -04:00
|
|
|
! Sequence of elements on a diagonal. Positive indices are above
|
|
|
|
! and negative indices are below the main diagonal. Only for
|
|
|
|
! square matrices.
|
|
|
|
TUPLE: diagonal index matrix ;
|
|
|
|
: >diagonal< dup diagonal-index swap diagonal-matrix ;
|
|
|
|
M: diagonal length ( daig -- n )
|
|
|
|
>diagonal< matrix-rows swap abs - ;
|
|
|
|
M: diagonal nth ( n diag -- n )
|
|
|
|
>diagonal< >r [ neg 0 max over + ] keep 0 max rot + r>
|
|
|
|
matrix-get ;
|
|
|
|
|
|
|
|
: trace ( matrix -- tr )
|
|
|
|
#! Product of diagonal elements.
|
|
|
|
0 swap <diagonal> product ;
|
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: +check ( matrix matrix -- )
|
2005-04-13 20:44:06 -04:00
|
|
|
#! Check if the two matrices have dimensions compatible
|
|
|
|
#! for being added or subtracted.
|
|
|
|
over matrix-rows over matrix-rows = >r
|
2005-05-03 04:40:13 -04:00
|
|
|
swap matrix-cols swap matrix-cols = r> and [
|
|
|
|
"Matrix dimensions do not equal" throw
|
2005-04-13 20:44:06 -04:00
|
|
|
] unless ;
|
|
|
|
|
2005-05-05 15:31:57 -04:00
|
|
|
: element-wise ( m m -- rows cols v v )
|
2005-05-03 04:40:13 -04:00
|
|
|
2dup +check >r >matrix< r> matrix-sequence ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
! Matrix operations
|
|
|
|
: m+ ( m m -- m ) element-wise v+ <matrix> ;
|
|
|
|
: m- ( m m -- m ) element-wise v- <matrix> ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: m* ( m m -- m )
|
|
|
|
#! Multiply two matrices element-wise. This is NOT matrix
|
|
|
|
#! multiplication in the usual mathematical sense. For that,
|
|
|
|
#! see the m. word.
|
|
|
|
element-wise v* <matrix> ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: *check ( matrix matrix -- )
|
|
|
|
swap matrix-cols swap matrix-rows = [
|
2005-04-14 19:37:13 -04:00
|
|
|
"Matrix dimensions inappropriate for composition" throw
|
2005-04-13 20:44:06 -04:00
|
|
|
] unless ;
|
|
|
|
|
|
|
|
: *dimensions ( m m -- rows cols )
|
|
|
|
swap matrix-rows swap matrix-cols ;
|
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: m. ( m1 m2 -- m )
|
|
|
|
#! Composition of two matrices.
|
|
|
|
2dup *check 2dup *dimensions [
|
|
|
|
( m1 m2 row col -- m1 m2 )
|
2005-05-23 00:25:52 -04:00
|
|
|
pick <col> >r pick <row> r> v.
|
2005-04-13 20:44:06 -04:00
|
|
|
] make-matrix 2nip ;
|
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: n*m ( n m -- m )
|
|
|
|
#! Multiply a matrix by a scalar.
|
|
|
|
>matrix< >r rot r> n*v <matrix> ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: m.v ( m v -- v )
|
|
|
|
#! Multiply a matrix by a column vector.
|
2005-05-28 20:52:23 -04:00
|
|
|
<col-matrix> m. matrix-sequence ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
2005-05-03 04:40:13 -04:00
|
|
|
: v.m ( v m -- v )
|
|
|
|
#! Multiply a row vector by a matrix.
|
2005-05-28 20:52:23 -04:00
|
|
|
>r <row-matrix> r> m. matrix-sequence ;
|
2005-04-13 20:44:06 -04:00
|
|
|
|
|
|
|
: row-list ( matrix -- list )
|
|
|
|
#! A list of lists, where each sublist is a row of the
|
|
|
|
#! matrix.
|
2005-05-03 04:40:13 -04:00
|
|
|
dup matrix-rows [ swap <row> >list ] project-with ;
|