|
|
1.1 ! root 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 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 ! 55: ir0: movl *4(ap),_randx # new seed ! 56: movl *4(ap),_nitval # new seed ! 57: ir1: 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 ! 70: irn1: movl r3,r0 ! 71: irn2: ret ! 72: .align 1 ! 73: ! 74: # compute the next 32 bit random number ! 75: rand: 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) ! 83: dr0: 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 ! 89: dr1: movl *4(ap),_randx # new seed ! 90: movl *4(ap),_nitval # new seed ! 91: dr2: 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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.