ari@sylvester (Ari Gross) (08/31/89)
Here is some Fortran code for generating random numbers.
subroutine randu(ix, iy, yfl)
iy = ix*65539
if (iy) 5, 6, 6
5 iy = iy + 2147483647 + 1
6 yfl = iy
yfl = yfl*(.4656613e-9)
return
The number 65539 is prime and is used to continue generating
random numbers, since when this number overflows the lower
order are taken (which should be fairly random).
It ran on the IBM mainframe, but when tried on the PC
it fails to deliver. The difference seems to be in the way
the machine handle an overflow condition, which naturally
happens when this code is called repeatedly (to say generate
a gauss distribution random number).
Any suggestions as to how the code can be modified, or does anyone have
some fortran code that DOES generate random numbers for the PC
(uniformly or gaussian distribution) ???
Ari Gross
Columbia University CS
ari@cs.columbia.edujohnl@esegue.uucp (John R. Levine) (08/31/89)
In article <6474@columbia.edu> ari@sylvester.UUCP (Ari Gross) writes: >Here is some Fortran code for generating random numbers. > >subroutine randu(ix, iy, yfl) >... The RANDU routine from the old IBM SSP package is notorious as the worst alleged random number generator ever written. See the October 1988 CACM for an article describing the pitfalls of pseudo-random number generation and an algorithm known to be robust. If you want normally distributed pseudo-randoms, see Knuth's Seminumerical Algorithms for some ways to do it, since the SSP GAUSS routine is no better than RANDU. Most of the SSP routines were pretty good, but they really blew it in the random number department. -- John R. Levine, Segue Software, POB 349, Cambridge MA 02238, +1 617 492 3869 {ima|lotus}!esegue!johnl, johnl@ima.isc.com, Levine@YALE.something Massachusetts has 64 licensed drivers who are over 100 years old. -The Globe
ttw@lanl.gov (Tony Warnock) (09/07/89)
Here is a portable random number generator. It gives the same
sequence of numbers on PC's, VAX's, CRAY's, etc.
FUNCTION fran ()
c
c...portable random number generator
c
c A..........temporary partial product
c B..........temporary partial product
c GL.........left hand part of multiplier
c GR.........right hand part of multiplier
c P..........2**24
c Q..........2**(-24)
c SEED.......random number internal value
c SMALL......2**(-49)
c
DOUBLE PRECISION p, q, a, b, seed, gl, gr, small
PARAMETER (p=2.**24, q=2.**(-24), gl=1136868.d0, gr=6328637.d0)
PARAMETER (small=2.**(-49))
REAL fran
SAVE /rhyzome/
COMMON /rhyzome/ seed(2)
c
a=gr*seed(2)+1.0d0
b=gl*seed(2)+gr*seed(1)+dint(a*q)
seed(1)=b-DINT(b*q)*p
seed(2)=a-DINT(a*q)*p
fran=((seed(1)*p+seed(2))*q*q)+small
RETURN
ENDbarkdoll@cattell.psych.upenn.edu (Edwin Barkdoll) (09/12/89)
I don't have a copy of the article on hand but Brent (1974) has a routine that returns a pseudo-random (which is all you're gonna get with a computer) number from a normal distribution. I'm pretty sure the code is in fortran but I won't swear to it. R. P. Brent (1974) Algorithm 488. _Comm. ACM_ sorry I don't have the volume or page numbers but the algorithm number should make it easy to find.