unmaintained/adsoda/solution2/solution2.factor

   1 USING: kernel
   2 sequences
   3 namespaces
   4
   5 math
   6 math.vectors
   7 math.matrices
   8 ;
   9 IN: adsoda.solution2
  10
  11 ! -------------------
  12 ! correctif solution
  13 ! ---------------
  14 SYMBOL: matrix
  15 : MIN-VAL-adsoda ( -- x ) 0.00000001
  16 ! 0.000000000001
  17 ;
  18
  19 : zero? ( x -- ? )
  20     abs MIN-VAL-adsoda <
  21 ;
  22
  23 ! [ number>string string>number ] map
  24
  25 : with-matrix ( matrix quot -- )
  26     [ swap matrix set call matrix get ] with-scope ; inline
  27
  28 : nth-row ( row# -- seq ) matrix get nth ;
  29
  30 : change-row ( row# quot -- seq ) ! row# quot -- | quot: seq -- seq )
  31     matrix get swap change-nth ; inline
  32
  33 : exchange-rows ( row# row# -- ) matrix get exchange ;
  34
  35 : rows ( -- n ) matrix get length ;
  36
  37 : cols ( -- n ) 0 nth-row length ;
  38
  39 : skip ( i seq quot -- n )
  40     over [ find-from drop ] dip length or ; inline
  41
  42 : first-col ( row# -- n )
  43     ! First non-zero column
  44     0 swap nth-row [ zero? not ] skip ;
  45
  46 : clear-scale ( col# pivot-row i-row -- n )
  47     [ over ] dip nth dup zero? [
  48         3drop 0
  49     ] [
  50         [ nth dup zero? ] dip swap [
  51             2drop 0
  52         ] [
  53             swap / neg
  54         ] if
  55     ] if ;
  56
  57 : (clear-col) ( col# pivot-row i -- )
  58     [ [ clear-scale ] 2keep [ n*v ] dip v+ ] change-row ;
  59
  60 : rows-from ( row# -- slice )
  61     rows dup <slice> ;
  62
  63 : clear-col ( col# row# rows -- )
  64     [ nth-row ] dip [ [ 2dup ] dip (clear-col) ] each 2drop ;
  65
  66 : do-row ( exchange-with row# -- )
  67     [ exchange-rows ] keep
  68     [ first-col ] keep
  69     dup 1 + rows-from clear-col ;
  70
  71 : find-row ( row# quot -- i elt )
  72     [ rows-from ] dip find ; inline
  73
  74 : pivot-row ( col# row# -- n )
  75     [ dupd nth-row nth zero? not ] find-row 2nip ;
  76
  77 : (echelon) ( col# row# -- )
  78     over cols < over rows < and [
  79         2dup pivot-row [ over do-row 1 + ] when*
  80         [ 1 + ] dip (echelon)
  81     ] [
  82         2drop
  83     ] if ;
  84
  85 : echelon ( matrix -- matrix' )
  86     [ 0 0 (echelon) ] with-matrix ;
  87
  88 : nonzero-rows ( matrix -- matrix' )
  89     [ [ zero? ] all? ] reject ;
  90
  91 : null/rank ( matrix -- null rank )
  92     echelon dup length swap nonzero-rows length [ - ] keep ;
  93
  94 : leading ( seq -- n elt ) [ zero? not ] find ;
  95
  96 : reduced ( matrix' -- matrix'' )
  97     [
  98         rows <reversed> [
  99             dup nth-row leading drop
 100             dup [ swap dup clear-col ] [ 2drop ] if
 101         ] each
 102     ] with-matrix ;
 103
 104 : basis-vector ( row col# -- )
 105     [ clone ] dip
 106     [ swap nth neg recip ] 2keep
 107     [ 0 spin set-nth ] 2keep
 108     [ n*v ] dip
 109     matrix get set-nth ;
 110
 111 : nullspace ( matrix -- seq )
 112     echelon reduced dup empty? [
 113         dup first length identity-matrix [
 114             [
 115                 dup leading drop
 116                 dup [ basis-vector ] [ 2drop ] if
 117             ] each
 118         ] with-matrix flip nonzero-rows
 119     ] unless ;
 120
 121 : 1-pivots ( matrix -- matrix )
 122     [ dup leading nip [ recip v*n ] when* ] map ;
 123
 124 : solution ( matrix -- matrix )
 125     echelon nonzero-rows reduced 1-pivots ;
 126