Want to generate random number in between large specific range

[Pratik Patel]

Hello,

I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:

call random_number(u)
n+(m+1-n)*u (where, n = 0 , m =10^-6)

But it didn't work.
can anyone help me with this issue?

Thanks.

Pratik

[Thomas Koenig]

Pratik Patel <pratikm...@gmail.com> schrieb:

If you have a problem like that, it is always good to share
the whole code, which can answer questions like the following:
Did you declare your variables? Undelcared variables in Fortran
are integer if their first letter is in the range of i to n.
It is also a convention to use variables that way; "n" is almost
always a number of entities.

As a general rule, you should always use IMPLICIT NONE to detect
any you left out.

Also, your formula seems off. If u=0, you get 0, and if u=1, you
get n+m+1-n = 0+1e-6+1 = 1.000001.

This works for me (with variables renamed to what I would use):

program memain
implicit none
real :: lower, upper
real :: a
lower = 0.
upper = 1e-6
call random_number(a)
print *,lower+a*(upper-lower)
end program memain

[gah4]

On Wednesday, November 23, 2022 at 9:19:04 PM UTC-8, Pratik Patel wrote:

> I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:

> call random_number(u)
> n+(m+1-n)*u (where, n = 0 , m =10^-6)

I believe this formula is right for m and n integers, and an integer range.
And more specifically, the result is truncated to an integer value.

In addition, note that for a continuous probability distribution,
the probability of any given (real) value is zero.

There is, then, no difference between that and [0 10^-6), which you would
get from m*u (assuming n=0). And especially, be sure than m and n are not
integer variables.

[Phillip Helbig (undress to reply)]

In article <c6a073b5-17df-45a1...@googlegroups.com>,

Pratik Patel <pratikm...@gmail.com> writes:

> I want to generate a real random number from [0 10^-6].

Use a standard program to generate them in [0,1] then divide by 10**6.

[Robin Vowels]

On Thursday, November 24, 2022 at 8:06:37 PM UTC+11, gah4 wrote:
> On Wednesday, November 23, 2022 at 9:19:04 PM UTC-8, Pratik Patel wrote:
>
> > I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:
>
> > call random_number(u)
> > n+(m+1-n)*u (where, n = 0 , m =10^-6)
> I believe this formula is right for m and n integers,

.
m = 10^-6 yields zero.
.

[pehache]

call random_number(u)
u = u * 1.0e-6

Am I missing something ??

--
"...sois ouvert aux idées des autres pour peu qu'elles aillent dans le
même sens que les tiennes.", ST sur fr.bio.medecine
ST passe le mur du çon : <j3nn2h...@mid.individual.net>

[gah4]

On Tuesday, November 29, 2022 at 11:49:41 PM UTC-8, pehache wrote:

(snip)

> > I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:

(snip)

> call random_number(u)
> u = u * 1.0e-6

> Am I missing something ??

The OP may, or may not, have wanted a non-zero lower value.

But also, the OP wanted an inclusive interval.

Mostly I believe that the one shouldn't worry about that.
In the case of a continuous distribution, the probability of any
individual value is zero. (Unless there is a delta function
in the distribution.)

In a digital computer, distributions aren't continuous, but they are close.
Close enough.

First note that for a range that isn't a half open interval that is a power
of two, the actual probability of a given representable value isn't uniform.

There are about 2**52 IEEE doubles between 0 and 1, so you might have
to see 2**52 values before seeing any specific value.

You would need many more that that, to know that a specific value
never occurred.

So, yes, u*1e-6.

[Robin Vowels]

On Wednesday, November 30, 2022 at 8:03:50 PM UTC+11, gah4 wrote:
> On Tuesday, November 29, 2022 at 11:49:41 PM UTC-8, pehache wrote:
>
> (snip)
> > > I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:
> (snip)
> > call random_number(u)
> > u = u * 1.0e-6
>
> > Am I missing something ??

.

> The OP may, or may not, have wanted a non-zero lower value.
>
> But also, the OP wanted an inclusive interval.
>
> Mostly I believe that the one shouldn't worry about that.
> In the case of a continuous distribution, the probability of any
> individual value is zero. (Unless there is a delta function
> in the distribution.)
>
> In a digital computer, distributions aren't continuous, but they are close.
> Close enough.
>
> First note that for a range that isn't a half open interval that is a power
> of two, the actual probability of a given representable value isn't uniform.
>
> There are about 2**52 IEEE doubles between 0 and 1, so you might have
> to see 2**52 values before seeing any specific value.
>
> You would need many more that that, to know that a specific value
> never occurred.
>
> So, yes, u*1e-6.

.
try u*1d-6

[gah4]

On Wednesday, November 30, 2022 at 7:19:19 PM UTC-8, Robin Vowels wrote:

(snip, someone wrote)

> > > > I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:

(and I wrote)

> > So, yes, u*1e-6.

> try u*1d-6

Maybe, or maybe not.

This whole problem depends on rounding more than usual.

First, neither 1e-6 or 1d-6 has an exact binary representation.
(But maybe the OP has a machine with DFP?)

First, we don't know the type of u.
Second, we don't know the type needed for the final value.
And finally, we don't know the exact range of values needed.

There will be rounding after the multiply, and possibly again
on assignment to a result variable.

Either or both of 1e-6 and 1d-6 might be more than the actual
desired value.

It might be important that the result possibly be larger
than either the true or represented value of the constant,
or that it never be larger.

[Robin Vowels]

On Thursday, December 1, 2022 at 7:44:09 PM UTC+11, gah4 wrote:
> On Wednesday, November 30, 2022 at 7:19:19 PM UTC-8, Robin Vowels wrote:
>
> (snip, someone wrote)
> > > > > I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:
> (and I wrote)
>
> > > So, yes, u*1e-6.
>
> > try u*1d-6
>
> Maybe, or maybe not.
>
> This whole problem depends on rounding more than usual.
>
> First, neither 1e-6 or 1d-6 has an exact binary representation.

.
If that bothers you, use u/1d6

>
> First, we don't know the type of u.

.
unimportant.

.
> Second, we don't know the type needed for the final value.

.
immaterial.

.
> And finally, we don't know the exact range of values needed.

.
you do. The OP told you.

.
> There will be rounding after the multiply, and possibly again
> on assignment to a result variable.
>
> Either or both of 1e-6 and 1d-6 might be more than the actual
> desired value.

.
see above.

[gah4]

On Thursday, December 1, 2022 at 5:52:52 PM UTC-8, Robin Vowels wrote:

(snip, I wrote)

> > Either or both of 1e-6 and 1d-6 might be more than the actual
> > desired value.

> see above.

So you know which one is more, or less, than the actual decimal value?

Would you like to tell the rest of us?

[gah4]

write(*,'(2E30.20)') 1e-6, 1d-6
end

[Ron Shepard]

On 12/1/22 2:44 AM, gah4 wrote:
> On Wednesday, November 30, 2022 at 7:19:19 PM UTC-8, Robin Vowels wrote:
>
> (snip, someone wrote)
>
>>>>> I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:
>
> (and I wrote)
>
>>> So, yes, u*1e-6.
>
>> try u*1d-6

These days, I would never recommend to anyone to use a D exponent. One
should define the appropriate KIND parameters, and use those.

u*1.0e-6_wp

> Maybe, or maybe not.
>
> This whole problem depends on rounding more than usual.
>
> First, neither 1e-6 or 1d-6 has an exact binary representation.

> (But maybe the OP has a machine with DFP?)

When a programmer uses 1.0e-6_wp (or the other versions), he is
specifying an exact floating point number that has a precise binary
representation. The problem in this case is that it is up to the
compiler exactly which number that is because, as you say, it could
round up to slightly larger or down to slightly smaller than the
mathematical number.

> First, we don't know the type of u.

> Second, we don't know the type needed for the final value.

> And finally, we don't know the exact range of values needed.

The OP told us that the range was [0,10^-6], but that may not help much
because, first, that is equivalent to [0,0] because of the integer
expressions, and second, as you point out, the floating point version of
that expression is lacking in detail about the real kind. So, lets make
a few assumptions and explore this issue a little more.

First, lets assume that the OP meant to use the range [0.0_wp, 1.e-6_wp]
for some yet to be specified working precision. That lower bound has no
ambiguity because the mathematical expression maps exactly to the
floating point value. The upper bound has some mapping to an exact
floating point value, so lets assume that that value is the one intended
to be used, and that he wants to return that upper bound value.

The PRNG will return a value u in the range [0.0_wp,1.0_wp). I think
that upper bound is reliable, meaning that u<1.0_wp will always be
satisfied. That means the product relation u*1.0e-6_wp<1.0e-6_wp will
always be satisfied if all the bits are computed correctly in the fp
multiplication. That means that the multiplication would never return
the desired upper bound. Is that acceptable? Who knows? Fortran does not
guarantee correct fp multiplication, of course, but if the upper bound
is not satisfied, then every multiplication must be tested against the
upper bound anyway. If that test is really necessary, then an expression
like

min( u*1.0e-6_wp, 1.0e-6_wp )

could be used. To my eye, that test looks like it would be overkill, but
it is possible that it is really critical to the program that the bound
be satisfied, so that kind of expression would work.

But what if the OP was a little imprecise, and he really wanted the
range [0.0_wp,1.e-6_wp)? That is, he does not want the upper bound value
to ever be generated. In this case, one could remember that we really
aren't working with real numbers, we are working with just a finite
subset of rational numbers. Standard fortran gives us an easy way to
specify what is the correct upper bound to the desired range, and it is

ub = nearest( 1.0e-6_wp -1.0_wp )

Now, the question is what values to use in the clipped expression if we
want to rigorously guard against returning the upper bound. I'm not sure
what is the answer. Two possibilities are

min( u*1.0e-6_wp, ub )
min( u*ub, ub )

I think that first expression is probably correct. I think the second
expression might always miss the desired upper bound by two spacing
units with exactly rounded multiplication, not just one spacing unit.

$.02 -Ron Shepard

[Ron Shepard]

On 11/30/22 3:03 AM, gah4 wrote:
[...]

> There are about 2**52 IEEE doubles between 0 and 1, so you might have
> to see 2**52 values before seeing any specific value.
>
> You would need many more that that, to know that a specific value
> never occurred.

There are 2**52-1 IEEE doubles between .5 and 1.

Most PRNGs will produce numbers that are equally distributed in the
interval [0,1). They do this by taking random integers with the
appropriate number of bits and scaling them by the floating point output
range. So this means that some subset of the floating point numbers in
that range will be selected, presumably with equal probability, and the
remaining subset will never be selected, no matter how long you wait for
them to occur.

Take for example, the values in the range [.5,1.0). All of those values
have the same exponent and therefore the same spacing between them. Then
consider the range [.25,.5). Those are the same mantissa bits with a
different exponent, so the size of the two subsets are the same, but the
spacing is half the size as that of the larger values. But if the PRNG
selects them by simply dividing, then only 1/2 of the [.25,0.5) set will
be chosen compared to the [.5,1.0) set. The subinterval [.125,.25) will
only contribute 1/4 of its values to the PRNG, the subinterval
[.0625,.125) will only contribute 1/8 of its values, and so on. Once the
exponent is small enough, i.e. below the minimum PRNG value produced,
then no values at all will be contributed to the PRNG set from those
exponent ranges.

I have always thought it would be interesting if instead PRNGs would
cycle over all possible floating point values. All the fp values in each
exponent range would eventually be selected, but with probabilities
associated with the exponent value. That is elements from the [.5,1.0)
range would occur with twice the frequency as from the [.25,.5) range,
and so on. I do not keep up with the PRNG literature, so maybe this
approach has been implemented somewhere already. Anyone know?

$.02 -Ron Shepard

[gah4]

On Friday, December 2, 2022 at 8:24:33 AM UTC-8, Ron Shepard wrote:

(snip)

> First, lets assume that the OP meant to use the range [0.0_wp, 1.e-6_wp]
> for some yet to be specified working precision. That lower bound has no
> ambiguity because the mathematical expression maps exactly to the
> floating point value. The upper bound has some mapping to an exact
> floating point value, so lets assume that that value is the one intended
> to be used, and that he wants to return that upper bound value.

> The PRNG will return a value u in the range [0.0_wp,1.0_wp). I think
> that upper bound is reliable, meaning that u<1.0_wp will always be
> satisfied.

> That means the product relation u*1.0e-6_wp<1.0e-6_wp will
> always be satisfied if all the bits are computed correctly in the fp
> multiplication.

I don't know that is so obvious. If the multiply rounds, it might be
able to round up. Before IEEE 754, each system decided how it
would do rounding. IEEE 754 allows one to choose the rounding,
but even so, I am not sure about this case.

> That means that the multiplication would never return
> the desired upper bound. Is that acceptable? Who knows? Fortran does not
> guarantee correct fp multiplication, of course, but if the upper bound
> is not satisfied, then every multiplication must be tested against the
> upper bound anyway. If that test is really necessary, then an expression
> like

(snip)

> But what if the OP was a little imprecise, and he really wanted the
> range [0.0_wp,1.e-6_wp)? That is, he does not want the upper bound value
> to ever be generated. In this case, one could remember that we really
> aren't working with real numbers, we are working with just a finite
> subset of rational numbers. Standard fortran gives us an easy way to
> specify what is the correct upper bound to the desired range, and it is

Since the OP was using a formula designed to give closed interval
for integers, I suspect the goal was to find the closed interval for
floating point values.

(snip)

[gah4]

On Friday, December 2, 2022 at 8:35:28 AM UTC-8, Ron Shepard wrote:

(snip)

> I have always thought it would be interesting if instead PRNGs would
> cycle over all possible floating point values. All the fp values in each
> exponent range would eventually be selected, but with probabilities
> associated with the exponent value. That is elements from the [.5,1.0)
> range would occur with twice the frequency as from the [.25,.5) range,
> and so on. I do not keep up with the PRNG literature, so maybe this
> approach has been implemented somewhere already. Anyone know?

My first thought was that you were asking about a log distribution,
until I figured out that was not what it does.

A fairly common use for random numbers want Gaussian distribution.
There is a formula to convert a uniform distribution to a Gaussian.

I haven't thought about that for a while, but one complication is that
it doesn't always get the tails right. Gaussians should have a small
probability of getting very far out. Given the usual [0, 1) values,
they don't get all that far, but it might be that they could do better.

[Ron Shepard]

On 12/2/22 2:43 PM, gah4 wrote:
> I don't know that is so obvious. If the multiply rounds, it might be
> able to round up. Before IEEE 754, each system decided how it
> would do rounding. IEEE 754 allows one to choose the rounding,
> but even so, I am not sure about this case.

When I posted that, I was convinced it was correct, but now it does seem
like rounding could violate that bound. At the time I posted, I was
thinking about the case of inaccurate multiplications and divisions,
like the CRAY hardware, and thus the min() expressions using nearest().

$.02 -Ron Shepard

[Phillip Helbig (undress to reply)]

In article <fPpiL.8424$7E38...@fx01.iad>, Ron Shepard

<nos...@nowhere.org> writes:

> Most PRNGs will produce numbers that are equally distributed in the
> interval [0,1). They do this by taking random integers with the
> appropriate number of bits and scaling them by the floating point output
> range. So this means that some subset of the floating point numbers in
> that range will be selected, presumably with equal probability, and the
> remaining subset will never be selected, no matter how long you wait for
> them to occur.

True of many, probably most, RNGs. However, RANLUX not only produces
all possible values, but the period is MUCH longer than the number of
possible values.

> I have always thought it would be interesting if instead PRNGs would
> cycle over all possible floating point values.

RANLUX does.

To get back to the OP: If you want random numbers between 0 and x where
x<1, why not just generate random deviates (equally spaced---on
average---between 0 and 1), which is what most RNGs do, then divide by
x?

[Phillip Helbig (undress to reply)]

In article <f2dde822-6a92-40ba...@googlegroups.com>,

gah4 <ga...@u.washington.edu> writes:

> A fairly common use for random numbers want Gaussian distribution.
> There is a formula to convert a uniform distribution to a Gaussian.
>
> I haven't thought about that for a while, but one complication is that
> it doesn't always get the tails right. Gaussians should have a small
> probability of getting very far out. Given the usual [0, 1) values,
> they don't get all that far, but it might be that they could do better.

There are formulae for the conversion.

For Gaussian distributions, generate N uniform-deviate random numbers,
add them and divide by N. That gives you your first Gaussian random
number. Repeat as often as you need to get more Gaussian random
numbers. Efficient? No.

[gah4]

On Saturday, December 3, 2022 at 9:30:55 AM UTC-8, Phillip Helbig (undress to reply) wrote:
> In article <f2dde822-6a92-40ba...@googlegroups.com>,

(snip about Gaussians)

> There are formulae for the conversion.

> For Gaussian distributions, generate N uniform-deviate random numbers,
> add them and divide by N. That gives you your first Gaussian random
> number. Repeat as often as you need to get more Gaussian random
> numbers. Efficient? No.

https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

It depends in a certain way about the statistics of the smallest input values.

I suspect with your distribution, it would give better results.

[Ron Shepard]

On 12/3/22 11:30 AM, Phillip Helbig (undress to reply) wrote:
> For Gaussian distributions, generate N uniform-deviate random numbers,
> add them and divide by N. That gives you your first Gaussian random
> number.

This would actually be a binomial distribution. While it is true that
binomial distributions approach a gaussian shape for large N, they have
some important differences, e.g. in the tail regions, as mentioned in
some previous posts.

$.02 -Ron Shepard

[Phillip Helbig (undress to reply)]

In article <oy4jL.61538$pem1....@fx10.iad>, Ron Shepard

Yes, in the limit of large N.

It is also the explanation why measurement errors are Gaussian: they are
the sum of random errors.

[David Jones]

Not a binomial distribution. For example, starting from an exact
uniform distribtion at N=1, for N=2 you get a triangular distribution.
You could find the eaxct formulae for larger N in various books, but
there would be little point in that. To get approximate standard
Gaussian, you would take the sum, subract N/2 divide the result by the
square root of the variance, which is N/12 ... which leads to using the
value N=12 as a common choice when this approach was in use.

[gah4]

On Monday, December 5, 2022 at 7:52:37 AM UTC-8, David Jones wrote:

(snip)

> Not a binomial distribution. For example, starting from an exact
> uniform distribtion at N=1, for N=2 you get a triangular distribution.
> You could find the eaxct formulae for larger N in various books, but
> there would be little point in that. To get approximate standard
> Gaussian, you would take the sum, subract N/2 divide the result by the
> square root of the variance, which is N/12 ... which leads to using the
> value N=12 as a common choice when this approach was in use.

I remember this one, from the subroutine GAUSS in the IBM SSP
(Scientific Subroutine Package), from before I knew about
Gaussian distributions.

I was remembering the 6 (number subtracted) instead of the 12,
but still remembering, at the time, not knowing what was
special about 12, other than that it was special.

Note that SSP is also where RANDU comes from, everyone's
least favorite random number generator.

I used it for some high school projects, when I didn't have
any other one to use.

[gah4]

On Wednesday, November 23, 2022 at 9:19:04 PM UTC-8, Pratik Patel wrote:

> I want to generate a real random number from [0 10^-6]. I have been unable to do so using the following logic:

I just noticed that the subject says "large specific range".

Is [0 10^-6] really supposed to be a large range?