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: float-arrays kernel math math.functions math.vectors
   4 sequences sequences.private prettyprint words tools.time hints ;
   5 IN: benchmark.spectral-norm
   6
   7 : fast-truncate >fixnum >float ; inline
   8
   9 : eval-A ( i j -- n )
  10     [ >float ] bi@
  11     dupd + dup 1+ * 2 /f fast-truncate + 1+
  12     recip ; inline
  13
  14 : (eval-A-times-u) ( u i j -- x )
  15     tuck eval-A >r swap nth-unsafe r> * ; inline
  16
  17 : eval-A-times-u ( n u -- seq )
  18     over [
  19         pick 0.0 [
  20             swap >r >r 2dup r> (eval-A-times-u) r> +
  21         ] reduce nip
  22     ] F{ } map-as { float-array } declare 2nip ; inline
  23
  24 : (eval-At-times-u) ( u i j -- x )
  25     tuck swap eval-A >r swap nth-unsafe r> * ; inline
  26
  27 : eval-At-times-u ( n u -- seq )
  28     over [
  29         pick 0.0 [
  30             swap >r >r 2dup r> (eval-At-times-u) r> +
  31         ] reduce nip
  32     ] F{ } map-as { float-array } declare 2nip ; inline
  33
  34 : eval-AtA-times-u ( n u -- seq )
  35     dupd eval-A-times-u eval-At-times-u ; inline
  36
  37 : u/v ( n -- u v )
  38     dup 1.0 <float-array> dup
  39     10 [
  40         drop
  41         dupd eval-AtA-times-u
  42         2dup eval-AtA-times-u
  43         swap
  44     ] times
  45     rot drop ; inline
  46
  47 : spectral-norm ( n -- norm )
  48     u/v [ v. ] keep norm-sq /f sqrt ;
  49
  50 HINTS: spectral-norm fixnum ;
  51
  52 : spectral-norm-main ( -- )
  53     2000 spectral-norm . ;
  54
  55 MAIN: spectral-norm-main