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