extra/random/cmwc/cmwc.factor

   1 ! Copyright (C) 2009 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors arrays fry kernel locals math math.bitwise
   4 random sequences ;
   5 IN: random.cmwc
   6
   7 ! Multiply-with-carry RNG
   8
   9 TUPLE: cmwc Q a b c i r mod ;
  10
  11 TUPLE: cmwc-seed Q c ;
  12
  13 : <cmwc> ( length a b c -- cmwc )
  14     cmwc new
  15         swap >>c
  16         swap >>b
  17         swap >>a
  18         swap [ 1 - >>i ] [ 0 <array> >>Q ] bi
  19         dup b>> 1 - >>r
  20         dup Q>> length 1 - >>mod ;
  21
  22 : <cmwc-seed> ( Q c -- cmwc-seed )
  23     cmwc-seed new
  24         swap >>c
  25         swap >>Q ; inline
  26
  27 M: cmwc seed-random
  28     [ Q>> >>Q ]
  29     [ Q>> length 1 - >>i ]
  30     [ c>> >>c ] tri ;
  31
  32 M:: cmwc random-32* ( cmwc -- n )
  33     cmwc dup mod>> '[ 1 + _ bitand ] change-i
  34     [ a>> ]
  35     [ [ i>> ] [ Q>> ] bi nth * ]
  36     [ c>> + ] tri :> t!
  37
  38     t -32 shift cmwc (>>c)
  39
  40     t cmwc [ b>> bitand ] [ c>> + ] bi 64 bits t!
  41     t cmwc r>> > [
  42         cmwc [ 1 + ] change-c drop
  43         t cmwc b>> - 64 bits t!
  44     ] when
  45
  46     cmwc [ r>> t - 32 bits dup ] [ i>> ] [ Q>> ] tri set-nth ;
  47
  48 : cmwc-4096 ( -- cmwc )
  49     4096
  50     [ 18782 4294967295 362436 <cmwc> ]
  51     [
  52         '[ [ random-32 ] replicate ] with-system-random
  53         362436 <cmwc-seed> seed-random
  54     ] bi ;