extra/math/matrices/matrices.factor

   1 ! Copyright (C) 2005, 2008 Slava Pestov.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays kernel sequences math math.functions
   4 math.vectors math.order ;
   5 IN: math.matrices
   6
   7 ! Matrices
   8 : zero-matrix ( m n -- matrix )
   9     [ nip 0 <array> ] curry map ;
  10
  11 : identity-matrix ( n -- matrix )
  12     #! Make a nxn identity matrix.
  13     dup [ [ = 1 0 ? ] with map ] curry map ;
  14
  15 ! Matrix operations
  16 : mneg ( m -- m ) [ vneg ] map ;
  17
  18 : n*m ( n m -- m ) [ n*v ] with map ;
  19 : m*n ( m n -- m ) [ v*n ] curry map ;
  20 : n/m ( n m -- m ) [ n/v ] with map ;
  21 : m/n ( m n -- m ) [ v/n ] curry map ;
  22
  23 : m+   ( m m -- m ) [ v+ ] 2map ;
  24 : m-   ( m m -- m ) [ v- ] 2map ;
  25 : m*   ( m m -- m ) [ v* ] 2map ;
  26 : m/   ( m m -- m ) [ v/ ] 2map ;
  27
  28 : v.m ( v m -- v ) flip [ v. ] with map ;
  29 : m.v ( m v -- v ) [ v. ] curry map ;
  30 : m.  ( m m -- m ) flip [ swap m.v ] curry map ;
  31
  32 : mmin ( m -- n ) >r 1/0. r> [ [ min ] each ] each ;
  33 : mmax ( m -- n ) >r -1/0. r> [ [ max ] each ] each ;
  34 : mnorm ( m -- n ) dup mmax abs m/n ;
  35
  36 <PRIVATE
  37
  38 : x ( seq -- elt ) first ; inline
  39 : y ( seq -- elt ) second ; inline
  40 : z ( seq -- elt ) third ; inline
  41
  42 : i ( seq1 seq2 -- n ) [ [ y ] [ z ] bi* * ] [ [ z ] [ y ] bi* * ] 2bi - ;
  43 : j ( seq1 seq2 -- n ) [ [ z ] [ x ] bi* * ] [ [ x ] [ z ] bi* * ] 2bi - ;
  44 : k ( seq1 seq2 -- n ) [ [ y ] [ x ] bi* * ] [ [ x ] [ y ] bi* * ] 2bi - ;
  45
  46 PRIVATE>
  47
  48 : cross ( vec1 vec2 -- vec3 ) [ i ] [ j ] [ k ] 2tri 3array ;
  49
  50 : proj ( v u -- w )
  51     [ [ v. ] [ norm-sq ] bi / ] keep n*v ;
  52
  53 : (gram-schmidt) ( v seq -- newseq )
  54     [ dupd proj v- ] each ;
  55
  56 : gram-schmidt ( seq -- orthogonal )
  57     V{ } clone [ over (gram-schmidt) over push ] reduce ;
  58
  59 : norm-gram-schmidt ( seq -- orthonormal )
  60     gram-schmidt [ normalize ] map ;