Linear Congruential Pseudo-Random Number Generator Routines

by Bruce Clark, 7 Jun 2004

Random numbers are a complex and highly debated subject. The routines presented here use a method known as the linear congruential method. While no method of generating random numbers is perfect, the linear congruential method is widely considered to be a reasonable method. Only the basics of this method will be covered. For further information, the book "The Art Of Computer Programming, Volume 2" by Donald Knuth is highly recommended. Half of the book (chapter 3) covers the topic of random numbers. The strengths and weakness of the common random number generators are discussed, and numerous empirical and theoretical tests for randomness are described. In addition, it suggests several guidelines for random number generators. The book is quite mathematically intense, but it covers random numbers thoroughly.

A pseudo-random number generator does not return an unpredictable number, like rolling a die would. The idea is to take a long sequence of numbers and return them, one by one, in an order that seems random. The key is to start at a random point in this sequence. This starting point is called the seed. The idea is to make sequence long enough so that only a fraction of the sequence will be used when generating random numbers. Then, by starting a random place is the sequence, the fraction of the sequence used will be unlikely to overlap from one run to the next.

Note that by starting with the same seed, the same sequence of numbers can be generated repeatedly. This can be a useful feature.

On the 6502, 8-bit or 16-bit random numbers will be the most common, so the generator will return a 32-bit number. Each of the 2^32 possible values will be returned by the generator, in a seemingly random order, of course.

The term "linear congruential" sounds scary, but there are only two steps involved: multplication and addition. To get the new seed (i.e. the next number in the sequence), the old seed is multiplied by a constant, a, called the multiplier, then a second constant, c, is added to the product. The new seed is the lower 32 bits (remember, this generator returns 32-bit numbers) of the result. (Actually, keeping only the lower 32 bits is really a third step, and in fact, it is an extremely important step.) The computation performed is:

seed = a * seed + c

It turns out that the real key to making a linear congruential generator work is selecting a good value for a, the multiplier. The constant, c, is not as critical as the multiplier, but it must be an odd number. In fact, c can be chosen to be one, and this is what the routines here use, since it will simplify the calculation.

Note that an odd number (the old seed) times an odd number (the multiplier, a) is an odd number, and an odd number (the product) plus an odd number (the constant, c) will be an even number. Likewise, an even number times an odd number, plus an odd number will be an odd number. So the seed will alternate between odd and even numbers. This means that the upper bits of the seed must be used for generating 8-bit and 16-bit random numbers.

The first subroutine is called RAND and it performs the seed = a * seed + 1 (remember, c = 1 was chosen) calculation. There are three versions, a fast version, a short (in terms of number of bytes) version, and a middle version, which is between the other two versions in terms of space and speed, yet is reasonably short and reasonably fast.

The multiplier 1664525 (in hex, $19660D) is used in all three versions. This multiplier was chosen from line 16 of the table on p. 106 of "The Art Of Computer Programming, Volume 2". (That table lists various multipliers and their result in the spectral test, a critical theoretical test for measuring randomness.)

The seed is stored in SEED3 (high byte) through SEED0 (low byte). Putting SEED3 through SEED0 on the zero page will be faster and shorter (i.e. fewer bytes), but it is not necessary.

The fast version works by using four pages of pre-computed tables. There are four 256-byte tables, which are named T3, T2, T1, T0.

T3,X = byte X of T3 = bits 31-24 of 1664525 * X, (X = 0 to 255)
T2,X = byte X of T2 = bits 23-16 of 1664525 * X, (X = 0 to 255)
T1,X = byte X of T1 = bits 15-8  of 1664525 * X, (X = 0 to 255)
T0,X = byte X of T0 = bits 7-0   of 1664525 * X, (X = 0 to 255)

The product of 1664525 * SEED is:

Since only the lowest 4 bytes are kept, the multiplication is: (for brevity, SEED3 to SEED0 are called S3 to S0)

; Linear congruential pseudo-random number generator
;
; Calculate SEED = 1664525 * SEED + 1
; 
; Enter with:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; Returns:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; TMP is overwritten
;
; For maximum speed, locate each table on a page boundary
;
; Assuming that (a) SEED0 to SEED3 and TMP are located on page zero, and (b)
; all four tables start on a page boundary:
;
;   Space: 58 bytes for the routine
;          1024 bytes for the tables
;   Speed: JSR RAND takes 94 cycles
;
RAND     CLC       ; compute lower 32 bits of:
         LDX SEED0 ; 1664525 * ($100 * SEED1 + SEED0) + 1
         LDY SEED1
         LDA T0,X
         ADC #1
         STA SEED0
         LDA T1,X
         ADC T0,Y
         STA SEED1
         LDA T2,X
         ADC T1,Y
         STA TMP
         LDA T3,X
         ADC T2,Y
         TAY       ; keep byte 3 in Y for now (for speed)
         CLC       ; add lower 32 bits of:
         LDX SEED2 ; 1664525 * ($10000 * SEED2)
         LDA TMP
         ADC T0,X
         STA SEED2
         TYA
         ADC T1,X
         CLC
         LDX SEED3 ; add lower 32 bits of:
         ADC T0,X  ; 1664525 * ($1000000 * SEED3)
         STA SEED3
         RTS
;
; Generate T0, T1, T2 and T3 tables
;
; A different multiplier can be used by simply replacing the four bytes
; that are commented below
;
; To speed this routine up (which will make the routine one byte longer):
; 1. Delete the first INX instruction
; 2. Replace LDA Tn-1,X with LDA Tn,X (n = 0 to 3)
; 3. Replace STA Tn,X with STA Tn+1,X (n = 0 to 3)
; 4. Insert CPX #$FF between the INX and BNE GT1
;
GENTBLS  LDX #0      ; 1664525 * 0 = 0
         STX T0
         STX T1
         STX T2
         STX T3
         INX
         CLC
GT1      LDA T0-1,X  ; add 1664525 to previous entry to get next entry
         ADC #$0D    ; byte 0 of multiplier
         STA T0,X
         LDA T1-1,X
         ADC #$66    ; byte 1 of multiplier
         STA T1,X
         LDA T2-1,X
         ADC #$19    ; byte 2 of multiplier
         STA T2,X
         LDA T3-1,X
         ADC #$00    ; byte 3 of multiplier
         STA T3,X
         INX         ; note: carry will be clear here
         BNE GT1
         RTS

In the RAND routine, the ADC #1 instruction can be eliminated by calculating 1664525 * X + 1 and storing the result in tables T3I through T0I. SEED0 will then use TnI (n = 0 to 3) and SEED1 through SEED3 will use Tn. Since table T3 was only used by SEED0, it (T3) will no longer be necessary. Table T2 will be identical to T2I, so table T2 can also be eliminated. This means only 6 pages of tables (T0, T1, T0I, T1I, T2I, and T3I) are needed and that JSR RAND will take 92 cycles.

The short version is just an ordinary 32-bit * 32-bit multiplication routine. The DB pseudo-op is called .BYTE on some assemblers. Consult the assembler documentation for the pseudo-op name it expects.

; Linear congruential pseudo-random number generator
;
; Calculate SEED = 1664525 * SEED + 1
; 
; Enter with:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; Returns:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; TMP, TMP+1, TMP+2 and TMP+3 are overwritten
;
; Note that TMP to TMP+3 and RAND6 are high byte first, low byte last
;
; Assuming that (a) SEED0 to SEED3 and TMP+0 to TMP+3 are all located on page
; zero, and (b) none of the branches cross a page boundary:
;
;   Space: 53 bytes
;   Speed: JSR RAND takes 2744 cycles, on average (1624 to 3864 cycles)
;          specifically, JSR RAND takes 1624 + 70 * N cycles, where
;          N = number of bits of SEED that are 1
;
RAND     LDA #1              ; store 1 in TMP
         LDX #3
RAND1    STA TMP,X
         LSR
         DEX
         BPL RAND1
         LDY #$20            ; calculate SEED = SEED * RAND4 + TMP
         BNE RAND5           ; branch always
RAND2    BCC RAND4           ; branch if a zero was shifted out
         CLC                 ; add multiplier to product
         LDX #3
RAND3    LDA TMP,X
         ADC RAND4,X
         STA TMP,X
         DEX
         BPL RAND3
RAND4    ROR TMP             ; shift result right
         ROR TMP+1
         ROR TMP+2
         ROR TMP+3
RAND5    ROR SEED3           ; shift out old seed, and shift in new seed
         ROR SEED2
         ROR SEED1
         ROR SEED0
         DEY
         BPL RAND2
         RTS
RAND6    DB  $00,$19,$66,$0D ; multiplier (high byte first!)

Since the carry is known to be set after the BCC RAND4, the CLC can be eliminated by adding $0019660C instead of $0019660D, saving one byte.

The routine will take one fewer byte if RAND6 (the multiplier) is located on the zero page. However, since the zero page is usually RAM, RAND6 would have to be initialized somehow, which would likely nullify the one byte savings.

The middle version is similar to the typical 6502 non-looping multiply-by-10 routines, which are tailored to multiply by a specific number. Unfortunately, this routine must be rewritten if a different multiplier is to be used. The idea to add seed, $100 * seed, $10000 * seed, and/or $1000000 * seed when necessary (as determined by the multiplier), then shift seed (left), and repeat. Specifically,

Bit 0 of $0D (multiplier byte 0) is 1
Bit 0 of $66 (multiplier byte 1) is 0
Bit 0 of $19 (multiplier byte 2) is 1
Bit 0 of $00 (multiplier byte 3) is 0

so, seed is added to $10000 * seed, and seed is shifted left (multiplied by 2)

Bit 1 of $0D is 0
Bit 1 of $66 is 1
Bit 1 of $19 is 0
Bit 1 of $00 is 0

so, $100 * seed is added to the previous result, and seed is shifted left

Bit 2 of $0D is 1
Bit 2 of $66 is 1
Bit 2 of $19 is 0
Bit 2 of $00 is 0

so $100 * seed + seed is added to the previous result, and so on.

; Linear congruential pseudo-random number generator
;
; Calculate SEED = 1664525 * SEED + 1
; 
; Enter with:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; Returns:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; TMP, TMP+1, TMP+2 and TMP+3 are overwritten
;
; Assuming that (a) SEED0 to SEED3 and TMP+0 to TMP+3 are all located on page
; zero, and (b) none of the branches cross a page boundary:
;
;   Space: 106 bytes
;   Speed: JSR RAND takes 517 cycles
;
RAND     CLC         ; copy SEED into TMP
         LDA SEED0   ; and compute SEED = SEED * $10000 + SEED + 1
         STA TMP
         ADC #1
         STA SEED0
         LDA SEED1
         STA TMP+1
         ADC #0
         STA SEED1
         LDA SEED2
         STA TMP+2
         ADC TMP
         STA SEED2
         LDA SEED3
         STA TMP+3
         ADC TMP+1
         STA SEED3
;
; Bit 7 of $00, $19, $66, and $0D is 0, so only 6 shifts are necessary
;
         LDY #5
RAND1    ASL TMP     ; shift TMP (old seed) left
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
;
; Get X from the RAND4 table.  When:
;
; X = $00, SEED = SEED + $10000 * TMP
; X = $01, SEED = SEED + $100 * TMP
; X = $FE, SEED = SEED + $10000 * TMP + TMP
; X = $FF, SEED = SEED + $100 * TMP + TMP
;
         LDX RAND4,Y
         BPL RAND2   ; branch if X = $00 or X = $01
         CLC         ; SEED = SEED + TMP
         LDA SEED0
         ADC TMP
         STA SEED0
         LDA SEED1
         ADC TMP+1
         STA SEED1
         LDA SEED2
         ADC TMP+2
         STA SEED2
         LDA SEED3
         ADC TMP+3
         STA SEED3
         INX         ; $FE -> $00, $FF -> $01
         INX
RAND2    CLC
         BEQ RAND3   ; if X = $00, SEED = SEED + TMP * $10000
         LDA SEED1   ; SEED = SEED + TMP * $100
         ADC TMP
         STA SEED1
RAND3    LDA SEED2
         ADC TMP,X
         STA SEED2
         LDA SEED3
         ADC TMP+1,X
         STA SEED3
         DEY
         BPL RAND1
         RTS
RAND4    DB  $01,$01,$00,$FE,$FF,$01

6 cycles (and 1 byte) can be saved by locating the RAND4 table on the zero page. Since the zero page is usually RAM, RAND4 would have to be initialized.

Here is an example of the middle routine rewritten for the multiplier 69069 ($10DCD). This multiplier was taken from line 15 of the table on p. 106 of "The Art Of Computer Programming, Volume 2". (According to that table, this multiplier does not do as well on the spectral test as the multiplier 1664525.)

; Linear congruential pseudo-random number generator
;
; Calculate SEED = SEED * 69069 + 1
; 
; Enter with:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; Returns:
;
;   SEED0 = byte 0 of seed
;   SEED1 = byte 1 of seed
;   SEED2 = byte 2 of seed
;   SEED3 = byte 3 of seed
;
; TMP, TMP+1, TMP+2 and TMP+3 are overwritten
;
; Assuming that SEED0 to SEED3 and TMP+0 to TMP+3 are all located on page
; zero:
;
;   Space: 173 bytes
;   Speed: JSR RAND takes 326 cycles
;
RAND     LDA SEED0 ; TMP = SEED * 2
         ASL
         STA TMP
         LDA SEED1
         ROL
         STA TMP+1
         LDA SEED2
         ROL
         STA TMP+2
         LDA SEED3
         ROL
         STA TMP+3
         CLC       ; TMP = TMP + SEED (= SEED * 3)
         LDA SEED0
         ADC TMP
         STA TMP
         LDA SEED1
         ADC TMP+1
         STA TMP+1
         LDA SEED2
         ADC TMP+2
         STA TMP+2
         LDA SEED3
         ADC TMP+3
         STA TMP+3
         CLC       ; SEED = SEED + $10000 * SEED
         LDA SEED2
         ADC SEED0
         TAX       ; keep byte 2 in X for now (for speed)
         LDA SEED3
         ADC SEED1
         TAY       ; keep byte 3 in Y for now
         CLC       ; SEED = SEED + $100 * SEED
         LDA SEED1
         ADC SEED0
         PHA       ; push byte 1 onto stack
         TXA
         ADC SEED1
         TAX
         TYA
         ADC SEED2
         TAY
         LDA TMP   ; TMP = TMP * 4 (= old seed * $0C)
         ASL
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
         ASL
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
         STA TMP
         CLC       ; SEED = SEED + TMP
         ADC SEED0
         STA SEED0
         PLA       ; pull byte 1 from stack
         ADC TMP+1
         STA SEED1
         TXA
         ADC TMP+2
         TAX
         TYA
         ADC TMP+3
         TAY
         CLC
         LDA TMP   ; SEED = SEED + TMP * $100
         ADC SEED1
         STA SEED1
         TXA
         ADC TMP+1
         TAX
         TYA
         ADC TMP+2
         TAY
         LDA TMP   ; TMP = TMP * $10 (= old seed * $C0)
         ASL       ; put byte 0 of TMP in the accumulator
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
         ASL
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
         ASL
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
         ASL
         ROL TMP+1
         ROL TMP+2
         ROL TMP+3
         SEC       ; SEED = SEED + TMP + 1
         ADC SEED0
         STA SEED0
         LDA TMP+1
         ADC SEED1
         STA SEED1
         TXA
         ADC TMP+2
         STA SEED2
         TYA
         ADC TMP+3
         STA SEED3
         RTS

Having selected a RAND routine, the next step is to turn the new seed into a (pseudo-)random number. For a random number between 0 and 255 ($FF) inclusive, simply use SEED3. For a random number between 0 and 65535 ($FFFF), use SEED2 for the low byte and SEED3 for the high byte. For a random number between 0 and 1, use bit 7 of SEED3. For a random number between 0 and 3, use bits 7 through 6 of SEED3. For a random number between 0 and 7, use bits 7 through 5 of SEED3, and so on. Any random number between 0 and one less than a power of 2 is easily obtained, but what about a random number between 0 and 5, say? The solution is to multiply SEED by 6 (5 plus one) and keep the upper 8 bits of the 40-bit product. In this case, 6 is known as the modulus.

The result is the RANDOM subroutine, which returns a random number, RND, where 0 <= RND < MOD, and MOD is the modulus. There are two versions of RANDOM, one for 8-bit RND and MOD, and one for 16-bit RND and MOD. The 8-bit version, called RANDOM8, is up first, and it uses the number in the accumulator as the modulus, and returns the random number in the accumulator.

; Linear congruential pseudo-random number generator
;
; Get the next SEED and obtain an 8-bit random number from it
;
; Requires the RAND subroutine
;
; Enter with:
;
;   accumulator = modulus
;
; Exit with:
;
;   accumulator = random number, 0 <= accumulator < modulus
;
; MOD, TMP, TMP+1, and TMP+2 are overwritten
;
; Note that TMP to TMP+2 are only used after RAND is called.
;
RANDOM8  STA MOD    ; store modulus in MOD
         JSR RAND   ; get next seed
         LDA #0     ; multiply SEED by MOD
         STA TMP+2
         STA TMP+1
         STA TMP
         SEC
         ROR MOD    ; shift out modulus, shifting in a 1 (will loop 8 times)
R8A      BCC R8B    ; branch if a zero was shifted out
         CLC        ; add SEED, keep upper 8 bits of product in accumulator
         TAX
         LDA TMP
         ADC SEED0
         STA TMP
         LDA TMP+1
         ADC SEED1
         STA TMP+1
         LDA TMP+2
         ADC SEED2
         STA TMP+2
         TXA
         ADC SEED3
R8B      ROR        ; shift product right
         ROR TMP+2
         ROR TMP+1
         ROR TMP
         LSR MOD    ; loop until all 8 bits of MOD have been shifted out
         BNE R8A
         RTS

Here is the 16-bit version of RANDOM, called RANDOM16.

; Linear congruential pseudo-random number generator
;
; Get the next SEED and obtain an 16-bit random number from it
;
; Requires the RAND subroutine
;
; Enter with:
;
;   MOD = modulus
;
; Exit with:
;
;   RND = random number, 0 <= RND < MOD
;
; TMP is overwritten, but only after RND is called.
;
RANDOM16 JSR RAND  ; get next seed
         LDA #0    ; multiply SEED by MOD
         STA RND+1
         STA RND
         STA TMP
         LDY #16
R16A     LSR MOD+1 ; shift out modulus
         ROR MOD
         BCC R16B  ; branch if a zero was shifted out
         CLC       ; add SEED, keep upper 16 bits of product in RND
         ADC SEED0
         TAX
         LDA TMP
         ADC SEED1
         STA TMP
         LDA RND
         ADC SEED2
         STA RND
         LDA RND+1
         ADC SEED3
         STA RND+1
         TXA
R16B     ROR RND+1 ; shift product right
         ROR RND
         ROR TMP
         ROR
         DEY
         BNE R16A
         RTS

There are 2^32 possible values for SEED. When using a modulus that is a power of 2, each possible value (of RND) returned by RANDOM is equally likely. However, when modulus of 6 is used, for example, some of the possible values are a teensy (a very, very teensy) bit more likely than others, since 2^32 is not a multiple of 6. But 2^32 - 4 is a multiple of 6, so one solution is to reject four possible values of seed by simply calling RAND again to get the next seed when one of those four possibilities is encountered. Since only four values out of 2^32 get rejected, it is highly unlikely that two consecutive seed values will be rejected, so at most, RAND will be called twice. The right sequence can ensure that this is so. The question is, which four values should be rejected?

To answer this question, a simpler example is in order, specifically a 4-bit seed, and a modulus of 7. Since 16-2 is multiple of 7, two values must be rejected. Since there are only 16 possible seed values, all 16 possibilites will be listed, along with the corresponding seed * modulus product.

seed seed * 7
---- --------
0000 000 0000
0001 000 0111
0010 000 1110
0011 001 0101
0100 001 1100
0101 010 0011
0110 010 1010
0111 011 0001 
1000 011 1000 
1001 011 1111 
1010 100 0110 
1011 100 1101
1100 101 0100
1101 101 1011
1110 110 0010
1111 110 1001

The random numbers (the upper 3 bits of seed * modulus) 0 and 3 occur three times, and the rest occur twice. Notice the lower 4 bits of seed * modulus. When the lower 4 bits are less than 2 (the number of seed values to reject), the upper 3 bits are 000 or 011. Also, when the lower 4 bits are greater than or equal to 16-2, the upper three bits are 000 or 011. This means than there are two ways to determine whether to get the next seed (a) if the lower bits of seed * modulus are less than the number of seed values to reject, or (b) if the lower bits are greater than or equal to 16 minus the number of seed values to reject.

The latter method will be used because it can be implemented little faster. For the 32-bit case, this means checking when the lower 32 bits of the seed * modulus product is less than 2^32 minus the number of seed values to reject. This will be checked by adding the number of seed values to reject to the lower 32 bits of the product and looking for a carry from the 32-bit addition.

The only remaining question is, how many seed values should be rejected? This can be calculated by dividing 2^32 by the modulus. The remainder from this division is the number of seed values to reject.

The subroutine that does all of this is called URANDOM, because it returns a random number with a uniform distribution. Like RANDOM, it returns a random number, RND, where 0 <= RND < MOD. There are also two versions, one for 8-bit RND and MOD, and one for 16-bit RND and MOD. The 8-bit version, called URANDOM8, is once again up first, and it uses the number in the accumulator as the modulus, and returns the random number in the accumulator.

; Linear congruential pseudo-random number generator
;
; Get the next SEED and obtain an 8-bit uniform random number from it
;
; Requires the RAND subroutine
;
; Enter with:
;
;   accumulator = modulus
;
; Exit with:
;
;   accumulator = random number, 0 <= accumulator < modulus
;
; MOD, REM, TMP, TMP+1, TMP+2, and TMP+3 are overwritten
;
; Note that TMP to TMP+3 are only used after RAND is called.
;
URANDOM8 STA MOD   ; store modulus in MOD
         LDX #32   ; calculate remainder of 2^32 / MOD
         LDA #1
         BNE UR8B
UR8A     ASL       ; shift dividend left
         BCS UR8C  ; branch if a one was shifted out
UR8B     CMP MOD
         BCC UR8D  ; branch if partial dividend < MOD
UR8C     SBC MOD   ; subtract MOD from partial dividend
UR8D     DEX
         BPL UR8A
         STA REM   ; store remainder in REM
UR8E     JSR RAND
         LDA #0    ; multiply SEED by MOD
         STA TMP+3
         STA TMP+2
         STA TMP+1
         STA TMP
         LDY MOD   ; save MOD in Y
         SEC
         ROR MOD   ; shift out modulus, shifting in a 1 (will loop 8 times)
UR8F     BCC UR8G  ; branch if a zero was shifted out
         CLC       ; add SEED to TMP
         TAX
         LDA TMP
         ADC SEED0
         STA TMP
         LDA TMP+1
         ADC SEED1
         STA TMP+1
         LDA TMP+2
         ADC SEED2
         STA TMP+2
         LDA TMP+3
         ADC SEED3
         STA TMP+3
         TXA
UR8G     ROR TMP+3 ; shift product right
         ROR TMP+2
         ROR TMP+1
         ROR TMP
         ROR
         LSR MOD   ; loop until all 8 bits of MOD have been shifted out
         BNE UR8F
         CLC       ; add remainder to product
         ADC REM
         BCC UR8H  ; branch if no 8-bit carry
         INC TMP   ; carry a one to byte 1 of product
         BNE UR8H  ; branch if no 16-bit carry
         INC TMP+1 ; carry a one to byte 2 of product
         BNE UR8H  ; branch if no 24-bit carry
         INC TMP+2 ; carry a one to byte 3 of product
         STY MOD   ; restore MOD (does not affect Z flag!)
         BEQ UR8E  ; branch if 32-bit carry
UR8H     LDA TMP+3 ; return upper 8 bits of product in accumulator
         RTS

Here is the 16-bit version of URANDOM, called URANDOM16.

; Linear congruential pseudo-random number generator
;
; Get the next SEED and obtain an 16-bit random number from it
;
; Requires the RAND subroutine
;
; Enter with:
;
;   MOD = modulus
;
; Exit with:
;
;   RND = random number, 0 <= RND < MOD
;
; REMH, REML, TEMP, TMP, TMP+1, TMP+2, and TMP+3 are overwritten
;
; Note that TMP to TMP+3 are only used after RAND is called.
;
URANDOM16
         LDX #32   ; calculate 2^32 / MOD
         LDA #0
         STA REMH
         SEC       ; prepare to shift in a 1
UR16A    ROL       ; shift dividend left
         ROL REMH
         BCS UR16B ; branch if a one was shifted out
         TAY       ; compare partial dividend to MOD
         CMP MOD
         LDA REMH
         SBC MOD+1
         TYA
         BCC UR16C ; branch if partial dividend < MOD
UR16B    SBC MOD   ; subtract MOD from paritial dividend (computes remainder)
         TAY
         LDA REMH
         SBC MOD+1
         STA REMH
         TYA
         CLC
UR16C    DEX
         BPL UR16A
         STA REML  ; store low byte of remainder in REM
         LDA MOD+1 ; save high byte of MOD in TEMP
         STA TEMP
UR16D    JSR RAND
         LDA MOD   ; save low byte of MOD in TMP+3
         STA TMP+3
         LDA #0    ; multiply SEED by MOD
         STA RND+1
         STA RND
         STA TMP+2
         STA TMP+1
         LDY #16
UR16E    LSR MOD+1 ; shift out modulus
         ROR MOD
         BCC UR16F ; branch if a zero was shifted out
         CLC       ; add SEED, keep upper 16 bits of product in RND
         TAX
         LDA TMP+1
         ADC SEED0
         STA TMP+1
         LDA TMP+2
         ADC SEED1
         STA TMP+2
         LDA RND
         ADC SEED2
         STA RND
         LDA RND+1
         ADC SEED3
         STA RND+1
         TXA
UR16F    ROR RND+1 ; shift product right
         ROR RND
         ROR TMP+2
         ROR TMP+1
         ROR TMP
         ROR
         DEY
         BNE UR16E
         CLC       ; add remainder to product
         ADC REML
         LDA TMP
         ADC REMH
         BCC UR16G ; branch if no 16-bit carry
         INC TMP+1 ; carry a one to byte 2 of product
         BNE UR16G ; branch if no 24-bit carry
         LDA TMP+3 ; restore MOD
         STA MOD
         LDA TEMP
         STA MOD+1
         INC TMP+2 ; carry a one to byte 3 of product
         BEQ UR16D ; branch if 32-bit carry
UR16G    RTS

URANDOM16 could be sped up slightly by replacing the SEC instruction with a LDA #1 instruction followed by a branch (BNE) to the first TAY instruction. Then, the ROL at UR16A can be changed to an ASL, and the CLC just before UR16C can be deleted.

Finally, here is a routine that allows the RAND routine to be accessed via the USR function in EhBASIC. Rbyte4 through RByte1 (the seed for the RND function) could be used for SEED3 through SEED0 as follows:

SEED0 = Rbyte4
SEED3 = SEED0+1
SEED2 = SEED0+2
SEED1 = SEED0+3

Note that Rbyte4 through Rbyte1 are not in order in memory! The downside to using Rbyte is if RND and USR are both used. The RND seed (Rbyte) cannot be all zeros, or RND will always return 0, but an all zero seed is perfectly fine for a linear congruential generator. (EhBASIC cleverly avoids an all zeros seed by setting the seed only on a non-zero argument to RND.) SEED3 through SEED0 could use unused zero page locations (such as $DE to $E1) or other unused RAM elsewhere instead.

FAC1_1 is used for TMP, since it and the next 3 locations (FAC1_2, FAC1_3, and FAC1_s) are overwritten after the JSR RAND anyway. FAC1_e could also be used since the next 3 locations are FAC1_1, FAC1_2, and FAC1_3.

Since the usage of USR is identical to that of RND, RND can simply be replaced by USR. Remember, before using USR, the statement DOKE 11,address (where address is the address of the routine below) must be used to set the USR vector.

; USR-based linear congruential pseudo-random number generator
;
; Usage is identical to RND
;
; Requires the RAND subroutine
;
TMP     =       FAC1_1
	LDA	FAC1_e		; get FAC1 exponent
        BEQ     USRNG_1         ; do next random # if zero
        LDX     #<SEED0         ; set PRNG pointer low byte
        LDY     #>SEED0         ; set PRNG pointer high byte
        JSR     LAB_2778        ; pack FAC1 into (XY)
USRNG_1
        JSR     RAND            ; get the next random number
        LDX     #$02            ; three bytes to copy
USRNG_2
        LDA     SEED3,X         ; get PRNG byte
        STA     FAC1_1,X        ; save FAC1 byte
	DEX
        BPL     USRNG_2         ; loop if not complete
        LDA     SEED0           ; set the rounding byte
        STA     FAC1_r
        LDA     #$80            ; set the exponent
	STA	FAC1_e		; save FAC1 exponent
        ASL                     ; clear A
	STA	FAC1_s		; save FAC1 sign
        JMP     LAB_24D5        ; normalize FAC1 & return

Last page update: January 30th, 2005.