## Tuesday, August 31, 2010

### 16 bit xorshift rng (now with more period)

So, in Marsaglia's paper, the basic xorshift rng the code is:
`yˆ=(y<<a); y^=(y>>b); return (yˆ=(y<<c));`

I looked at this in my last post, and had a brief look at the (a,b,c) triplets to use. However, this only yields a period length of the word size used. In my case, that's 2^16-1 -- not much. However, further on in the paper, section 3.1, a method is described for longer periods through additional state and promotion:
`t=(xˆ(x<<a)); x=y; return y=(yˆ(y>>c))ˆ(tˆ(t>>b));`

This is listed as 2^64-1 using 32-bit words, but there's no reason we can't use 16-bit words. This reduces the period to 2^32-1, and requires new triplets. I scanned (1,1,1) through till (15,15,15), here's the list:

 (1,1,7) (1,1,12) (1,1,13) (2,5,8) (2,5,13) (2,13,15) (2,15,13) (3,7,6) (5,3,1)* (5,3,8) (5,3,13)* (5,7,4)* (6,3,8)* (7,1,6) (7,1,15) (7,2,1) (8,3,9)* (9,14,5) (11,8,5)* (13,12,3) (14,1,15) (15,10,1)

Values marked with an asterisk score well on the majority of diehard tests. It's probably not a good idea to use the other values, even though they do have a complete period. Notably, all values fail binary matrices tests for 31x31 and 32x32 - I think this is expected, and mentioned in Marsaglia2003.

So, the routine in C:
`uint16_t rnd_xorshift_32() {  static uint16_t x=1,y=1;  uint16_t t=(x^(x<<5));   x=y;   return y=(y^(y>>1))^(t^(t>>3));}`

and in assembly:
`; returns a 16-bit value >0, has a 2^32-1 periodrandom:                push cx        push dx        db 0b8h                 ; t=(x^(x shl a))random_x: dw 1        mov dx,ax        mov cl,5        shl dx,cl        xor ax,dx        mov dx,ax               ; y=(y^(y shr c))^(t^(t shr b))        mov cl,3                ;                  ^^^^^^^^^^^        shr dx,cl        xor ax,dx        push ax                 ; save t^(t shr b)        db 0b8hrandom_y: dw 1        mov cs:[random_x],ax    ; x=y        mov dx,ax               ; y=(y^(y shr c))^(t^(t shr b))        shr dx,1                ;    ^^^^^^^^^^^        xor ax,dx        pop dx        xor ax,dx        mov cs:[random_y],ax        pop dx        pop cx        ret`

Ben Jamin' said...

Which triplet passes the most tests? (5,3,1)?

From memory the asterisked triplets performed equally well.

Russell Borogove said...

Thanks for this - I was looking for a suitable RNG to run on an AVR microcontroller (8-bit, but the register file is generous enough that the compiler can handle wider operations reasonably) and 64-bit xorshift seemed overkill.

To avoid register churn you can use 8-bit sizes and the form mentioned in 3.1:

static uint8_t Q[4]={1,1,1,1};
uint8_t g(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
{
uint8_t t;
t = (Q[0] ^ (Q[0] << i)) ^
(Q[1] ^ (Q[1] >> j)) ^
(Q[2] ^ (Q[2] << k)) ^
(Q[3] ^ (Q[3] << l));
Q[0] = Q[1]; Q[1] = Q[2]; Q[2] = Q[3]; Q[3] = t;
return t;
}

where (i,j,k,l) is one of:

(1,3,1,2)
(1,3,2,3)
(1,4,2,5)
(1,6,2,3)
(1,6,7,3)
(1,7,1,3)
(1,7,3,5)
(2,1,1,1)
(3,1,1,1)
(3,5,2,5)
(3,5,4,5)
(3,5,7,5)
(3,6,2,1)
(3,6,5,1)
(4,5,1,6)
(5,3,2,3)
(6,2,1,5)
(6,3,3,1)
(6,3,7,1)
(6,3,7,4)
(7,1,3,5)
(7,3,4,3)
(7,5,3,2)
(7,7,2,1)
(7,7,4,1)
(7,7,6,1)

hopper said...

Is it necessary to start with x = y = 1 or can I use any non-zero values?

Thanks

M

hopper said...

I'm also curious how to test 16 bit values with Diehard or other test suite. As far as I knew, they mostly expect 32 bit values?

Hi M,

For initialization you'll need to ensure `x != 0 || y != 0`.

I used the Diehard test suite at the time, but nowadays you can use `PractRand stdin16`.