extra/random/lagged-fibonacci/lagged-fibonacci.factor

   1 ! Copyright (C) 2009 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors alien.c-types fry kernel literals locals math
   4 random sequences specialized-arrays namespaces ;
   5 SPECIALIZED-ARRAY: double
   6 IN: random.lagged-fibonacci
   7
   8 TUPLE: lagged-fibonacci u pt0 pt1 ;
   9
  10 <PRIVATE
  11
  12 CONSTANT: p-r 1278
  13 CONSTANT: q-r 417
  14
  15 CONSTANT: lagged-fibonacci 899999963
  16 CONSTANT: lagged-fibonacci-max-seed 900000000
  17 CONSTANT: lagged-fibonacci-sig-bits 24
  18
  19 : normalize-seed ( seed -- seed' )
  20     abs lagged-fibonacci-max-seed mod ;
  21
  22 : adjust-ptr ( ptr -- ptr' )
  23     1 - dup 0 < [ drop p-r ] when ;
  24
  25 PRIVATE>
  26
  27 M:: lagged-fibonacci seed-random ( lagged-fibonacci seed! -- lagged-fibonacci )
  28     seed normalize-seed seed!
  29     seed 30082 /i :> ij
  30     seed 30082 ij * - :> kl
  31     ij 177 /i 177 mod 2 + :> i!
  32     ij 177 mod 2 + :> j!
  33     kl 169 /i 178 mod 1 + :> k!
  34     kl 169 mod :> l!
  35
  36     lagged-fibonacci u>> [
  37         drop
  38         0.0 :> s!
  39         0.5 :> t!
  40         0.0 :> m!
  41         lagged-fibonacci-sig-bits [
  42             i j * 179 mod k * 179 mod m!
  43             j i!
  44             k j!
  45             m k!
  46             53 l * 1 + 169 mod l!
  47             l m * 64 mod 31 > [ s t + s! ] when
  48             t 0.5 * t!
  49         ] times
  50         s
  51     ] change-each
  52     lagged-fibonacci p-r >>pt0
  53         q-r >>pt1 ;
  54
  55 : <lagged-fibonacci> ( seed -- lagged-fibonacci )
  56     lagged-fibonacci new
  57         p-r 1 + <double-array> >>u
  58         swap seed-random ;
  59
  60 GENERIC: random-float* ( tuple -- r )
  61
  62 : random-float ( -- n ) random-generator get random-float* ;
  63
  64 M:: lagged-fibonacci random-float* ( lagged-fibonacci -- x )
  65     lagged-fibonacci [ pt0>> ] [ u>> ] bi nth
  66     lagged-fibonacci [ pt1>> ] [ u>> ] bi nth - :> uni!
  67     uni 0.0 < [ uni 1.0 + uni! ] when
  68     uni lagged-fibonacci [ pt0>> ] [ u>> ] bi set-nth
  69     lagged-fibonacci [ adjust-ptr ] change-pt0 drop
  70     lagged-fibonacci [ adjust-ptr ] change-pt1 drop
  71     uni ; inline