extra/adsoda/solution2/solution2.factor

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