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.