Home > Archive > Fortran > November 2005 > random_number
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
|
|
| Vladimir Gantovnik 2005-07-24, 8:54 pm |
| random_number() produces the uniform distribution within the range 0
<= x < 1.
How to get random number from the range 0 <= x <= 1, including 1?
Thank you!
| |
| meek@skyway.usask.ca 2005-07-24, 8:54 pm |
| In a previous article, "Vladimir Gantovnik" <gantovnik@yahoo.com> wrote:
>random_number() produces the uniform distribution within the range 0
><= x < 1.
>How to get random number from the range 0 <= x <= 1, including 1?
>Thank you!
Why does it matter? The probability of getting exact 1.0
is (theoretically) zero - or in practice, 1/ number_of_RVs.
Chris
| |
|
| DO
call random_number(x)
x=1.0001*x
if (x<=1.0) EXIT
ENDDO
and x the number you're after.
of course, the probability of getting 1.0 exactly is 0.0 (neglecting
machine precision issues) so this doesn't matter.
| |
| Richard E Maine 2005-07-24, 8:54 pm |
| In article <1122060491.002171.102850@z14g2000cwz.googlegroups.com>,
"Joost" <jv244@cam.ac.uk> wrote:
> DO
> call random_number(x)
> x=1.0001*x
> if (x<=1.0) EXIT
> ENDDO
> and x the number you're after.
>
> of course, the probability of getting 1.0 exactly is 0.0 (neglecting
> machine precision issues) so this doesn't matter.
And indeed, I'd argue that the above probably does more "damage" to the
randomness (of low-order bits) by the multiplication than it helps get
the requested range. Wouldn't surprise me at all to find that 1.0001* x
could never be exactly 1.0 and thus that this achieves nothing at all
(but I didn't check).
Elaborating on why having the "=" in x<=1.0 for the range shouldn't
matter anyway: It can be important in some applications that the result
is strictly <1.0. For example, it would be simple to show a practical
situation where a value of 1.0 could cause array bounds to be exceeded
or other similarly bad things to happen. However, it is hard for me to
imagine it being important that a value of exactly 1.0 could be returned.
--
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
| |
| Dr Ivan D. Reid 2005-07-24, 8:54 pm |
| On Fri, 22 Jul 2005 15:08:41 -0400, Vladimir Gantovnik <gantovnik@yahoo.com>
wrote in <dbrg8d$9bj$1@solaris.cc.vt.edu>:
> random_number() produces the uniform distribution within the range 0
><= x < 1.
> How to get random number from the range 0 <= x <= 1, including 1?
> Thank you!
There was a discussion on this (initiated by me) about September
2003. I suggest you use google groups to search for it -- I was finally
disabused of certain preconceptions by careful explanations from others.
--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan.Reid@[brunel.ac.uk|cern.ch] Room 40-1-B12, CERN
KotPT -- "for stupidity above and beyond the call of duty".
| |
| glen herrmannsfeldt 2005-07-24, 8:54 pm |
| Vladimir Gantovnik wrote:
> random_number() produces the uniform distribution within the range 0
> <= x < 1.
> How to get random number from the range 0 <= x <= 1, including 1?
RANDOM_NUMBER(X)
IF(X.EQ.0.5) X=1.0
I think you need to explain the problem in more detail.
Random number generators that I know of generate an integer between
0 and 2**N-1 for some N, and add the appropriate exponent to make
a floating point value.
For single precision it is likely N=24 and double N=53.
That is, one in 2**24 or one in 2**53 is the probability of
any specific value being returned. If you want one of the possible
random values to be 1.0 you will have to remove one of the others.
For symmetry I suggest 0.5, as above.
-- glen
| |
| Phillip Helbig---remove CLOTHES to reply 2005-07-24, 8:54 pm |
| In article <dbrg8d$9bj$1@solaris.cc.vt.edu>, "Vladimir Gantovnik"
<gantovnik@yahoo.com> writes:
> random_number() produces the uniform distribution within the range 0
> <= x < 1.
> How to get random number from the range 0 <= x <= 1, including 1?
> Thank you!
Find the maximum value produced by your RNG, and use the quotient of a
random number and that instead of the random number itself.
HOWEVER, why do you need to distinguish between 1 and the largest number
produced by your RNG?
| |
| Ingo Thies 2005-07-24, 8:54 pm |
| Hi,
gah@ugcs.caltech.edu (glen herrmannsfeldt) wrote:
> Random number generators that I know of generate an integer between
> 0 and 2**N-1 for some N, and add the appropriate exponent to make
> a floating point value.
>
> For single precision it is likely N=24 and double N=53.
> That is, one in 2**24 or one in 2**53 is the probability of
> any specific value being returned. If you want one of the possible
> random values to be 1.0 you will have to remove one of the others.
> For symmetry I suggest 0.5, as above.
I suppose it would be enough just to convert the integer to a real/double
precision value and divide it by 2**N-1 (instead of 2**N in the case of
exlusion of 1). This is what the real version of the Mersenne Twister
does (but it uses 32 as exponent, so its resolution is less of that of
double precision).
BTW: As far as I know, integers in F77 are restricted to the interval [-
2147483648,2147483647] = [-2**31,2**31-1]. Is there any simple
possibility to enhance the resolution of the Mersenne Twister on a normal
PC to fit that of double precision numbers? AFAIK are the double integers
(integer**8 and so on) only available on higher precision SUN computers
(which also can handle quad precision real numbers), aren't they?
Would it help just to replace integers by double precisions and use
aint(*) where it is necessary? Oder could I just add a downscaled random
perturbation to the first one, i.e. (grnd() is the MT function)
r1 = grnd()
r2 = grnd()/2**21
r = (r1+r2)*(2**21+1)/2**21
?
The exponent of 21 is just to get 2**53 from 2**32, and the division in
the last line to scale the sum to the interval [0,1].
Best regards
Ingo
| |
| Steven G. Kargl 2005-07-24, 8:54 pm |
| In article <9aTEQyLTe$B@thies-101284.user.dfncis.de>,
ingo.thies@gmx.de (Ingo Thies) writes:
>
> BTW: As far as I know, integers in F77 are restricted to the interval [-
> 2147483648,2147483647] = [-2**31,2**31-1].
AFAIK, this is incorrect. A quick scan of ansi-x3dot9-1978-Fortran77.pdf
shows that the Fortran 77 standard did not define the range of its
data types, not does it define the "model numbers".
> Is there any simple possibility to enhance the resolution of the
> Mersenne Twister on a normal PC to fit that of double precision
> numbers?
See the implementation of random_number in gfortran. The default
prng is KISS, but you can select to use MT. random_number is set
up to use the full precision of 24, 53, or 64 bits depending on
the kind type of the variable and the underlying OS.
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| Tim Prince 2005-07-24, 8:54 pm |
| Steven G. Kargl wrote:
> In article <9aTEQyLTe$B@thies-101284.user.dfncis.de>,
> ingo.thies@gmx.de (Ingo Thies) writes:
>
>
>
> AFAIK, this is incorrect. A quick scan of ansi-x3dot9-1978-Fortran77.pdf
> shows that the Fortran 77 standard did not define the range of its
> data types, not does it define the "model numbers".
Maybe what he means is that f77 does not have a standard way of
declaring larger integers (although most compilers have a non-standard
extension). The usual default size is as quoted; some compilers have a
way to make that larger by command line option, without requiring
non-f77 syntax.
>
>
>
>
> See the implementation of random_number in gfortran. The default
> prng is KISS, but you can select to use MT. random_number is set
> up to use the full precision of 24, 53, or 64 bits depending on
> the kind type of the variable and the underlying OS.
>
| |
| Gordon Sande 2005-07-24, 8:54 pm |
|
Tim Prince wrote:
> Steven G. Kargl wrote:
>
>
> Maybe what he means is that f77 does not have a standard way of
> declaring larger integers (although most compilers have a non-standard
> extension). The usual default size is as quoted; some compilers have a
> way to make that larger by command line option, without requiring
> non-f77 syntax.
A different reading is that Ingo Thies intended to say that many F77
implementations use a 4 byte or 32 bit ones complement integer
arithmetic but said that the F77 standard defined the integer arithmetic
that way. Steven G. Kargl then pointed out that the F77 standard did not
define the range of values for its integer data type.
There have been Fortrans on machines of other word sizes in the
past. 36 bits was "the standard" at one point in history. Some
machines even used sign magnitude (twos complement in other jargons)
for symmetric ranges.
This is common confusion between "the compiler I use" and "the
standard". There have been several threads over how small the range of
integers could be and still meet the standard.
[color=darkred]
>
| |
| Ron Shepard 2005-07-24, 8:54 pm |
| In article <D9REe.150995$tt5.54510@edtnps90>,
Gordon Sande <g.sande@worldnet.att.net> wrote:
> Some
> machines even used sign magnitude (twos complement in other jargons)
> for symmetric ranges.
I think the terminology is mixed here. See
http://en.wikipedia.org/wiki/Twos_Compliment, for example. I think
"sign magnitude" is the representation where there are two integer
zero representations and a symmetric range (like normal decimal
written notation), whereas "twos complement" has just a single zero
representation and has one more negative value than positive values.
In twos complement arithmetic, the zero maps to itself upon
negation, and the smallest negative number maps to itself upon
negation. I do not know how this terminology originated, it has
always seemed odd to me.
$.02 -Ron Shepard
| |
| Gordon Sande 2005-07-24, 8:54 pm |
|
Ron Shepard wrote:
> In article <D9REe.150995$tt5.54510@edtnps90>,
> Gordon Sande <g.sande@worldnet.att.net> wrote:
>
>
>
>
> I think the terminology is mixed here. See
> http://en.wikipedia.org/wiki/Twos_Compliment, for example. I think
> "sign magnitude" is the representation where there are two integer
> zero representations and a symmetric range (like normal decimal
> written notation), whereas "twos complement" has just a single zero
> representation and has one more negative value than positive values.
> In twos complement arithmetic, the zero maps to itself upon
> negation, and the smallest negative number maps to itself upon
> negation. I do not know how this terminology originated, it has
> always seemed odd to me.
>
Perfect description of what I thought I learned to call ones complement.
At best the terminology is strange although at one time I thought I
understood why it made technical sense.
One-of-the-complements has one zero and an asymmetric range and the
other-of-the-complements has two zeros (plus and minus) and a
symmetric range.
By the way I had used ones complement earlier as the asymmetric range so
you missed the chance to correct both errors, if I have the terminology
exactly backwards.
> $.02 -Ron Shepard
| |
| Ron Shepard 2005-07-24, 8:54 pm |
| In article <VVSEe.172059$on1.151760@clgrps13>,
Gordon Sande <g.sande@worldnet.att.net> wrote:
> Perfect description of what I thought I learned to call ones complement.
> At best the terminology is strange although at one time I thought I
> understood why it made technical sense.
I think "ones complement" is when negation consists of reversing all
the bits. I think CDC machines used this convention.
$.02 -Ron Shepard
| |
| glen herrmannsfeldt 2005-07-25, 4:25 am |
| Gordon Sande wrote:
(snip)
> There have been Fortrans on machines of other word sizes in the
> past. 36 bits was "the standard" at one point in history. Some
> machines even used sign magnitude (twos complement in other jargons)
> for symmetric ranges.
> This is common confusion between "the compiler I use" and "the
> standard". There have been several threads over how small the range of
> integers could be and still meet the standard.
I am pretty sure that Fortran still allows other than 32 bits, and
twos complement, ones complement, or sign magnitude, and probably
even allows decimal (sign magnitude, nines complement or tens
complement, or even more).
Historically, sign magnitude was used more, as it is more like
the way people do it on paper. The IBM 704 is 36 bit sign magnitude.
Like sign magnitude, ones complement is also symmetric with positive and
negative zero, but easier for doing addition and subtraction.
CDC made ones complement machines, and Univac still does.
A parallel adder is just about as easy to build in ones complement,
the difference is that the carry out from the MSB is used as a carry in
to the LSB.
C still allows twos complement, ones complement and sign magnitude for
signed arithmetic. Unsigned is unsigned. C does not allow decimal
representation unless it is made to look like binary.
-- glen
| |
| e p chandler 2005-07-25, 4:25 am |
| glen herrmannsfeldt wrote:
> Gordon Sande wrote:
> (snip)
>
>
>
[snip]
> Like sign magnitude, ones complement is also symmetric with positive and
> negative zero, but easier for doing addition and subtraction.
> CDC made ones complement machines, and Univac still does.
> A parallel adder is just about as easy to build in ones complement,
> the difference is that the carry out from the MSB is used as a carry in
> to the LSB.
>
[snip]
1. Just for fun, I spent some time this w end working with an old
8080 Fortran compiler (under CP/M emulation). It behaved as though it
had 48 bit words and packed binary coded decimal floating point
arithmetic. Very odd! (Integers implemented as reals, no character
variables, no double precision and lots of entertaining bugs. But
integers and reals did take the same amount of storage space - 6 eight
bit bytes!)
2. Back to the OP's question. He desired random numbers on (0,1] rather
than the usual [0,1). IMO, there are times when having strictly (0,1)
is desirable. I once was burned by a simulation program running on a
micro-computer which bombed for no apparent reason. Being naive, I
simply started the program with a different seed, and the problem went
away. I never thought I would be able to run through enough numbers on
a micro to encounter 0, or have to worry about taking the log of 0, but
i did :-(.
| |
| Ingo Thies 2005-07-25, 4:25 am |
| kargl@troutmask.apl.washington.edu (Steven G. Kargl) wrote:
> AFAIK, this is incorrect. A quick scan of ansi-x3dot9-1978-
Fortran77.pdf
> shows that the Fortran 77 standard did not define the range of its
> data types, not does it define the "model numbers".
Sorry, you are right. What I mean is, that in Fortran 77 there is no
standard way to change the accuracy, and on standard PCs this will result
in the precision I mentioned.
Thus I am looking for a simple way to enhance the MT routine to get the
full resolution without destroying the random distribution due to
numerical problems.
Best regards,
Ingo
| |
|
|
| Ingo Thies 2005-07-25, 5:34 pm |
| Hi,
> Hi, I use this version of MT (
> www.lce.hut.fi/research/polymer/softsimu2002/mt19937.f ), and it has as
> output a number in [0,1].
I use the same version, and it simply converts a 4 byte integer (the
minimum and maximum values are explicitely given in the declaration) into
a double precision number. Therefore the minimum possible difference
between two numbers is 2^-32 = 2.3*10^-10 which is well above the
accuracy limit of double precision (app. 10^-16).
Best regards,
Ingo
| |
| glen herrmannsfeldt 2005-07-25, 5:34 pm |
| Ingo Thies wrote:
(snip)
> I use the same version, and it simply converts a 4 byte integer (the
> minimum and maximum values are explicitely given in the declaration) into
> a double precision number. Therefore the minimum possible difference
> between two numbers is 2^-32 = 2.3*10^-10 which is well above the
> accuracy limit of double precision (app. 10^-16).
A C library routine, I am not sure it is part of the standard, but it
is in most unix implementations, called rand48, has a 48 bit linear
congruential generator. The one returning a double, drand48, I believe
uses the 48 bit value for the double. A little less than the usual 53
bit fraction, but not so much less.
-- glen
| |
| Ingo Thies 2005-07-25, 5:34 pm |
| Hi glen,
at 25 Jul 05 you wrote:
> A C library routine, I am not sure it is part of the standard, but it
> is in most unix implementations, called rand48, has a 48 bit linear
> congruential generator. The one returning a double, drand48, I believe
> uses the 48 bit value for the double. A little less than the usual 53
> bit fraction, but not so much less.
Again my question: Would it be a suitable solution to add a small random
perturbation to the 32 bit integer based number as I suggestet before?
Just to remember, I thought about something like this
r1 = grnd()
r2 = grnd()/2**21
r = (r1+r2)*(2**21+1)/2**21
r1 is just the course standard output of the Mersenne Twister function,
and r2 is a second one. r2 is scaled down to fill the spaces the 32 bit
representation leaves. Mathematically, the final random number r should
now be randomly distributed within [0,1) with a 53 bit resolution. But I
am not shure about numerical problems (especially rounding errors) which
could interfere here. I am not experienced enough in testing such a self-
fudged random number generator in a scientifically satisfying way...
Best regards,
Ingo
| |
| Steven G. Kargl 2005-07-25, 5:34 pm |
| In article <9aXMHbdDe$B@thies-101284.user.dfncis.de>,
ingo.thies@gmx.de (Ingo Thies) writes:
> Again my question: Would it be a suitable solution to add a small random
> perturbation to the 32 bit integer based number as I suggestet before?
> Just to remember, I thought about something like this
>
> r1 = grnd()
> r2 = grnd()/2**21
> r = (r1+r2)*(2**21+1)/2**21
>
> r1 is just the course standard output of the Mersenne Twister function,
> and r2 is a second one. r2 is scaled down to fill the spaces the 32 bit
> representation leaves. Mathematically, the final random number r should
> now be randomly distributed within [0,1) with a 53 bit resolution. But I
> am not shure about numerical problems (especially rounding errors) which
> could interfere here. I am not experienced enough in testing such a self-
> fudged random number generator in a scientifically satisfying way...
Your last 2 sentence sum up why you should not do the above.
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| glen herrmannsfeldt 2005-07-25, 5:34 pm |
| Ingo Thies wrote:
(snip regarding the C rand48 routines)
> Again my question: Would it be a suitable solution to add a small random
> perturbation to the 32 bit integer based number as I suggestet before?
> Just to remember, I thought about something like this
> r1 = grnd()
> r2 = grnd()/2**21
> r = (r1+r2)*(2**21+1)/2**21
> r1 is just the course standard output of the Mersenne Twister function,
> and r2 is a second one. r2 is scaled down to fill the spaces the 32 bit
> representation leaves. Mathematically, the final random number r should
> now be randomly distributed within [0,1) with a 53 bit resolution. But I
> am not shure about numerical problems (especially rounding errors) which
> could interfere here. I am not experienced enough in testing such a self-
> fudged random number generator in a scientifically satisfying way...
A good 32 bit random number generator has a period of 2**32. A poor
one has even less, though Knuth explains in volume 2 how to select a
good one.
If you take two successive numbers from a series with a period of 2**32
you now have a series of 64 bit numbers with a period of 2**31.
Personally, I prefer to use integer arithmetic as much as possible in
doing random number problems, but it depends much on what you are
actually trying to do.
As far as:
r1 = grnd()
r2 = grnd()/2**21
r = (r1+r2)*(2**21+1)/2**21
The result is very sensitive to any rounding done by the floating
point arithmetic. If I understand what you mean, you actually meant
r = (r1+r2)*(2**21)/(2**21+1)
but consider what
r = (r1+r2)*(2**21+1)/2**21
does. It takes (r1+r2) and adds (r1+r2) shifted right 21 bits.
That doesn't do much to improve the randomness, and may or may not
result in the desired range.
If you want to combine two random number generators to generate a larger
random number you need two separate generators with incommensurate
periods. That is, the GCD of the periods should be 1.
At that point I believe it is necessary to know in much more detail what
you are trying to do with the numbers. If you are at the point where a
period of 2**32 isn't enough, if you are generating a significant
fraction of 2**32 numbers, there are many other details to worry about.
-- glen
| |
| James Giles 2005-07-25, 10:09 pm |
| glen herrmannsfeldt wrote:
....
> A good 32 bit random number generator has a period of 2**32. A poor
> one has even less, though Knuth explains in volume 2 how to select a
> good one.
Actually, there are cycle free random number generators. And, many
generators have longer cycles, though still finite. All these require that
the seed is not the same as one of the samples. The gererator I use
most has a cycle of 2^64 - 2^32. It's the modulo 2^32 sum of a linear
congruential generator and a shift register sequence. The first has a
cycle of 2^32 and the second has a cycle of 2^32-1. Since they're
independent, the cycle of their sum (or XOR) is the product of their
separate cycles.
If I need double precision samples, I usually use a pair (again a linear
congruential and a shift register) that are each 64-bit. This has a cycle
of 2^128 - 2^64.
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
| |
| glen herrmannsfeldt 2005-07-25, 10:09 pm |
| James Giles wrote:
> glen herrmannsfeldt wrote:
[color=darkred]
> Actually, there are cycle free random number generators. And, many
> generators have longer cycles, though still finite. All these require that
> the seed is not the same as one of the samples.
OK. Well, the cycle can't be longer than the possible combinations of
the state register. So, yes, when I said 32 bit I meant ones with 32
bits of state, not necessarily ones that returned only 32 bits.
If you put in enough bits it won't be possible to see the cycle in less
than the lifetime of the universe. Some Xilinx chips have built in
logic for LFSR (linear feedback shift registers). It isn't hard to
combine them for a really long register, thousands or even millions of
bits. Finding good coefficients for such a register might be hard,
though, but there would be no worry for the cycle to ever repeat.
> The gererator I use
> most has a cycle of 2^64 - 2^32. It's the modulo 2^32 sum of a linear
> congruential generator and a shift register sequence. The first has a
> cycle of 2^32 and the second has a cycle of 2^32-1. Since they're
> independent, the cycle of their sum (or XOR) is the product of their
> separate cycles.
I would have called that 64 bits, as it has 64 bits of state registers,
even though it only returns 32 bits.
> If I need double precision samples, I usually use a pair (again a linear
> congruential and a shift register) that are each 64-bit. This has a cycle
> of 2^128 - 2^64.
-- glen
| |
| Gordon Sande 2005-07-26, 9:05 am |
|
glen herrmannsfeldt wrote:
> Ingo Thies wrote:
>
> (snip regarding the C rand48 routines)
>
>
>
>
>
>
>
> A good 32 bit random number generator has a period of 2**32. A poor
> one has even less, though Knuth explains in volume 2 how to select a
> good one.
There is clearly a terminology problem here. If the recurrence is of
order 1 like the typical multiplicative congruence pseudo random
number generator then the statement stands. The internal state (seed)
will have 2**32 vales. BUT if the recurrence is of higher order then
this is trivially wrong. For a second order recurrence the period would
be (2**32)**2 {neglecting a variety of technical fine points) and the
internal state would have 2**64 values, or be a two word seed in plain
language. And so on for higher order recurrences. The modern long
period PRNGs have rather high order recurrences and rather large
number of words of internal state. Knuth has some example of order 2
recurrences to illustrate the point. The Mersenne Twister uses slightly
different number theory as it is a GF(2**16000+) rather than a GF(p**k)
with p around 2**30 and more moderate k. The details are the sort of
thing that only a number theorist could love as it as it is technical
and arcane. (It is a fun read if you like that sort of thing.) The
quick and dirty is that the period length matches the seed size and
can be much larger than the reported out value as shown by the multiple
word seeds of the long period generators.
A higher order recurrence will be exactly equi-distributed to its order
with MT being exact to in excess of 600 dimensions, if I recall the
specs correctly. So the classical order one PRNG is exactly uniform
and can be pretty good up to six or so dimensions. A second order PRNG
will be exactly uniform on a square and can be pretty good through more.
The interest in very long periods was driven by both curiosity and the
practical need to shuffle cards (52 factorial is a big number) correctly
for the gaming folks.
Knuth is a good introduction but there is quite a bit more available.
Random number generation and cryptography share a fair bit of ground.
> If you take two successive numbers from a series with a period of 2**32
> you now have a series of 64 bit numbers with a period of 2**31.
>
> Personally, I prefer to use integer arithmetic as much as possible in
> doing random number problems, but it depends much on what you are
> actually trying to do.
>
> As far as:
>
> r1 = grnd()
> r2 = grnd()/2**21
> r = (r1+r2)*(2**21+1)/2**21
>
> The result is very sensitive to any rounding done by the floating
> point arithmetic. If I understand what you mean, you actually meant
>
> r = (r1+r2)*(2**21)/(2**21+1)
>
> but consider what
>
> r = (r1+r2)*(2**21+1)/2**21
>
> does. It takes (r1+r2) and adds (r1+r2) shifted right 21 bits.
> That doesn't do much to improve the randomness, and may or may not
> result in the desired range.
>
> If you want to combine two random number generators to generate a larger
> random number you need two separate generators with incommensurate
> periods. That is, the GCD of the periods should be 1.
>
> At that point I believe it is necessary to know in much more detail what
> you are trying to do with the numbers. If you are at the point where a
> period of 2**32 isn't enough, if you are generating a significant
> fraction of 2**32 numbers, there are many other details to worry about.
>
> -- glen
>
| |
| Ingo Thies 2005-07-26, 5:03 pm |
| glen, at 25 Jul 05 you wrote to :
> A good 32 bit random number generator has a period of 2**32. A poor
> one has even less, though Knuth explains in volume 2 how to select a
> good one.
MT19937 should have a period of 2**19937-1
<http://www.math.sci.hiroshima-u.ac....ARTICLES/mt.pdf>
This should not be the problem... ;-)
> Personally, I prefer to use integer arithmetic as much as possible in
> doing random number problems, but it depends much on what you are
> actually trying to do.
Normally I need a random generator for Monte Carlo experiments (e.g.
random samples of stars mapped with the stellar mass function of the
cluster) and for random distributions of particles (e.g. in N-body
simulations). But other purposes are possible (even if I don't know them
yet).
> r1 = grnd()
> r2 = grnd()/2**21
> r = (r1+r2)*(2**21+1)/2**21
>
> The result is very sensitive to any rounding done by the floating
> point arithmetic. If I understand what you mean, you actually meant
>
> r = (r1+r2)*(2**21)/(2**21+1)
Of course, now it is correct.
> does. It takes (r1+r2) and adds (r1+r2) shifted right 21 bits.
> That doesn't do much to improve the randomness, and may or may not
> result in the desired range.
The idea is to fill the holes the coarser 32 bit distribution leaves.
Without this added second random number the smallest possible difference
between two numbers is about 10^-10 compared to 10^-16 in double
precision floating point. The rounding errors should be in the 10^-16
order, so I suppose that there should be an improvement if many random
numbers are needed. However, the improvement may be somewhat less then
the full 16 decimals, maybe 15 or ven 14 decs, but this would still be a
much higher resolution (by at least four decimals) than the original
resolution, if I didn't forget anything important.
The best method, however, would be to enlarge the possible range of
integer number.
-- Ingo
| |
| Dr Ivan D. Reid 2005-07-27, 5:06 pm |
| On 24 Jul 2005 23:17:15 -0700, e p chandler <epc8@juno.com>
wrote in <1122272235.189087.172890@f14g2000cwb.googlegroups.com>:
> 2. Back to the OP's question. He desired random numbers on (0,1] rather
> than the usual [0,1).
No, he wanted it on [0,1]; (0,1] is trivially derivable from [0,1).
> IMO, there are times when having strictly (0,1)
> is desirable.
Not usually if you look at the underlying physics.
> I once was burned by a simulation program running on a
> micro-computer which bombed for no apparent reason. Being naive, I
> simply started the program with a different seed, and the problem went
> away. I never thought I would be able to run through enough numbers on
> a micro to encounter 0, or have to worry about taking the log of 0, but
> i did :-(.
You and me both, but as I said it's trivial to invert the RNG
so that you never take log(0) but do take log(1).
I originally pointed at a discussion on the topic almost two years
ago. The URL is (hopefully):
http://groups-beta.google.com/group...98386aeb1ea0082
--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan.Reid@[brunel.ac.uk|cern.ch] Room 40-1-B12, CERN
GSX600F, RG250WD, DT175MX "You Porsche. Me pass!" DoD #484 JKLO# 003, 005
WP7# 3000 LC Unit #2368 (tinlc) UKMC#00009 BOTAFOT#16 UKRMMA#7 (Hon)
KotPT -- "for stupidity above and beyond the call of duty".
| |
|
| > HOWEVER, why do you need to distinguish between 1 and the largest
number
> produced by your RNG?
I am working with genetic algorithm with continuous variables, the
range of changing variables is [0.0,1.0], and for my test problem I
know that exact optimum solution is x=(1.0,1.0,...,1.0). I understand
that the solution x=(0.9999,0.9999,.....) will be ok, but...
| |
| Random Programmer 2005-07-31, 5:17 pm |
| Hmm, how about:
--------------------
10000 continue
Call rannum(random,idum)
output = random*1.1
output = output - 0.05
if (output < 0.0) then
goto 10000
endif
if (output > 1.0) then
goto 10000
endif
----------------------------
Assuming "rannum" is your random number generation subroutine,
otherwise use Fortran's implicit function.
Vladimir, would you happen to be the same guy who wrote the 2002 paper
with Gurdal & Watson?
| |
| James Giles 2005-07-31, 5:17 pm |
| I've not had time to address this before, but the problem
of adjusting the range or the random generator from [0.0, 1.0)
(half-open interval) so that you get the range [0.0, 1.0] (closed
interval) is ill-defined. For the half-open interval, uniform
distribution of possible return values seems quite appropriate
(especially given the usual applications that use it). But, if
you want the closed interval, I think that decision needs to be
reconsidered. If you think of each possible result as being
the center of an interval of equally probable REAL outcomes,
the intervals containing 0.0 and 1.0 are half as large as the
intervals surrounding each other possibility.
What about generating a random number as usual [0.0,1.0).
Then, if the result is exactly 0.0, you generate independently
a single random bit such that if this bit is one you return 1.0
and if it's zero you return 0.0. Now you have the appropriate
probabilities for these end points.
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
| |
|
| > Vladimir, would you happen to be the same guy who wrote the 2002
paper
> with Gurdal & Watson?
Thank you! Yes, I did it! :)
| |
| David Jones 2005-08-01, 4:06 am |
| VBG wrote:
> I am working with genetic algorithm with continuous variables, the
> range of changing variables is [0.0,1.0], and for my test problem I
> know that exact optimum solution is x=(1.0,1.0,...,1.0). I
understand
> that the solution x=(0.9999,0.9999,.....) will be ok, but...
Sounds as if you really need to move to a model for your genetic
transitions that copes with what you want to do more explicitly and
controlably. Thus, have a state model that distinguishes explicitly,
between exacly zero, exactly one and intermediate and define
conditional transition probabilities in a reasonable way:
(i) probability of staying in exactly zero or exactly one moderately
high, possibly increasing as the number of iterations increases.
(ii) probability of moving to exactly zero or exactly one from an
intermediate value increasing with closeness to zero or one.
Of course the problem is then not one of programming but of designing
a suitable model for the transitions.
David Jones
| |
|
| >In a previous article, "Vladimir Gantovnik" <gantovnik@yahoo.com> wrote:
>random_number() produces the uniform distribution within the range 0
><= x < 1.
>How to get random number from the range 0 <= x <= 1, including 1?
>Thank you!
You could try multiplying your random number by
NEAREST(1.0, 2.0)
and, just for safety, checking that the product is
not greater than 1.0. If it is, replace it with 1.0.
| |
| Rich Townsend 2005-11-15, 7:01 pm |
| robin wrote:
>
>
> You could try multiplying your random number by
> NEAREST(1.0, 2.0)
> and, just for safety, checking that the product is
> not greater than 1.0. If it is, replace it with 1.0.
>
>
No, this is wrong. If you are mapping numbers > 1.0 back to 1.0, then your
distribution is not going to be uniform over 0. <= x <= 1.
cheers,
Rich
| |
| Dr Ivan D. Reid 2005-11-15, 7:01 pm |
| On Tue, 15 Nov 2005 14:11:23 GMT, robin <robin_v@bigpond.com>
wrote in <f0mef.18282$Hj2.7323@news-server.bigpond.net.au>:[color=darkred]
We had this discussion a wee while ago, and some reasonable
solutions were proposed. However, I'm still not convinced whether many
applications require both endpoints to be included[1] (the one that
started the discussion was a HEP problem that came up in a summer-school
I attended in 2003). What's the application?
[1] E.g. often a logarithm is involved so one endpoint cannot be exact!
--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan.Reid@[brunel.ac.uk|cern.ch] Room 40-1-B12, CERN
| |
|
| Rich Townsend wrote in message ...
>robin wrote:
>
>No, this is wrong.
It is not wrong.
> If you are mapping numbers > 1.0 back to 1.0, then your
>distribution is not going to be uniform over 0. <= x <= 1.
It does NOT map numbers back into 1.0.
Numerically, it cannot map any number <1.0 to a number
greater than 1.0.
But, just in case a number should be rounded up to just
over 1.0, the IF will reset it to 1.0.
>cheers,
>
>Rich
| |
| Rich Townsend 2005-11-15, 9:57 pm |
| robin wrote:
> Rich Townsend wrote in message ...
>
>
>
> It is not wrong.
>
>
>
>
> It does NOT map numbers back into 1.0.
> Numerically, it cannot map any number <1.0 to a number
> greater than 1.0.
So a number <1.0 cannot end up >1.0
> But, just in case a number should be rounded up to just
> over 1.0, the IF will reset it to 1.0.
>
So a number <1.0 can end up >1.0
What, is this Schrodinger's Logic?
cheers,
Rich
| |
|
|
<22JUL05.13240739@skyway.usask.ca>...
In a previous article, "Vladimir Gantovnik" <gantovnik@yahoo.com> wrote:
>random_number() produces the uniform distribution within the range 0
><= x < 1.
>How to get random number from the range 0 <= x <= 1, including 1?
>Thank you!
One further point is that if using NEAREST,
the kind of 1.0 and 2.0 used will need to be the same
as that of your random number.
| |
| Rich Townsend 2005-11-16, 7:02 pm |
| robin wrote:
> <22JUL05.13240739@skyway.usask.ca>...
> In a previous article, "Vladimir Gantovnik" <gantovnik@yahoo.com> wrote:
>
>
>
> One further point is that if using NEAREST,
> the kind of 1.0 and 2.0 used will need to be the same
> as that of your random number.
>
>
Why can't you use NEAREST(1., 1.)? Why does the second argument have to be 2.?
| |
| Rich Townsend 2005-11-16, 7:02 pm |
| robin wrote:
> <22JUL05.13240739@skyway.usask.ca>...
> In a previous article, "Vladimir Gantovnik" <gantovnik@yahoo.com> wrote:
>
>
>
> One further point is that if using NEAREST,
> the kind of 1.0 and 2.0 used will need to be the same
> as that of your random number.
Sure. But your approach is still fundamentally flawed and broken -- for exactly
the reason that I have already stated. Take the following test program:
program foo
implicit none
real(KIND(0.D0)) :: a,b,scale
a = NEAREST(1.D0, -1.D0)
b = NEAREST(a, -1.D0)
scale = NEAREST(1.D0, 1.D0)
print *,'a < 1? : ', a<1.D0
print *,'b < 1? : ', b<1.D0
print *,'a == b? : ', a==b
print *,'a maps to 1? : ', a*scale == 1.D0
print *,'b maps to 1? : ', b*scale == 1.D0
end program foo
I've tried this program on both Lahey lf95 Express and Intel Fortran (Linux;
x86). I get the following output:
a < 1? : T
b < 1? : T
a == b? : F
a maps to 1? : T
b maps to 1? : T
So, we have two *distinct* numbers, both less than unity, that map to unity
using your approach.
As a consequence, your naive attempt at generating a random number 0 <= x <= 1
produces a non-uniform distribution, and is therefore wholly worthless to the OP.
cheers,
Rich
| |
| Rich Townsend 2005-11-16, 7:02 pm |
| Rich Townsend wrote:
> robin wrote:
>
>
>
> Sure. But your approach is still fundamentally flawed and broken -- for
> exactly the reason that I have already stated. Take the following test
> program:
>
> program foo
>
> implicit none
>
> real(KIND(0.D0)) :: a,b,scale
>
> a = NEAREST(1.D0, -1.D0)
> b = NEAREST(a, -1.D0)
>
> scale = NEAREST(1.D0, 1.D0)
>
> print *,'a < 1? : ', a<1.D0
> print *,'b < 1? : ', b<1.D0
>
> print *,'a == b? : ', a==b
>
> print *,'a maps to 1? : ', a*scale == 1.D0
> print *,'b maps to 1? : ', b*scale == 1.D0
>
> end program foo
>
> I've tried this program on both Lahey lf95 Express and Intel Fortran
> (Linux; x86). I get the following output:
>
> a < 1? : T
> b < 1? : T
> a == b? : F
> a maps to 1? : T
> b maps to 1? : T
>
> So, we have two *distinct* numbers, both less than unity, that map to
> unity using your approach.
>
> As a consequence, your naive attempt at generating a random number 0 <=
> x <= 1 produces a non-uniform distribution, and is therefore wholly
> worthless to the OP.
Unless, of course, you can guarantee that the random number generator will only
return one of a or b. Can such guarantees be made? How?
cheers,
Rich
| |
|
| Rich Townsend wrote in message ...
>robin wrote:
>So a number <1.0 cannot end up >1.0
>So a number <1.0 can end up >1.0
>
>What, is this Schrodinger's Logic?
No; I thought tht I made it perfectly clear for you.
Here is is again:
"But, just in case a number should be rounded up to just
"over 1.0, the IF will reset it to 1.0."
Never forget Robert's Law.
>cheers,
>
>Rich
| |
|
| Rich Townsend wrote in message ...
>robin wrote:
>
>Why can't you use NEAREST(1., 1.)? Why does the second argument have to be 2.?
The manual will tell you. (Hint: Where is the second 1.0 going?)
| |
|
| Rich Townsend wrote in message ...
>robin wrote:
>
>Sure. But your approach is still fundamentally flawed and broken -- for exactly
>the reason that I have already stated. Take the following test program:
>
>program foo
>
> implicit none
>
> real(KIND(0.D0)) :: a,b,scale
>
> a = NEAREST(1.D0, -1.D0)
> b = NEAREST(a, -1.D0)
>
> scale = NEAREST(1.D0, 1.D0)
>
> print *,'a < 1? : ', a<1.D0
> print *,'b < 1? : ', b<1.D0
>
> print *,'a == b? : ', a==b
>
> print *,'a maps to 1? : ', a*scale == 1.D0
> print *,'b maps to 1? : ', b*scale == 1.D0
>
>end program foo
>
>I've tried this program on both Lahey lf95 Express and Intel Fortran (Linux;
>x86). I get the following output:
>
> a < 1? : T
> b < 1? : T
> a == b? : F
> a maps to 1? : T
> b maps to 1? : T
>
>So, we have two *distinct* numbers, both less than unity, that map to unity
>using your approach.
What value for scale do you think you have?
>As a consequence, your naive attempt at generating a random number 0 <= x <= 1
>produces a non-uniform distribution, and is therefore wholly worthless to the OP.
Apart from 'scale' just mentioned, you may have to specify
a compiler option to prevent the f.p. processor from rounding.
And as for non-uniformity, you haven't demonstrated anything,
and as for 'worthless to the OP', you don't know that either.
>cheers,
>
>Rich
| |
| James Giles 2005-11-16, 9:58 pm |
| robin wrote:
> Rich Townsend wrote in message ...
....
>
> The manual will tell you. (Hint: Where is the second 1.0 going?)
Once again, Robin doesn't read the standard but elects to
lecture others on it's content. NEAREST produces the next
representable neighbor of its first argument in the direction
of the infinity that has the same sign as it's second argument.
So, NEAREST(1.0, DIR) is the same value regardless
of the magnitude of the value DIR (provided DIR is not
equal to zero, which is unecessarily prohibited by the
standard - probably to accomodate whatever implementations
remain that don't have signed zeros). The return value of
NEAREST only depends on the *sign* of the second
argument.
NEAREST(1.0, 10000000.0) == NEAREST(1.0,0.00000001)
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
| |
| Rich Townsend 2005-11-16, 9:58 pm |
| robin wrote:
> Rich Townsend wrote in message ...
>
>
>
> The manual will tell you. (Hint: Where is the second 1.0 going?)
>
>
From the Fortran 95 standard (ISO/IEC 1539-1997):
----
13.14.76 NEAREST (X, S)
Description. Returns the nearest different machine representable number in a
given direction.
Class. Elemental function.
Arguments.
X shall be of type real.
S shall be of type real and not equal to zero.
Result Characteristics. Same as X.
Result Value. The result has a value equal to the machine representable number
distinct from X and nearest to it in the direction of the infinity with the same
sign as S.
Example. NEAREST (3.0, 2.0) has the value on a machine whose representation is
that of the model at the end of 13.7.1.
----
Hence, there is absolutely no difference between NEAREST(1., 1.) and NEAREST(1.,
2.). So do you want to have another bite at explaining why you use 2. rather
than 1. as the second argument?
cheers,
Rich
| |
| Rich Townsend 2005-11-16, 9:58 pm |
| robin wrote:
> Rich Townsend wrote in message ...
>
>
>
> What value for scale do you think you have?
>
Why, the same as NEAREST(1.D0, 2.D0), of course!
>
>
>
> Apart from 'scale' just mentioned, you may have to specify
> a compiler option to prevent the f.p. processor from rounding.
>
So your procedure for generating random numbers on 0 <= x <= 1 now relies on
compiler options? Hmmm, doesn't sound very portable to me.
cheers,
Rich
| |
|
| James Giles wrote in message <4SRef.113941$zb5.89578@bgtnsc04-news.ops.worldnet.att.net>...
>robin wrote:
>...
>
>Once again, Robin doesn't read the standard but elects to
>lecture others on it's content.
Don't talk rot.
You really are jumping to conclusions.
| |
|
| Rich Townsend wrote in message ...
>robin wrote:
>
>Why, the same as NEAREST(1.D0, 2.D0), of course!
That's right.
>
>So your procedure for generating random numbers on 0 <= x <= 1 now relies on
>compiler options? Hmmm, doesn't sound very portable to me.
It might be more worthwhile your investigating the reason
that your results are inexact.
>cheers,
>
>Rich
| |
|
| Rich Townsend wrote in message ...
>robin wrote:
>
>Sure. But your approach is still fundamentally flawed and broken --
There's nothing "fundamentally flawed" about the method,
as the following proves. Below are printed the first ten values below 1.0
(left-hand column) and their products (right-hand column) with
the nearest value above 1.0 :- Only the least-significant 8 bits of each
mantissa are printed.
11111111 00000000
11111110 11111111
11111101 11111110
11111100 11111101
11111011 11111100
11111010 11111011
11111001 11111010
11111000 11111001
11110111 11111000
>So, we have two *distinct* numbers, both less than unity, that map to unity
>using your approach.
Not here, nor in general.
>cheers,
>
>Rich
| |
|
| Rich Townsend wrote in message ...
>As a consequence, your naive attempt at generating a random number 0 <= x <= 1
Did I attempt? I'm sure that you will agree that I didn't.
I SUGGESTED that the OP try something. Here it is again, in case that you missed it:-
"You could try multiplying your random number by
"NEAREST(1.0, 2.0)
"and, just for safety, checking that the product is
"not greater than 1.0. If it is, replace it with 1.0."
>produces a non-uniform distribution, and is therefore wholly worthless to the OP.
Let the OP decide that. It is not your call.
>cheers,
>Rich
| |
| glen herrmannsfeldt 2005-11-17, 9:57 pm |
| robin wrote:
(snip)
> "You could try multiplying your random number by
> "NEAREST(1.0, 2.0)
> "and, just for safety, checking that the product is
> "not greater than 1.0. If it is, replace it with 1.0."
if(x.eq.0.5) x=1.0
is somewhat easier, and has the same result.
There are only a finite number of return values from the random number
generator. You can move them around, but you can't increase the number
by applying mathematical operations.
Many of the suggested methods give the 1.0 result while removing the
possible 0.5 result.
-- glen
| |
| James Giles 2005-11-18, 3:57 am |
| robin wrote:
....
> There's nothing "fundamentally flawed" about the method,
> as the following proves. Below are printed the first ten values
> below 1.0 (left-hand column) and their products (right-hand column)
> with
> the nearest value above 1.0 :- Only the least-significant 8 bits of
> each mantissa are printed.
>
> 11111111 00000000
> 11111110 11111111
> 11111101 11111110
> 11111100 11111101
> 11111011 11111100
> 11111010 11111011
> 11111001 11111010
> 11111000 11111001
> 11110111 11111000
Most random number generators return values in the
interval 0 <= r < 1 that are uniformly spaced in that
interval and are equiprobable. Both of those last properties
are very desirable. Your "solution" will not return values
on the interval 0 <= r <= 1 that are uniformly spaced.
For example, on IEEE single precision, you get 2^24
values in the interval 0 <= r < 1 and the spacing between
the possible results is 1/2^24. There are (2^24)+1
values with that spacing in the interval 0 <= r <= 1.
Most of the values generated by your "solution"
will have that spacing. But there will be at least one such
value in that interval which your algorithm will *not*
generate. (How could it? There are only 2^24 values
generated to cover (2^24)+1 positions.) Values not
equat to an integer multiple of 1/2^24 will probably be
generated as well.
There are several really useful methods to generate numbers
on the interval 0 <= r <= 1. One is to generate numbers
in the usual way, reject any that are greater than 0.5 (roll
again) and double the result. This gives uniform spacing and
equiprobable results.
Another is to have an independent generator that creates either
zero or one. Generate a number in the usual way (the interval
is 0 <= r < 1), if the value is exactly zero, call the second
generator, the second value is zero return 0, otherwise return 1.
All values between zero and one have the same probability,
zero and one have half that probability. This is appropriate
if the samples are viewed as the mid-points of a continuous
distribution between zero and one (and the endpoints
correspond to half inervals on thiat model). These numbers
are still equally spaced.
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
| |
| Gordon Sande 2005-11-18, 7:56 am |
| On 2005-11-15 19:33:20 -0400, "Dr Ivan D. Reid" <Ivan.Reid@brunel.ac.uk> said:
> On Tue, 15 Nov 2005 14:11:23 GMT, robin <robin_v@bigpond.com>
> wrote in <f0mef.18282$Hj2.7323@news-server.bigpond.net.au>:
>
> We had this discussion a wee while ago, and some reasonable
> solutions were proposed. However, I'm still not convinced whether many
> applications require both endpoints to be included[1] (the one that
> started the discussion was a HEP problem that came up in a summer-school
> I attended in 2003). What's the application?
>
> [1] E.g. often a logarithm is involved so one endpoint cannot be exact!
In the jargon of measure theory this is a discussion about a set of
measure zero. In probability jargon it is about things which are
almost surely equal. Either way one gets rather curious about why
it matters. The convenience of not having to check for logarithm of
zero is quite real but that is a property of a PSEUDO random number.
If you are just trying to get past a "pointy haired boss" then
accept the psuedo random numbers on the range 0 <= x <= 1/2 and double.
Otherwise the suggestions, some polite and others a bit more direct, are
to understand what it is you are asking for and why it matters as it
will very often not make much sense.
| |
|
| James Giles wrote in message ...
>robin wrote:
>...
>
>Most random number generators return values in the
>interval 0 <= r < 1 that are uniformly spaced in that
>interval and are equiprobable.
Ideally, yes, but in practice, they are pseudo (note, pseudo)
and aren't perfect.
> Both of those last properties
>are very desirable. Your "solution"
It was a suggestion. Please note the wording:
"You could try multiplying your random number by
"NEAREST(1.0, 2.0)
"and, just for safety, checking that the product is
"not greater than 1.0. If it is, replace it with 1.0."
There are many options available to the OP, including
adapting and/or using one of the numerous PNRGs
available in computer science literature and on the internet.
> will not return values
>on the interval 0 <= r <= 1 that are uniformly spaced.
I am sure that everyone is aware of that. However,
virtually all are [uniformly spaced].
>For example, on IEEE single precision, you get 2^24
>values in the interval 0 <= r < 1 and the spacing between
>the possible results is 1/2^24.
Do you really think that anyone might notice the one in
some 16,000,000?
In a typical use, the value might never be generated anyway.
BTW, I have been writing pseudo-random number
generators and using PRNGs for 30+ years, so the
remainder is
<snipped>
>--
>J. Giles
| |
| Rich Townsend 2005-11-21, 7:01 pm |
| robin wrote:
> James Giles wrote in message ...
>
>
>
> Ideally, yes, but in practice, they are pseudo (note, pseudo)
> and aren't perfect.
>
>
>
>
> It was a suggestion. Please note the wording:
> "You could try multiplying your random number by
> "NEAREST(1.0, 2.0)
> "and, just for safety, checking that the product is
> "not greater than 1.0. If it is, replace it with 1.0."
>
> There are many options available to the OP, including
> adapting and/or using one of the numerous PNRGs
> available in computer science literature and on the internet.
>
>
>
>
> I am sure that everyone is aware of that. However,
> virtually all are [uniformly spaced].
>
>
>
>
> Do you really think that anyone might notice the one in
> some 16,000,000?
>
> In a typical use, the value might never be generated anyway.
>
Welcome to the school of faith-based programming!
Seriously, though, I've recently been bitten by assumptions that 1-in-10^6
occurences in a RNG-driven simulation are 'impossible'. I've recently been
working on a Monte-Carlo simulation of photon diffusion through a porous slab.
In tracking the photons from cell to cell, there are some situations where a
photon will propagate through the exact cell corner. However, such situations
are so rare that I decided they weren't going to actually arise in simulations,
and I didn't bother writing code to handle them.
Of couse, this came back to bite me on the arse. In some of the simulations, I
end up firing well over 10^6 photons through the slab -- enough such that the
probability of a cell-corner event in a single simulation becomes significant.
Hence, my code kept on falling over because I had designated a certain class of
events as 'improbable, hence impossible'.
cheers,
Rich
| |
| Dr Ivan D. Reid 2005-11-22, 7:01 pm |
| On Mon, 21 Nov 2005 10:45:36 -0500, Rich Townsend <rhdt@barVOIDtol.udel.edu>
wrote in <dlsq55$ajq$1@scrotar.nss.udel.edu>:
> robin wrote:
[color=darkred]
[color=darkred]
> Welcome to the school of faith-based programming!
> Seriously, though, I've recently been bitten by assumptions that 1-in-10^6
> occurences in a RNG-driven simulation are 'impossible'. I've recently been
> working on a Monte-Carlo simulation of photon diffusion through a porous slab.
> In tracking the photons from cell to cell, there are some situations where a
> photon will propagate through the exact cell corner. However, such situations
> are so rare that I decided they weren't going to actually arise in simulations,
> and I didn't bother writing code to handle them.
> Of couse, this came back to bite me on the arse. In some of the simulations, I
> end up firing well over 10^6 photons through the slab -- enough such that the
> probability of a cell-corner event in a single simulation becomes significant.
> Hence, my code kept on falling over because I had designated a certain class of
> events as 'improbable, hence impossible'.
It was something very similar that led to the discovery of ion
"channeling" in crystals, 30 odd years ago. Apparently it hadn't been
considered, until someone's simulation used a lot more CPU than expected
one night...
--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan.Reid@[brunel.ac.uk|cern.ch] Room 40-1-B12, CERN
KotPT -- "for stupidity above and beyond the call of duty".
|
|
|
|
|