Random number generators

The random number generators given here originally appeared as posts by Patrik Rak on WoSF and are used here with permission.

The default built-in pseudo random number generator in the ZX Spectrum cycles through 65536 numbers (2^16 sequence of numbers). Patrik proposes two alternatives for better random number generation starting with the XOR shift method (based on George Marsaglia’s derivation), in assembly for the Z80:

Yesterday I had a bit of spare time so I put together a quick test for finding the a,b,c parameters for n=8 and four seeds (I first verified the Marsaglia’s 81 parameters to make sure I got it all right). There are six sets of parameters which provide the desired order of 2^32-1: (1,1,3), (3,6,1), (3,3,2), (5,3,2), (1,7,2), (6,7,1). Of these, the first four seem to do reasonably well on the Diehard test. Here is the Z80 code for the (1,1,3) variant, the result being the 8 bit in the A register:

rnd     ld  hl,0xA280   ; xz -> yw
        ld  de,0xC0DE   ; yw -> zt
        ld  (rnd+1),de  ; x = y, z = w
        ld  a,e         ; w = w ^ ( w << 3 )
        add a,a
        add a,a
        add a,a
        xor e
        ld  e,a
        ld  a,h         ; t = x ^ (x << 1)
        add a,a
        xor h
        ld  d,a
        rra             ; t = t ^ (t >> 1) ^ w
        xor d
        xor e
        ld  h,l         ; y = z
        ld  l,a         ; w = t
        ld  (rnd+4),hl
        ret

In order to return the 8 bit random number to BASIC, do a LD B,0 and LD C, A just before the RET in the end.

Seeing the interest regarding the Xor-Shift random number generator for Z80, I became curious how difficult it would be to implement CMWC (Complimentary-Multiply-With-Carry) RNG for Z80. These generators are based on a prime p of the form a*b^r + 1, carry c < a and a sequence of r x’s < b such that t = a * x(n-r) + c(n-1), c(n) = t div b, x(n) = (b-1) – t mod b.

To make it suitable for Z80, I have chosen base b=2^8. I was also interested how well it could do with fairly small state, so I have chosen r=8 rather than the tempting r=256. From the possible primes, I have chosen the one with a=253, for which the order of the generated sequence is (p-1)/2^5 = 253*2^59 (assuming I got all the math right ), which is almost 2^67, and more than 2^66 or 10^20.

I have tested the C version of the generator with diehard, and it seems to pass all the tests. Here is the corresponding Z80 variant, which returns the 8 bit result in A:

org 40000

    call rnd      ; BASIC driver
    ld   c,a
    ld   b,0
    ret

rnd ld   de,0     ; c,i

    ld   b,0
    ld   c,e
    ld   hl,table
    add  hl,bc

    ld   c,(hl)   ; y = q[i]

    push hl

    ld   a,e      ; i = ( i + 1 ) & 7
    inc  a
    and  7
    ld   e,a

    ld   h,c      ; t = 256 * y
    ld   l,b

    sbc  hl,bc    ; t = 255 * y
    sbc  hl,bc    ; t = 254 * y
    sbc  hl,bc    ; t = 253 * y

    ld   c,d
    add  hl,bc    ; t = 253 * y + c

    ld   d,h      ; c = t / 256

    ld   (rnd+1),de

    ld   a,l      ; x = t % 256
    cpl           ; x = (b-1) - x = -x - 1 = ~x + 1 - 1 = ~x

    pop  hl

    ld   (hl),a   ; q[i] = x

    ret

table
    db   82,97,120,111,102,116,20,12

As you can see, the table access is quite expensive compared to the computation itself. Having the table aligned to 256 byte boundary could be used to optimize this. Also, when generating more than 8 bits at once, most of this could be easily factored out and reused. But I leave these optimizations for anyone interested…

The original posts contain some more tidbits that may be of use including BASIC and ZX Boriel Compiler versions:
http://www.worldofspectrum.org/forums/showthread.php?t=23070

http://www.worldofspectrum.org/forums/showthread.php?t=39632