1# Date: 4 Mar 1982 20:54:18-PST 2# From: cithep!citcsv!kingsley at Berkeley 3# To: cithep!ucbvax!4bsd-bugs@Berkeley 4# Subject: random number generation. 5 6# I have done a little work with rand supplied with the system and I have 7# discovered that it is flawed. The manual page claims that it has a 8# period of 2^32 and returns numbers from 0 to 2^31-1. The code makes it 9# look like the author thought it was correct, but it is not. Instead of 10# masking out the most significant (and also most random) bit, you should 11# do an unsigned shift to throw out the least significant (and least 12# random) bit. I have also found a multiplier that passes Knuth's 13# spectral test very well. 14# I have written a new rand, along with randint(n) which returns 0 15# to n-1, and flat() which returns 0.0 to <1.0. I did it in assembler 16# (Mea Maxima Culpa!) to use the extended multiply and some bit fiddling. 17 18# Yes, I realize that the bottom bits aren't random. In fact, the bottom 19# n bits have a period of 2^n. The rng delivered, though, throws out the 20# most significant bit to produce a 31 bit number, and claims that it has 21# a period of 2^32. 22 23# The actual generator is the routine rand, the global routines just 24# do range conversion. 25 26# Chris Kingsley 27 28# Adapted to f77 by David Wasley, U.C.Berkeley 29# April 1983 30# @(#)rand_.s.vax 1.1 31 32 .data 33 .align 2 34_randx: .long 1 # current value 35_nitval:.long 1 # seed 36 .text 37 .align 1 38 39 .globl _isrand_ # set the random seed 40_isrand_: .word 0 # isrand(seed) int seed; 41 movl _nitval,r0 # return old seed 42 movl *4(ap),_randx 43 movl *4(ap),_nitval 44 ret 45 .align 1 46 47 .globl _irand_ # give a 31 bit random positive integer 48_irand_:.word 0 # integer rand(flag) int flag 49 tstl *4(ap) # 0 is normal 50 beql ir1 51 cmpl $1,*4(ap) # if arg is 1, restart 52 bneq ir0 53 movl _nitval,_randx 54 jbr ir1 55ir0: movl *4(ap),_randx # new seed 56 movl *4(ap),_nitval # new seed 57ir1: jsb rand 58 bicl2 $1,r0 59 rotl $-1,r0,r0 60 ret 61 62 .globl _irandn_ # give a random positive integer from 0 to n-1 63_irandn_:.word 0xc # integer irandn(n) int n; 64 jsb rand 65 emul *4(ap),r0,$0,r2 66 tstl r0 67 jgeq irn1 68 addl3 *4(ap),r3,r0 69 jbr irn2 70irn1: movl r3,r0 71irn2: ret 72 .align 1 73 74# compute the next 32 bit random number 75rand: mull3 $505360173,_randx,r0 76 addl2 $907633385,r0 77 movl r0,_randx 78 rsb 79 .align 1 80 81 .globl _drand_ # give a random double from 0. to <1. 82_drand_: .word 0xc # double precision drand(flag) 83dr0: tstl *4(ap) # 0 is normal 84 beql dr2 85 cmpl $1,*4(ap) # if arg is 1, restart 86 bneq dr1 87 movl _nitval,_randx 88 jbr dr2 89dr1: movl *4(ap),_randx # new seed 90 movl *4(ap),_nitval # new seed 91dr2: jsb rand 92 movl r0,r2 93 movf $0f1.0,r0 94 extzv $25,$7,r2,r3 95 insv r3,$0,$7,r0 96 extzv $9,$16,r2,r3 97 insv r3,$16,$16,r0 98 extzv $0,$9,r2,r1 99 subd2 $0d1.0,r0 100 ret 101 102 .globl _rand_ # fake entry for single precision rand 103_rand_: .word 0xc 104 jbr dr0 105