Random number generators
in [Fortran]
|
Prev: pointer and allocatable arrays
Next: f95 to windows dll
From: nmm1 on 19 Nov 2009 15:38 In article <%ShNm.55991$ze1.49888(a)news-server.bigpond.net.au>, robin <robin_v(a)bigpond.com> wrote: >"David Flower" <DavJFlower(a)AOL.COM> wrote in message >news:8b69b09f-fd08-44a2-a48b-03afab751b88(a)z41g2000yqz.googlegroups.com... > >>Posters may be interested in the following reference: > >>A.C.M. Trans. Math, Software, 5, #2, 132 (1979) by Linus Schrage > >It's a bit ancient. >Those by George Marsaglia are not only portable, >they also have extremely long periods. > >His RNGs include 32-bit generators and 64-bit generators. 32-bit generators should never be used in any simulation which uses more than a million numbers in total. To a first approximation, ALWAYS use 64-bit ones. Not all of Marsaglia's generators are free from 32-bit defects, even when they are used in 64-bit forms, though most of them are good or very good. Their actual code isn't always very portable, but that's easy to clean up. Regards, Nick Maclaren.
From: glen herrmannsfeldt on 19 Nov 2009 16:25 nmm1(a)cam.ac.uk wrote: (snip) > 32-bit generators should never be used in any simulation which > uses more than a million numbers in total. To a first approximation, > ALWAYS use 64-bit ones. All the ones I knew from 40 years ago were 32 bit (or 31 bit). For many problems, that is probably enough. I have noticed that the games that come with Windows have bad RNGs. Of course closed source so we really don't know how they do it. > Not all of Marsaglia's generators are free from 32-bit defects, even > when they are used in 64-bit forms, though most of them are good or > very good. Their actual code isn't always very portable, but that's > easy to clean up. -- glen
From: frank on 19 Nov 2009 19:40 On Wed, 18 Nov 2009 12:03:34 +0000, glen herrmannsfeldt wrote: > David Jones <dajxxxx(a)ceh.ac.uk> wrote: >> steve wrote: >>> On Nov 16, 6:04 pm, frank <fr...(a)example.invalid> wrote: > (snip) > >>>> Isn't there an asymmetry in the unit interval though as to which >>>> endpoint is included? So if there's N outcomes on one side of .5 >>>> there would be N +-1 on the other. > >>> Can you restate your question without a little more detail? >> <snip> > >> The point centres on the explicit statement of the range in the >> standard as 0<= x < 1. Thus an exact value of zero is allowed as an >> output. while an exact value of one is not. So there is asymmetry on >> this point. For some uses, the question of whether or not exact zeroes >> and ones are ever returned is important and is sometimes addressed by >> testing for such values and re-generating another as necessary. >> Obviously such testing should be avoided if the zero or one case is >> guaranteed not to occur. Random number generators are such as to >> naturally include zero as a possible output, and to avoid this would >> involve slowing the code down. > > Consider a true random number generator (not a pseudo-random number > generator, as is usually used) generating an IEEE double with a 53 bit > significand generating uniformly 2**53 possible values, 0<=x<1. The > expected standard deviation of the mean after averaging N values is > 1/sqrt(N). To notice a deviation of 1 in 2**53, the offset due to the > asymmetry, requires N of about 2**104, or about 1e31. At 1GHz, that > takes 1e22 seconds, or, if I did it right, about 6e14 years. Good PRNG > generate enough additional bits internally that it will be about as hard > to notice. That's on one end of the spectrum. On the other is me, PC Joe in single precision. I've worked a lot with this material, for example when I was working up different sorts. Things look pretty glompy down here: program rand implicit integer(a-t) integer, parameter :: r = 10 real, dimension(1:r) :: myarray real u integer m call init_seed do i =1,r call random_number(u) myarray(i) = u end do m = 2**30 * 1.9 print *, "myarray is ", myarray print *, "m is ", m contains subroutine init_seed() integer :: n, ival(8), v(3), i integer, allocatable :: seed(:) call date_and_time(values=ival) v(1) = ival(8) + 2048*ival(7) v(2) = ival(6) + 64*ival(5) ! value(4) isn't really 'random' v(3) = ival(3) + 32*ival(2) + 32*8*ival(1) call random_seed(size=n) allocate(seed(n)) call random_seed() ! Give the seed an implementation-dependent kick call random_seed(get=seed) do i=1, n seed(i) = seed(i) + v(mod(i-1, 3) + 1) enddo call random_seed(put=seed) deallocate(seed) end subroutine end program rand ! gfortran ran1.f90 -Wall -Wextra -o out dan(a)dan-desktop:~/source$ gfortran ran1.f90 -Wall -Wextra -o out dan(a)dan-desktop:~/source$ ./out myarray is 0.92712229 0.85859424 0.79266590 0.57997841 0.31713003 0.38777798 0.23500705 0.40011257 0.90164834 0.72985947 m is 2040109440 dan(a)dan-desktop:~/source$ Since 2**31 overflowed, I think m is close to maybe 2.1 billion outcomes to expect on the unit interval. If m =~ 2n +-1, then n is on the order of one billion, which was only a large number while the great William Proxmire was alive. > >> A description of a RNG should explicity state the range and >> discretisation of the values that might be reurned, so that it would be >> possible to see whether just rejecting zero values would leave the set >> of allowed values symmetric about zero. > > If processor speeds double every year for the next 100 years then you > might have a chance to notice. If by then the generators switch to a > longer return value then you won't. Glen, I would almost be scared if our processor speeds trended along the same curve for the next forty years. -- frank "Guns: yes, they are harmful."
From: frank on 19 Nov 2009 19:45 On Mon, 16 Nov 2009 18:48:57 -0800, steve wrote: >> Isn't there an asymmetry in the unit interval though as to which >> endpoint is included? So if there's N outcomes on one side of .5 there >> would be N +-1 on the other. > > Can you restate your question without a little more detail? For > example, > for N=4, random_number in gfortran returns 0.998, 0.567, 0.966, 0.748. > If I draw 10000 numbers and count the number in [0,0.5) and [0.5,1), the > outcome is 4926 and 5074. If I draw 10000000 numbers, I get 4999679 and > 50000321. If I reseed the rng with a different seed, these numbers > change. At least conceptually I think it's important to get the [ ) [ ) correct on this on the test, as you have it. So the test would be if (x.gt. (.5)). -- frank "Guns: yes, they are harmful."
From: steve on 19 Nov 2009 20:22
On Nov 19, 4:45 pm, frank <fr...(a)example.invalid> wrote: > On Mon, 16 Nov 2009 18:48:57 -0800, steve wrote: > >> Isn't there an asymmetry in the unit interval though as to which > >> endpoint is included? So if there's N outcomes on one side of .5 there > >> would be N +-1 on the other. > > > Can you restate your question without a little more detail? For > > example, > > for N=4, random_number in gfortran returns 0.998, 0.567, 0.966, 0.748.. > > If I draw 10000 numbers and count the number in [0,0.5) and [0.5,1), the > > outcome is 4926 and 5074. If I draw 10000000 numbers, I get 4999679 and > > 50000321. If I reseed the rng with a different seed, these numbers > > change. > > At least conceptually I think it's important to get the [ ) [ ) > correct on this on the test, as you have it. So the test would be if > (x.gt. (.5)). With your suggested test, the intervals are [0,0.5] and (0.5,1). See Ron Shepard's post with respect to finite-precision floating point numbers. -- steve |