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