extra/benchmark/spectral-norm/spectral-norm.factor

   1 ! Factor port of
   2 ! http://shootout.alioth.debian.org/gp4/benchmark.php?test=spectralnorm&lang=all
   3 USING: specialized-arrays.double kernel math math.functions
   4 math.vectors sequences sequences.private prettyprint words hints
   5 locals ;
   6 IN: benchmark.spectral-norm
   7
   8 :: inner-loop ( u n quot -- seq )
   9     n [| i |
  10         n 0.0 [| j |
  11             u i j quot call +
  12         ] reduce
  13     ] double-array{ } map-as ; inline
  14
  15 : eval-A ( i j -- n )
  16     [ >float ] bi@
  17     [ drop ] [ + [ ] [ 1 + ] bi * 0.5 * ] 2bi
  18     + 1 + recip ; inline
  19
  20 : (eval-A-times-u) ( u i j -- x )
  21     tuck [ swap nth-unsafe ] [ eval-A ] 2bi* * ; inline
  22
  23 : eval-A-times-u ( n u -- seq )
  24     [ (eval-A-times-u) ] inner-loop ; inline
  25
  26 : (eval-At-times-u) ( u i j -- x )
  27     tuck [ swap nth-unsafe ] [ swap eval-A ] 2bi* * ; inline
  28
  29 : eval-At-times-u ( u n -- seq )
  30     [ (eval-At-times-u) ] inner-loop ; inline
  31
  32 : eval-AtA-times-u ( u n -- seq )
  33     [ eval-A-times-u ] [ eval-At-times-u ] bi ; inline
  34
  35 : ones ( n -- seq ) [ 1.0 ] double-array{ } replicate-as ; inline
  36
  37 :: u/v ( n -- u v )
  38     n ones dup
  39     10 [
  40         drop
  41         n eval-AtA-times-u
  42         [ n eval-AtA-times-u ] keep
  43     ] times ; inline
  44
  45 : spectral-norm ( n -- norm )
  46     u/v [ v. ] [ norm-sq ] bi /f sqrt ;
  47
  48 HINTS: spectral-norm fixnum ;
  49
  50 : spectral-norm-main ( -- )
  51     2000 spectral-norm . ;
  52
  53 MAIN: spectral-norm-main