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