BBC BASIC Programmers' Reference

High Quality Random Number Generation

by Richard Russell, March 2010

Listed below are BBC BASIC implementations of two high performance Pseudo Random Number Generators: the WELL512A and WELL1024A algorithms developed in 2006 (see this Wikipedia article). They have state sizes of 512 bits and 1024 bits, and therefore sequence lengths of 2^512-1 and 2^1024-1, respectively. They both return 32-bit signed integer values, but these can easily be converted to other formats, for example using the code snippets listed at Notes on the use of RND.

These Well Equidistributed Long-period Linear algorithms have long sequences, excellent statistical properties and are computationally simpler than some alternative methods. They are ideally suited to applications such as Monte Carlo simulations.

In the listings below, BBC BASIC's RND function is used to provide the initialisation (seed). This might be adequate for some uses, but for more critical applications a seed with more entropy would be desirable. The 'state' information is contained in the array State%() and the variable State_i% which are passed by reference to the generator function. To 'replay' a sequence it is only necessary to save and restore these values.

        REM BBC BASIC version of the WELL512A Pseudo Random Number Generator
        REM By Richard Russell, http://www.rtrussell.co.uk/, v1.00, 30-Mar-2010
        REM Based on the original C code by François Panneton & Pierre L'Écuyer
        REM (Université de Montréal) and Makoto Matsumoto (Hiroshima University)
        REM http://en.wikipedia.org/wiki/Well_Equidistributed_Long-period_Linear
 
        DIM State%(15)
        State_i% = 0
 
        REM Initialisation:
        FOR j% = 0 TO 15
          State%(j%) = RND
        NEXT j%
 
        FOR j% = 1 TO 32
          PRINT RIGHT$("0000000"+STR$~(FN_WELLRNG512a(State_i%, State%())), 8)
        NEXT j%
        END
 
        DEF FN_WELLRNG512a(RETURN I%, S%())
        LOCAL X%, Y%, Z%
        X% = S%(I%) EOR S%(I%)<<16 EOR S%(I%+13AND15) EOR S%(I%+13AND15)<<15
        Y% = S%(I%+9AND15) EOR S%(I%+9AND15)>>>11
        S%(I%) = X% EOR Y%
        Z% = S%(I%) EOR S%(I%)<<5 AND &DA442D24
        I% = I%+15AND15
        S%(I%) = S%(I%) EOR S%(I%)<<2 EOR X% EOR X%<<18 EOR Y%<<28 EOR Z%
        = S%(I%)

        REM BBC BASIC version of the WELL1024A Pseudo Random Number Generator
        REM By Richard Russell, http://www.rtrussell.co.uk/, v1.00, 30-Mar-2010
        REM Based on the original C code by François Panneton & Pierre L'Écuyer
        REM (Université de Montréal) and Makoto Matsumoto (Hiroshima University)
        REM http://en.wikipedia.org/wiki/Well_Equidistributed_Long-period_Linear
 
        DIM State%(31)
        State_i% = 0
 
        REM Initialisation:
        FOR j% = 0 TO 31
          State%(j%) = RND
        NEXT j%
 
        FOR j% = 1 TO 32
          PRINT RIGHT$("0000000"+STR$~(FN_WELLRNG1024a(State_i%, State%())), 8)
        NEXT j%
        END
 
        DEF FN_WELLRNG1024a(RETURN I%, S%())
        LOCAL X%, Y%
        X% = S%(I%) EOR S%(I%+3AND31) EOR S%(I%+3AND31)>>>8
        Y% = S%(I%+24AND31) EOR S%(I%+24AND31)<<19 EOR S%(I%+10AND31) EOR S%(I%+10AND31)<<14
        S%(I%) = X% EOR Y%
        I% = I%+31AND31
        S%(I%) = S%(I%) EOR S%(I%)<<11 EOR X% EOR X%<<7 EOR Y% EOR Y%<<13
        = S%(I%)