checked in for jerry

[unix-history] / usr / src / usr.bin / f77 / libF77 / rand_.s.vax

# Date: 4 Mar 1982 20:54:18-PST
# From: cithep!citcsv!kingsley at Berkeley
# To: cithep!ucbvax!4bsd-bugs@Berkeley
# Subject: random number generation.

# I have done a little work with rand supplied with the system and I have
# discovered that it is flawed.  The manual page claims that it has a
# period of 2^32 and returns numbers from 0 to 2^31-1.  The code makes it
# look like the author thought it was correct, but it is not.  Instead of
# masking out the most significant (and also most random) bit, you should
# do an unsigned shift to throw out the least significant (and least
# random) bit.  I have also found a multiplier that passes Knuth's
# spectral test very well.
#       I have written a new rand, along with randint(n) which returns 0
# to n-1, and flat() which returns 0.0 to <1.0.  I did it in assembler
# (Mea Maxima Culpa!) to use the extended multiply and some bit fiddling.

# Yes, I realize that the bottom bits aren't random.  In fact, the bottom
# n bits have a period of 2^n.  The rng delivered, though, throws out the
# most significant bit to produce a 31 bit number, and claims that it has
# a period of 2^32.

# The actual generator is the routine rand, the global routines just 
# do range conversion.

# Chris Kingsley

# Adapted to f77 by David Wasley, U.C.Berkeley
# April 1983
#       @(#)rand_.s.vax 1.1

       .data
       .align  2
_randx: .long    1          # current value
_nitval:.long    1          # seed
       .text
       .align      1

       .globl  _isrand_    # set the random seed
_isrand_: .word   0         # isrand(seed) int seed;
       movl    _nitval,r0  # return old seed
       movl    *4(ap),_randx
       movl    *4(ap),_nitval
       ret
       .align  1

       .globl  _irand_     # give a 31 bit random positive integer
_irand_:.word   0           # integer rand(flag) int flag
       tstl    *4(ap)      # 0 is normal
       beql    ir1
       cmpl    $1,*4(ap)    # if arg is 1, restart
       bneq    ir0
       movl    _nitval,_randx
       jbr     ir1
ir0:    movl    *4(ap),_randx # new seed
       movl    *4(ap),_nitval # new seed
ir1:    jsb     rand
       bicl2   $1,r0
       rotl    $-1,r0,r0
       ret

       .globl  _irandn_    # give a random positive integer from 0 to n-1
_irandn_:.word  0xc         # integer irandn(n) int n;
       jsb     rand
       emul    *4(ap),r0,$0,r2
       tstl    r0
       jgeq    irn1
       addl3   *4(ap),r3,r0
       jbr     irn2
irn1:   movl    r3,r0
irn2:   ret
       .align  1

# compute the next 32 bit random number
rand:   mull3   $505360173,_randx,r0
       addl2   $907633385,r0
       movl    r0,_randx
       rsb
       .align 1

       .globl  _drand_     # give a random double from 0. to <1.
_drand_: .word   0xc        # double precision drand(flag)
dr0:    tstl    *4(ap)      # 0 is normal
       beql    dr2
       cmpl    $1,*4(ap)    # if arg is 1, restart
       bneq    dr1
       movl    _nitval,_randx
       jbr     dr2
dr1:    movl    *4(ap),_randx # new seed
       movl    *4(ap),_nitval # new seed
dr2:    jsb     rand
       movl    r0,r2
       movf    $0f1.0,r0
       extzv   $25,$7,r2,r3
       insv    r3,$0,$7,r0
       extzv   $9,$16,r2,r3
       insv    r3,$16,$16,r0
       extzv   $0,$9,r2,r1
       subd2   $0d1.0,r0
       ret

       .globl  _rand_      # fake entry for single precision rand
_rand_: .word   0xc
       jbr     dr0