extra/random/cmwc/cmwc.factor

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