extra/random/pcg/pcg.factor

   1 ! Copyright (C) 2022 Alexander Ilin.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors alien.c-types classes.struct kernel locals math
   4 math.bitwise random sets system ;
   5 IN: random.pcg
   6
   7 ! https://www.pcg-random.org/
   8 ! https://github.com/tkaitchuck/Mwc256XXA64/blob/main/impl/src/gen32.rs
   9 ! https://github.com/tkaitchuck/Mwc256XXA64/blob/main/impl/src/gen64.rs
  10
  11 ! Since we are only returning 32 bits per step of this 64-bit PRNG, the rest are
  12 ! saved in the rem field.
  13 TUPLE: Mwc256XXA64 x1 x2 x3 c rem ;
  14 TUPLE: Mwc128XXA32 x1 x2 x3 c ;
  15
  16 <PRIVATE
  17
  18 CONSTANT: MULTIPLIER 0xfeb3_4465_7c0a_f413
  19 CONSTANT: MULTIPLIER-32 3487286589
  20
  21 : big>d/d ( n -- low high )
  22     dup -64 shift [ 64 bits ] bi@ ; inline
  23
  24 : multiply-32 ( n -- low high )
  25     MULTIPLIER-32 * d>w/w ; inline
  26
  27 : multiply ( n -- low high )
  28     MULTIPLIER * big>d/d ; inline
  29
  30 : permute-32 ( high x1 x2 x3 -- n )
  31     [ bitxor ] 2bi@ w+ ; inline
  32
  33 : permute ( high x1 x2 x3 -- n )
  34     [ bitxor ] 2bi@ W+ ; inline
  35
  36 :: rot-state ( obj x1 c -- struct' )
  37     obj
  38         obj x2>> >>x3
  39         obj x1>> >>x2
  40         x1 >>x1
  41         c >>c ; inline
  42
  43 : update-state-32 ( obj low high -- )
  44     [ over c>> + d>w/w ] dip w+ rot-state drop ; inline
  45
  46 : update-state ( obj low high -- )
  47     [ over c>> + big>d/d ] dip W+ rot-state drop ; inline
  48
  49 : next-u32 ( obj -- n )
  50     dup x3>> multiply-32 [ pick [ x1>> ] [ x2>> ] [ x3>> ] tri permute-32 ] keep
  51     swap [ update-state-32 ] dip ; inline
  52
  53 : next-u64 ( obj -- n )
  54     dup x3>> multiply [ pick [ x1>> ] [ x2>> ] [ x3>> ] tri permute ] keep
  55     swap [ update-state ] dip ; inline
  56
  57 PRIVATE>
  58
  59 : <Mwc128XXA32> ( key1 key2 -- obj )
  60     0xcafef00d 0xd15ea5e5 Mwc128XXA32 boa
  61     6 [ dup next-u32 drop ] times ;
  62
  63 : <Mwc256XXA64> ( key1 key2 -- obj )
  64     0xcafef00dd15ea5e5 0x14057B7EF767814F f Mwc256XXA64 boa
  65     6 [ dup next-u64 drop ] times ;
  66
  67 : <pcg> ( key1 key2 -- obj )
  68     cpu { x86.32 ppc.32 arm.32 } in? [ <Mwc128XXA32> ] [ <Mwc256XXA64> ] if ;
  69
  70 M: Mwc128XXA32 random-32*
  71     next-u32 ;
  72
  73 M: Mwc256XXA64 random-32*
  74     dup '[ [ f ] [ _ next-u64 d>w/w ] if* ] change-rem drop ;
  75
  76 ! USING: random random.pcg ;
  77 ! gc 0 0 random.pcg:<Mwc256XXA64> [ 10,000,000 [ dup random-32* drop ] times ] time drop
  78 ! gc 0 0 random.pcg:<Mwc128XXA32> [ 10,000,000 [ dup random-32* drop ] times ] time drop