Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Newbie: problems with random_seed

205 views
Skip to first unread message

Kamaraju Kusumanchi

unread,
Sep 23, 2003, 3:49:33 PM9/23/03
to
Hi all
I am trying to write a program which needs to generate different random
numbers sequence from time to time. I am having problems with
randomizing the seed. I searched the archives for this. I only
understood that the problem is faced by many but no "complete solution"
was posted as of today.

Consider the program

program random3
implicit none
real, dimension(10) :: randArray

call random_seed()
call random_number(randArray)
write(*,*) randArray
call random_seed()
call random_number(randArray)
write(*,*) randArray
end program random3

with redhat linux 9, absoft f90 v 8.0 the result is

0.139087 0.633978 0.236929 0.430455 0.100457
0.475715 0.167549 0.448885 2.444941E-02 0.699719

0.538699 0.378019 0.417190 0.931743 0.545394
0.851934 0.217492 0.947982 0.202268 0.478387

with redhat linux 9, intel ifc Version 7.0 the result is
0.7732671 0.8232512 0.1481206 0.6567149 0.8420556
0.6285801 0.3304787 0.7688151 0.6810636 0.1044316

0.7732671 0.8232512 0.1481206 0.6567149 0.8420556
0.6285801 0.3304787 0.7688151 0.6810636 0.1044316

In words, ifc is not changing the seed and the above code thus becomes
non-portable. I am looking for some code in fortran which accomplishes
the portability. I know how to do this in other languages as C etc.,
(using srand(time(NULL)) say). But my question is how to do this in
fortran 90 in a less cumbersome way? Being the language of "numerical
computing" I would find it wierd if there no standardized way of doing
this simple task. What do you say?

thanks in advance
raju

Michael Metcalf

unread,
Sep 23, 2003, 4:05:57 PM9/23/03
to

"Kamaraju Kusumanchi" <kk...@cornell.edu> wrote in message
news:bkq84r$fn6$1...@news01.cit.cornell.edu...

>>
> In words, ifc is not changing the seed and the above code thus becomes
> non-portable. I am looking for some code in fortran which accomplishes
> the portability. I know how to do this in other languages as C etc.,
> (using srand(time(NULL)) say). But my question is how to do this in
> fortran 90 in a less cumbersome way? Being the language of "numerical
> computing" I would find it wierd if there no standardized way of doing
> this simple task. What do you say?
>

The results you have are both allowed by the standard. The following
approach allows you to restart the sequence from where you left off by
storing the seed. In that way you have just one sequence that you can
interrupt as you wish (see also "Fortran 90/95 Explained", Section 8.16.3).

Hope that helps,

Mike Metcalf

program random3
implicit none
real, dimension(10) :: randArray

integer, allocatable, dimension(:) :: seed
integer :: size

call random_seed(size=size) ! get size of seed
allocate(seed(size)) ! allocate its space


call random_seed()
call random_number(randArray)
write(*,*) randArray

call random_seed(get=seed) ! get last value of seed
call random_seed(put=seed) ! restore seed to last value

Kamaraju Kusumanchi

unread,
Sep 23, 2003, 4:40:11 PM9/23/03
to

Michael Metcalf wrote:
> "Kamaraju Kusumanchi" <kk...@cornell.edu> wrote in message
> news:bkq84r$fn6$1...@news01.cit.cornell.edu...
>

> The results you have are both allowed by the standard. The following
> approach allows you to restart the sequence from where you left off by
> storing the seed. In that way you have just one sequence that you can
> interrupt as you wish (see also "Fortran 90/95 Explained", Section 8.16.3).
>
> Hope that helps,
>
> Mike Metcalf
>
> program random3
> implicit none
> real, dimension(10) :: randArray
> integer, allocatable, dimension(:) :: seed
> integer :: size
>
> call random_seed(size=size) ! get size of seed
> allocate(seed(size)) ! allocate its space
> call random_seed()
> call random_number(randArray)
> write(*,*) randArray
> call random_seed(get=seed) ! get last value of seed
> call random_seed(put=seed) ! restore seed to last value
> call random_number(randArray)
> write(*,*) randArray
> end program random3

Thanks a lot. This code is working for me. I want to know one more
thing. If I want to change the seed, what would be a good way to do it?
This would be useful, for example, as I want to generate different
random series for each run of the program. It is too boring to see the
same sequence of numbers and calling it random! Any ideas on this also?
I looked into section 8.16 of "Fortran 90/95 explained", but could not
information on how to do this.

For eg if size in the above program is 47. Is initializing each element
of seed (from 1 to 47) to a value returned by date_and_time a good idea?
or is there any other efficient method?

Any code in this regard is very useful.
regards
raju

keith bierman

unread,
Sep 23, 2003, 5:10:00 PM9/23/03
to
Kamaraju Kusumanchi wrote:
Is initializing each element
> of seed (from 1 to 47) to a value returned by date_and_time a good idea?
> or is there any other efficient method?

Often people use a time of day clock to set their random seed.
However, this may not result in very random numbers (viz. if your
job turns out to run at the same time every day ...).

If you want to ensure quality of results, you may want to use the
technique that Michael suggested, always storing the "last" value
and restarting from there. Thus you will have the assurance that the
numbers meet whatever criteria the authors of the RNG proved ;>

Or if your system has some sort of "entropy" you can access; some
systems provide a /dev/random, or if your system permits read access
to the entire memory state of the system you can combine that with
time of day to get interesting results (unless it's a single user
system, and you reboot regularly ;>).

Without knowing how much and what sort of *quality* of pseudo random
number one is seeking, it's really hard to come up with a maximally
efficient scheme (after all, 42 is always a good starting spot if
you don't care about the quality of the result ;>).

Are your runs really so short that the cost of computing a suitable
seed impacts your runtime?


--
Keith H. Bierman keith....@Sun.COM| 650-352-4432 voice+fax
Sun Microsystems Laboratories | sun internal 68207
15 Network Circle UMPK 15-224 |
Menlo Park, California 94025 | kbie...@acm.org
<speaking for myself, not Sun*> Copyright 2003

Richard Maine

unread,
Sep 23, 2003, 5:27:49 PM9/23/03
to
Kamaraju Kusumanchi <kk...@cornell.edu> writes:
[on the subject of pseudo-random number generators]

> I am looking for some code in fortran which accomplishes
> the portability.

Others have addressed several aspects of yuor questions at least as
well as I would have. One salient point I haven't noticed made however.

If you want complete portability, you will need to use something
other than the intrinsic procedures provided by the compiler. There
are several excellent and free random generators regularly cited
here.

With the intrinsic procedures, you can make your source code portable.
However, what you cannot do with the intrinsic procedures is generate
the same sequence of pseudo-random numbers with different compilers.
Different compilers use different algorithms for their intrinsic
random number generators. Sometimes such cross-platform repeatability
is important for such things as verification and validation. If
you need portability in the sense of generating the same sequence
of numbers with different compilers, then you will need to use one of
the third-party generators.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Kamaraju Kusumanchi

unread,
Sep 23, 2003, 6:09:19 PM9/23/03
to
No I dont need to generate the same sequence of random numbers. The
results cited in my first message seem to be very non-intuitive, they
highlight another important aspect.

One compiler changes the seed every time random_seed() is called. The
other compiler hooks on to the same seed irrespective of howmany times
the random_seed() subroutine is called. I felt this to be very
confusing. Why is it so "loosely" defined? On a personal note, I cant
see any advantage by leaving "what random_seed() should do" in an
ambiguous state.

By portability I meant the same functionality of random_seed() on
different compilers.
--raju

Richard Maine

unread,
Sep 23, 2003, 6:32:19 PM9/23/03
to
Kamaraju Kusumanchi <kk...@cornell.edu> writes:

> Why is it [random seed] so "loosely" defined?

Such "why" questions are notoriously difficult to answer. The
"usual" and uninformative answer is "because that's the way the
vote came out". That answer is all to apt in this case.

My understanding is that the "real" reason was simply imprecise
technical writing. The author of the material had intended to
specify things more precisely, but didn't manage to express his
intent in the words that actually got written down and eventually
passed as the standard.

That author has subsequently commented, in my personal presense,
that the standard doesn't say what he meant for it to say. He
has also made proposals to change it retroactively.

However, the expressed intention of a single person, even the person
who authored the material, doesn't override the formal reviews and
votes that passed the standard. There is no guarantee, just because
that person says he had a different intent, that a majority would have
agreed with that intent. A majority passed the words that are in the
standard.

By now, there are (lots of) existing implementations, existing codes
that depend on those implementations, and resulting resistance to
changing it. At least one proposal to change it retroactively failed.

I take this as a sample of a lesson that many technical people
are slow to learn (I certainly was). It doesn't matter how
good a job you can do on technical matters if you can't manage to
communicate them. Yes, writing skills do matter in technical
work; they matter quite a lot.

All of this is based on my personal obervation only. It isn't
any kind of official answer. You aren't going to get an official
answer to such a "why" question about the standard any more than
you are going to get an official answer as to why some election
came out with the result that it did.

John Harper

unread,
Sep 23, 2003, 8:13:39 PM9/23/03
to
In article <ued6drc...@altair.dfrc.nasa.gov>,

Richard Maine <nos...@see.signature> wrote:
>
>However, the expressed intention of a single person, even the person
>who authored the material, doesn't override the formal reviews and
>votes that passed the standard. There is no guarantee, just because
>that person says he had a different intent, that a majority would have
>agreed with that intent. A majority passed the words that are in the
>standard.

This sad story also applies in places other than Fortran. A year or two
ago our Parliament passed an Act allowing the police to arrest someone
if they had done A or B. (The precise definitions of A and B don't
matter for my present purpose.) The intention had been A and B. But the
law now says A or B. The mistake was pointed out while Parliament was
still discussing the matter, but party discipline was deemed more
important than giving the police excessive powers.

John Harper, School of Mathematical and Computing Sciences,
Victoria University, PO Box 600, Wellington, New Zealand
e-mail john....@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045

Dr Chaos

unread,
Sep 23, 2003, 10:40:44 PM9/23/03
to
Kamaraju Kusumanchi <kk...@cornell.edu> wrote:
> No I dont need to generate the same sequence of random numbers. The
> results cited in my first message seem to be very non-intuitive, they
> highlight another important aspect.
>
> One compiler changes the seed every time random_seed() is called. The
> other compiler hooks on to the same seed irrespective of howmany times
> the random_seed() subroutine is called. I felt this to be very
> confusing. Why is it so "loosely" defined? On a personal note, I cant
> see any advantage by leaving "what random_seed() should do" in an
> ambiguous state.
>
> By portability I meant the same functionality of random_seed() on
> different compilers.
> --raju

Here is what I do for ifc. It is a hack. I use it in
scripts and hence it prints to the standard output.

program iseedfromair
! make a random seed out of thin air, different
! almost every time.
! no guarantees on coverage or uniformity.
! 0 .. HUGE(0)
real :: r
integer :: seed, count, countrate, countmax
integer :: i, iters
integer :: maxr = HUGE(0) - 1

call random_seed
call system_clock(count,countrate,countmax)
iters = mod(count,101)
do i=1,iters+1
call random_number(r)
end do

seed = int(r*maxr)

print '(I0)', seed

end program iseedfromair

Michel OLAGNON

unread,
Sep 24, 2003, 2:50:56 AM9/24/03
to

In article <bkqb3c$h4v$1...@news01.cit.cornell.edu>, Kamaraju Kusumanchi <kk...@cornell.edu> writes:
>
>Thanks a lot. This code is working for me. I want to know one more
>thing. If I want to change the seed, what would be a good way to do it?
>This would be useful, for example, as I want to generate different
>random series for each run of the program. It is too boring to see the
>same sequence of numbers and calling it random! Any ideas on this also?
>I looked into section 8.16 of "Fortran 90/95 explained", but could not
>information on how to do this.
>
>For eg if size in the above program is 47. Is initializing each element
>of seed (from 1 to 47) to a value returned by date_and_time a good idea?
>or is there any other efficient method?
>
>Any code in this regard is very useful.
>regards
>raju
>

Many people use date_and_time.
I like another method: write down the last value of the seed
at the end of a run, and read it to use it as the initial seed for
the next run. I feel that it is safer than using date_and_time,
because I don't know how the seed digits are used, and the output of
date_and_time is strongly correlated between successive calls.

Michael Metcalf

unread,
Sep 24, 2003, 3:24:08 AM9/24/03
to

"Richard Maine" <nos...@see.signature> wrote in message
news:ued6drc...@altair.dfrc.nasa.gov...

> Kamaraju Kusumanchi <kk...@cornell.edu> writes:
>
> > Why is it [random seed] so "loosely" defined?
>
> Such "why" questions are notoriously difficult to answer. The
> "usual" and uninformative answer is "because that's the way the
> vote came out". That answer is all to apt in this case.
>

For what it's worth, here is a copy of my notes taken at the 1986 Halifax
meeting where a RNG was added to the standard:

"A persistent request has been for Fortran to provide a random number
intrinsic function. As was clear from the discussion, it was not intended
that the results of the function be portable, only the interface. There were
two proposals, both containg a means to extract the seed, one being a
function, the other a subroutine. The function was initially favoured
20-5-7. There was a long discussion on the effect a function has on
optimization and on side-effects -- the fact that the value of the seed is
changed at each call is a clear side-effect and therefore inappropriate for
a function. Additionally, a scalar value returned as part of an array-valued
expression would broadcast its single value over the whole expression, which
would certainly surprise the naive user. Sentiment turned against the
function, 8-11-3, and towards the subroutine, 13-3-14. (I voted against both
of them as I had increasing misgivings about the wisdom of standardizing on
something which is inherently non-portable without appearing to be so.)
After more furious debate, the subroutine was finally passed 15-12. My
motion to change the range from 0.0 <= x <= 1.0 to 0.0 <= x < 1.0 then
passed 28-0."

Regards,

Mike Metcalf


Dr Ivan D. Reid

unread,
Sep 24, 2003, 6:28:28 AM9/24/03
to
On Wed, 24 Sep 2003 09:24:08 +0200, Michael Metcalf
<michael...@t-online.de>
wrote in <bkrgp2$9pt$02$1...@news.t-online.com>:
> [My] motion to change the range from 0.0 <= x <= 1.0 to 0.0 <= x < 1.0 then
> passed 28-0."

I had a brief chat with Keith Ellis about this a few weeks ago.
Having ported an RNG around various systems, including machine code on
Z80s and 808xen, I'm not sure I can come up with a method that will generate
a uniform distribution including both *exactly* 0.0 and 1.0. Is this done
anywhere and if so, how?

--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan...@brunel.ac.uk Room 40-1-B12, CERN

David Jones

unread,
Sep 24, 2003, 6:59:05 AM9/24/03
to
Dr Ivan D. Reid wrote:
> On Wed, 24 Sep 2003 09:24:08 +0200, Michael Metcalf
> <michael...@t-online.de>
> wrote in <bkrgp2$9pt$02$1...@news.t-online.com>:
>> [My] motion to change the range from 0.0 <= x <= 1.0 to 0.0 <= x <
>> 1.0 then passed 28-0."
>
> I had a brief chat with Keith Ellis about this a few weeks ago.
> Having ported an RNG around various systems, including machine code
on
> Z80s and 808xen, I'm not sure I can come up with a method that will
> generate a uniform distribution including both *exactly* 0.0 and
1.0.
> Is this done anywhere and if so, how?

If you really want this (I've always needed both 0 and 1 excluded),
then a slow way is to take a generator giving 0.0 <= x < 1.0 , giving
U say, then return U with probability one half, otherwise return 1-U.
If you can get into the underlying code of a congruential generator
then you should be able to adjust the final divisor so that an exact 1
is a possible answer (e.g. divide by M-1 instead of M), provided that
the other parameters are set so that the maximum number of different
values are generated (i.e. M-1 is a possible value for the underlying
integer sequence).

David Jones


me...@skyway.usask.ca

unread,
Sep 24, 2003, 4:44:02 AM9/24/03
to
In a previous article, "Dr Ivan D. Reid" <Ivan...@brunel.ac.uk> wrote:
>On Wed, 24 Sep 2003 09:24:08 +0200, Michael Metcalf
> <michael...@t-online.de>
> wrote in <bkrgp2$9pt$02$1...@news.t-online.com>:
>> [My] motion to change the range from 0.0 <= x <= 1.0 to 0.0 <= x < 1.0 then
>> passed 28-0."
>
> I had a brief chat with Keith Ellis about this a few weeks ago.
>Having ported an RNG around various systems, including machine code on
>Z80s and 808xen, I'm not sure I can come up with a method that will generate
>a uniform distribution including both *exactly* 0.0 and 1.0. Is this done
>anywhere and if so, how?
>
>--
Does "exactly" have a meaning in reference to a floating point number?
Maybe that's the answer.
Chris

Richard Maine

unread,
Sep 24, 2003, 11:12:17 AM9/24/03
to
me...@skyway.usask.ca writes:

> Does "exactly" have a meaning in reference to a floating point number?

Yes.

It's meaning is not appropriate for many applications. It is often
not what is really wanted. Its use is subject to lots of caveats,
enough so that some compilers will give you a warning whenever you
use it.

But yes, it has a meaning, and it can be used reliably if you
understand the caveats. Testing for equality of calculated
values is likely to be problematic in many cases. But there are
many, many codes that quite reliably test for specific flag values.

In the case at hand, it can make quite a big difference whether
the exact value 1.0 is allowed or not. As in, it can easily
cause array bounds to be exceeded if you get a value of exactly
1.0 when you didn't expect that to be possible and failed to
protect against it. Makes a real good candidate for a subtle
bug that doesn't show up on most runs, but randomly (well,
pseudo-randomly) causes failures.

Dr Ivan D. Reid

unread,
Sep 24, 2003, 11:28:45 AM9/24/03
to
On 24 SEP 03 08:44:02 GMT, me...@skyway.usask.ca <me...@skyway.usask.ca>
wrote in <24SEP03....@skyway.usask.ca>:

> In a previous article, "Dr Ivan D. Reid" <Ivan...@brunel.ac.uk> wrote:

>>Having ported an RNG around various systems, including machine code on
>>Z80s and 808xen, I'm not sure I can come up with a method that will generate
>>a uniform distribution including both *exactly* 0.0 and 1.0. Is this done
>>anywhere and if so, how?

> Does "exactly" have a meaning in reference to a floating point number?


> Maybe that's the answer.

Both 0.0 and 1.0 are exactly representable in floating-point. My
method (and many others') of generating pseudo-random numbers generates
integers between 0 and 2^N-1, then scales into [0,1) by effectively
dividing by 2^N. Obviously this can generate 0.0 but never 1.0, and if I
generate the integers uniformly the floats will also be uniform -- there
is a philosophical question of whether to return denormalised numbers but
this only pertains when your probabilities are close to this
representation point. The suggestion of returning f or 1-f with 50%
probability looks OK at first glance until you realise that in fact 0.0
and 1.0 will each be returned with half the probability of any other valid
number. How much you care about this depends on the granularity (i.e.
1.0/2^N) compared to the statistics of the system being modelled, and whether
exactly 1.0 or 0.0 generates exceptional behaviour.

David Jones

unread,
Sep 24, 2003, 12:24:09 PM9/24/03
to
Dr Ivan D. Reid wrote:
> The suggestion of returning f or 1-f with 50%
> probability looks OK at first glance until you realise that in fact
> 0.0
> and 1.0 will each be returned with half the probability of any other
> valid number.

No ... I meant f or 1-f with 50% in every case, not just those
instances where f=0. However there may be some other asymmetry arising
from 1-(small) being returned as exactly 1, with a different behaviour
for 1-(close to 1).

David Jones


John C. Bollinger

unread,
Sep 24, 2003, 12:37:28 PM9/24/03
to
David Jones wrote:

> Dr Ivan D. Reid wrote:
>
>> The suggestion of returning f or 1-f with 50%
>>probability looks OK at first glance until you realise that in fact
>>0.0
>>and 1.0 will each be returned with half the probability of any other
>>valid number.
>
>
> No ... I meant f or 1-f with 50% in every case, not just those
> instances where f=0.

Ivan's comment still applies. Consider: any number other than 1.0 and
0.0 is a possible result whether the f or the 1-f case is chosen. The
extremes, on the other hand, are each possible in one case only,
therefore their probability is scaled by the probability that the
appropriate case is used.


John Bollinger
jobo...@idniana.edu

Kamaraju Kusumanchi

unread,
Sep 24, 2003, 11:34:29 PM9/24/03
to
Following Mike Metcalf's suggestion I wrote two simple subroutines for
randomizing the seed. I thought it would be useful if someone out there
wants final code to plugin into their programs. The routines are written
with ifc in mind, but can easily be changed for other compilers. It
assumes that there is a file "latestIfcSeed.dat" which contains the seed
to start with.


subroutine ifcSeedRead()
implicit none
!!Reads the seed from a file
integer, allocatable, dimension(:) :: array

integer :: i, no
call random_seed(size=no)
allocate(array(no))
open(UNIT=1, FILE="latestIfcSeed.dat", STATUS="OLD", ACTION="READ")
do i=1,no
read (1,*) array(i)
end do
call random_seed(put=array)
!write(*,*) array
if (allocated(array)) deallocate(array)
close(1)
end subroutine ifcSeedRead

subroutine ifcSeedWrite()
implicit none
!!Writes the seed to a file
integer, allocatable, dimension(:) :: array

integer :: i, no
call random_seed(SIZE=no)
allocate(array(no))
call random_seed(GET=array)

open(UNIT=1, FILE="latestIfcSeed.dat", STATUS="REPLACE", ACTION="WRITE")
do i=1,no
write(1,*) array(i)
end do
!write(*,*) array
if (allocated(array)) deallocate(array)
close(1)
end subroutine ifcSeedWrite

program test
implicit none
real, dimension(10) :: randnumber
real :: anumber

!previous seed is stored in latestseed.dat; read (*,*) it from there
!This initializes the seed
call ifcSeedRead()

call random_number(anumber)
write(*,*) anumber
call ifcSeedWrite()
end program test

PS: Even if one person finds this useful, I consider the bandwidth is
not wasted.
--raju

Dan Nagle

unread,
Sep 25, 2003, 8:09:46 AM9/25/03
to
Hello,

On Wed, 24 Sep 2003 23:34:29 -0400, Kamaraju Kusumanchi
<kk...@cornell.edu> wrote:

<snip>

> call random_seed(size=no)
> allocate(array(no))
> open(UNIT=1, FILE="latestIfcSeed.dat", STATUS="OLD", ACTION="READ")
> do i=1,no
> read (1,*) array(i)

<snip>

Good, you've got the seed written somewhere so it's available.
That a big help when debugging a particular run. :-)

Remember to save your "latestIfcSeed.dat" as part
of the documentation of a particular run. That helps
when a colleague asks how you got the results. :-)

--
Cheers!

Dan Nagle
Purple Sage Computing Solutions, Inc.

Dick Hendrickson

unread,
Sep 25, 2003, 11:45:19 AM9/25/03
to

Kamaraju Kusumanchi wrote:
>
> Following Mike Metcalf's suggestion I wrote two simple subroutines for
> randomizing the seed. I thought it would be useful if someone out there
> wants final code to plugin into their programs. The routines are written
> with ifc in mind, but can easily be changed for other compilers. It
> assumes that there is a file "latestIfcSeed.dat" which contains the seed
> to start with.

You can relax this requirement if you want to. Do an inquire by file,
INQUIRE ( FILE="latestIfcSeed.dat", exist = logical_variable)
before the OPEN.

If logical_variable comes back .FALSE. you can create the file and
put in some default seed.

Dick Hendrickson

Brooks Moses

unread,
Sep 26, 2003, 1:01:16 AM9/26/03
to
"Dr Ivan D. Reid" wrote:
> On Wed, 24 Sep 2003 09:24:08 +0200, Michael Metcalf
> <michael...@t-online.de>
> wrote in <bkrgp2$9pt$02$1...@news.t-online.com>:
> > [My] motion to change the range from 0.0 <= x <= 1.0 to 0.0 <= x < 1.0 then
> > passed 28-0."
>
> I had a brief chat with Keith Ellis about this a few weeks ago.
> Having ported an RNG around various systems, including machine code on
> Z80s and 808xen, I'm not sure I can come up with a method that will generate
> a uniform distribution including both *exactly* 0.0 and 1.0. Is this done
> anywhere and if so, how?

It's actually trivial to do this: using a slight modification to the
usual routines, create a random number in the range -1.0 < x <= 1.0.
(For instance, generate 0.0 <= x < 1.0, multiply by 2.0, subtract from
1.0.) If the x you get is negative, try again. Since these are
pseudorandom numbers, rather than random numbers, the number of times
you will have to try again before getting a nonnegative answer is
bounded, and so you can guarantee that the method will return a result
in some finitely bounded time.

Doing this in an efficient manner may be a bit more difficult, though.

- Brooks


--
Remove "-usenet" from my address to reply; the bmoses-usenet address
is currently disabled due to an overload of W32.Gibe-F worm emails.

Brooks Moses

unread,
Sep 26, 2003, 1:10:15 AM9/26/03
to
"Dr Ivan D. Reid" wrote:
> Both 0.0 and 1.0 are exactly representable in floating-point. My
> method (and many others') of generating pseudo-random numbers generates
> integers between 0 and 2^N-1, then scales into [0,1) by effectively
> dividing by 2^N. Obviously this can generate 0.0 but never 1.0, and if I
> generate the integers uniformly the floats will also be uniform -- there
> is a philosophical question of whether to return denormalised numbers but
> this only pertains when your probabilities are close to this
> representation point. The suggestion of returning f or 1-f with 50%
> probability looks OK at first glance until you realise that in fact 0.0
> and 1.0 will each be returned with half the probability of any other valid
> number.

I find it interesting that one could claim that (in many
implementations) those are the probabilities that should happen. For
instance, If you have random real numbers between zero and one
inclusive, and round them to the nearest n*2^(-k) for fixed k, the
probability of getting zero or one is half that of the other 2^k-1
numbers.

For a different challenge, write a program that will give an equivalent
result to taking a (pseudo)random real number between zero and one, and
returning the nearest floating-point representation. The increased
resolution as one gets near zero makes this notably more difficult, I'd
think....

Dr Ivan D. Reid

unread,
Sep 27, 2003, 8:37:43 AM9/27/03
to
On Thu, 25 Sep 2003 22:10:15 -0700, Brooks Moses
<bmoses...@cits1.stanford.edu>
wrote in <3F73CA37...@cits1.stanford.edu>:

> "Dr Ivan D. Reid" wrote:
>> Both 0.0 and 1.0 are exactly representable in floating-point. My
>> method (and many others') of generating pseudo-random numbers generates
>> integers between 0 and 2^N-1, then scales into [0,1) by effectively
>> dividing by 2^N. Obviously this can generate 0.0 but never 1.0, and if I
>> generate the integers uniformly the floats will also be uniform -- there
>> is a philosophical question of whether to return denormalised numbers but
>> this only pertains when your probabilities are close to this
>> representation point. The suggestion of returning f or 1-f with 50%
>> probability looks OK at first glance until you realise that in fact 0.0
>> and 1.0 will each be returned with half the probability of any other valid
>> number.

> I find it interesting that one could claim that (in many
> implementations) those are the probabilities that should happen. For
> instance, If you have random real numbers between zero and one
> inclusive, and round them to the nearest n*2^(-k) for fixed k, the
> probability of getting zero or one is half that of the other 2^k-1
> numbers.

Good point! I hadn't looked at it from that angle. I was fixated
on "in a uniform distribution all numbers are equally probable", but
should have remembered that each discrete number represents an _interval_.
So my desk-top histogramming 20 billion numbers an hour is just wasting
the weekend...

> For a different challenge, write a program that will give an equivalent
> result to taking a (pseudo)random real number between zero and one, and
> returning the nearest floating-point representation. The increased
> resolution as one gets near zero makes this notably more difficult, I'd
> think....

--

Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan...@brunel.ac.uk Room 40-1-B12, CERN

KotPT -- "for stupidity above and beyond the call of duty".

0 new messages