CompuServe Messages

Caleb, the following is the algorithm of the random number generator used by Turbo Pascal for the IBM PC. I've run some statistical tests on it and found it not too bad; HOWEVER, I am not a statistician and the tests were based on my "common sense" ideas of how a random number generator should work, not on standardized tests. I've heard that it does ok, but not wonderfully, on such tests. Anyway: Assume that a 32 bit random number seed is stored in location SEED. To compute a 16-bit "random" integer [Random()]: Multiply SEED by $81 (ignore overflow) Add $361962E9 to SEED (ditto) Return the high 16 bits of the new SEED as the random number To compute a "random" number from 0 to N-1 [RandInt()]: Call Random() Return the 16-bit value modulus N There is a >definite< problem with RandInt(). Imagine that N is 50000. Then if a value from 0 to 49999 is returned from Random(), that value is passed back unchanged. If 50000-65535 is returned, the 0 to N-1 value will be in the range 0-15535. Over time, approximately twice as many random numbers of the range 0-15535 will be returned, as compared to numbers in the range 15536-49999. To solve this, I propose that if N is not a power of 2, Random() should be called repeatedly until the returned value is in the range 0-(N*M-1), where N*M is the largest multiple of N that is less than 65536. RandInt(50000) would reject values 50000-65535, RanInt(1000) would reject values 65000-65535. This should not be terribly hard to implement. Note that the multiplication and addition constants, $81 and $361962E9, are critical to proper operation of the generator — a few other pairs of constants may be equally good, but none will be better, and most will be MUCH worse. This pair of constants causes the seed to step through all possible 32-bit values before returning to the initial seed. Even distributions of the 65536 possible 16-bit return values are generated. The low word of the seed is of poor quality as a random value. The full 32 bit value >may< be of use as a mantissa for a random real number, but it would have to be thoroughly tested. Oh heck… I sat down and coded the routines. CAVEATS: first, I'm just a beginner are 68000 assembly; therefore, these routines may have bugs and/or be less efficient than possible. Second, I did NOT test them — I'm far from comfortable with the Amiga assembler yet; this code is mind-executed only. If it's wrong, well, at least it should provide a starting point for an experienced 68K programmer. Please someone let me know how I did. I tried to figure a way to avoid the first DIVU, but to no avail. I still suspect that it can be done. Discarding all values >=D1.W is not acceptable for two reasons: too slow (and too variable in speed), and it can eat too many cycles of the random number generator, spoiling its randomness. I provide no code to initialize Seed. It should probably be initialized from the system time & date, suitably munged to cram all the bits together into one longword. Caleb: if you're not comfortable with assembly linked to "C", don't panic! Anyone: would you apply any necessary fixes to my code, including whatever is necessary to make it linkable to "C", and upload a linkable object file? If applicable, one for Lattice and one for Aztec? Thanks… Seed DC.L 0 ; Seed value, Random()-private variable ; No register inputs ; Output register: D0.W = 16-bit random value Random: MOVE.L Seed,D0 ; D0=Seed LSL.L #7,D0 ; D0=Seed * $80 ADD.L Seed,D0 ; D0=Seed * $81 ADD.L #$361962E9,D0 ; D0=Seed * $81 + $361962E9 MOVE.L D0,Seed ; New Seed is saved SWAP D0 ; Return the high word RTS ; New random value returned in D0.W ; Input register: D1.W = 16-bit ceiling value ; Output register: D0.W = 16-bit random number in the range 0 to (D1.W-1) ; D2 is destroyed — add code to save it if appropriate RandInt:MOVE.L #$0000FFFF,D2 ; D2=65535 TST.W D1 ; RandInt(0)? (==RandInt(65536)) BEQ .1 ; If RandInt(0), max allowable value is 65535 ; (i.e., all values are allowable) (avoid 0-div) INC.L D2 ; D2=65536 DIVU D1,D2 ; D2=65536/D1 SWAP D2 ; D2.W=the remainder of 65536/D1 NOT.W D2 ; D2.W=the max allowable Random value .1: JSR Random ; Get a random value CMP.W D0,D2 ; Is it within the given bounds? BHI .1 ; If not, try Random() again SWAP D0 CLR.W D0 SWAP D0 ; Clear high word of D0, prepare for 32-bit div DIVU D1,D0 SWAP D0 ; D0.W=the remainder of Random()/D1 RTS ; Random value, 0 to (D1.W-1), is in D0.W – Bela