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 period
random:        
        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 0b8h
random_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