
123

On Fri, Sep 01, 2017 at 04:00:59AM 0700, Michael Truog wrote:
> On 09/01/2017 01:54 AM, Raimo Niskanen wrote:
> > On Thu, Aug 31, 2017 at 10:29:34PM 0700, Michael Truog wrote:
> > :
> >> I have some examples that can make this desire a bit clearer:
> >>
> >> https://github.com/CloudI/cloudi_core/blob/a1c10a02245f0f4284d701a2ee5f07aad17f6e51/src/cloudi_core_i_runtime_testing.erl#L139L149> >>
> >> % use BoxMuller transformation to generate Gaussian noise
> >> % (G. E. P. Box and Mervin E. Muller,
> >> % A Note on the Generation of Random Normal Deviates,
> >> % The Annals of Mathematical Statistics (1958),
> >> % Vol. 29, No. 2 pp. 610–611)
> >> X1 = random(),
> >> X2 = PI2 * random(),
> >> K = StdDev * math:sqrt(2.0 * math:log(X1)),
> >> Result1 = erlang:max(erlang:round(Mean + K * math:cos(X2)), 1),
> >> Result2 = erlang:max(erlang:round(Mean + K * math:sin(X2)), 1),
> >> sleep(Result2),
> > Why not use rand:normal/3?
> >
> > It uses the Ziggurat Method and is supposed to be much faster and
> > numerically more stable than the basic BoxMuller method.
> >
> The BoxMuller is simpler and producing 2 results instead of 1 . I believe I looked at the source code for rand:normal/3 and expected the BoxMuller to be faster only because it creates 2 results, though I should check that. I will have to investigate it more.
Simpler  yes.
The basic benchmark in rand_SUITE indicates that rand:normal() is only
about 50% slower than rand:uniform(1 bsl 58) (internal word size),
which I think is a very good number.
The BoxMuller transform method needs 4 calls to the 'math' module for
nontrivial floating point functions i.e log(), sqrt(), cos() and sin(),
which is why I think that "must" be slower.
But I have also not measured... :/
Looking forward to hear your results!

/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On Fri, Sep 01, 2017 at 09:36:37AM +0000, Hugo Mills wrote:
> On Fri, Sep 01, 2017 at 10:41:15AM +0200, Raimo Niskanen wrote:
> > On Thu, Aug 31, 2017 at 10:29:34PM 0700, Michael Truog wrote:
> > > On 08/31/2017 07:57 PM, Richard A. O'Keefe wrote:
> > > >
> > > > On 1/09/17 6:35 AM, Michael Truog wrote:
> > > >> As I argued in the original pull request for these recent 20.0 random
> > > >> number changes, a uniform distribution is much more intuitive if it is
> > > >> inclusive: [0,1]
> > > >
> > > > Intuitive? For integers, I'll grant that. For reals? Not so much.
> > > > I certainly cannot grant "MUCH more intuitive".
> > > > I've had occasion to do (random next * 2  1) arcTanh, which of
> > > > course breaks down if *either* 0 or 1 is returned by (random next).
> > >
> > > A uniform distribution should be uniformly distributed. I understand the woes of floatingpoint prevent perfect uniform distribution, but we could at least try to pay attention to the limits involved, and if we did, that would make the idea much more intuitive.
> >
> > If I try to be philosophical, picking a random number in the range
> > 0.0 to 1.0 of real numbers, the probability of getting a number exactly 0.0
> > (or exactly 1.0) is infinitely low. Therefore the range (0.0,1.0) is more
> > natural.
>
> Mathematically, there's a distinction. What you've just described
> is that in a random variable over the interval [0.0, 1.0], 0.0 and 1.0
> happen *almost never* (which is a very specific technical term), and
> that values in the open interval (0.0, 1.0) occur *almost surely*.
>
> Being discrete, the computer implementation based on floating point
> numbers ensures that the probability of getting 0.0 or 1.0 in that
> case is measurably nonzero, whereas in the ideal version over the
> reals, above, it is infinitesimally small. In that distinction lie
> most of the problems that people are talking about here, I think.
Precisely!
>
> > > My belief is that the [0,1) distribution is the most common
> > > because it is the easiest to implement given the IEEE floating
> > > point standard format. However, I would also like to be proven
> > > wrong, to have more faith in the current situation.
>
> > I think that is very possible.
>
> From my relatively limited practical experience, either I've wanted
> [0, 1) or I don't care. Example:
>
> Red = int(random() * 256)
>
> where I don't want the value 256, because it's out of range for my
> 8bit graphics mode, but I do want the probability of 255 to be the
> same as every other value. So I want [0, 1) as my range.
>
> Alternatively:
>
> P = random(),
> if
> P =< 0.3 > ...;
> P =< 0.7 > ...;
> P > 0.7 > ...
> end
>
> where, in general, I don't care if I could get 0.0 or 1.0 or not,
> because the differences are immeasurably small for all practical
> purposes.
Especially since decimal numbers below 1.0 have no exact representation as
IEEE floating point numbers.
>
> I think it's clear to me that _several_ functions are needed for
> different cases, with fullyclosed, fullyopen and halfopen
> intervals. IMO, the fullyclosed and halfopen are probably the most
> useful (and, modulo any floatingpoint issues which I'm not qualified
> to talk about, [0,1) can be turned into (0,1] with
> 1random_halfopen()).
There should be no issues in this case with the updated algorithms in
the 'rand' module that produce numbers on the form N * 2^53, and due to
the fact that IEEE floating point arithmetics is defined to produce
correctly rounded results for the simple operations +, , *, /, ...
>
> With a fullyclosed interval, it should be possible to write
> helpers for generating the other three by simply calling
> random_closed() again if you get an undesirable endpoint. You can't
> easily extend the range of the halfopen or open intervals to give you
> the closed ones. So I'd say at minimum, there should be a function
> giving the closed interval.
The closed interval [0.0,1.0] may be the most general for the reasons you
mention. Unfortunately it is cumbersome to implement "efficiently".
You would have to generate integers [0, 1+2^53] and to do that with a 58
bit generator (the default) you generate one number and have roughly a
1/32 chance of landing in the high interval where you need to retry.
The amortized impact of this is only about 3% slower, but the code has to
do all the testing for these corner cases which adds up to maybe 20% slower.
But once in a million it will take 4 full attempts or more to get a number,
which is rather annoying, and the worst case is infinity (but will never
happen ;).
Therefore the half open interval [0.0,1.0) is probably the most useful one,
and the question is if the closed interval [0.0,1.0] (and the open interval
(0.0,1.0) is worth to implement...
>
> Whether the "test and retry" approach is the best implementation
> or not is a matter for discussion, as is the question of whether all
I have never heard of a better alternative than "test and retry" when you
are limiting the interval, except possibly when the probability for retry is
high (approaching 50%)  then you may consider to generate twice as many
bits and do a simple 'mod' operation for the result; the skew would be
impossible to notice. The "test and retry" method is still expected to be
as fast or faster, amortized, but the "generate double and 'mod'" method
has more or less fixed execution time.
> or some of these functions should be in the standard lib, or they are
> expected to be hacked together on an asneeded basis.
Yes. That is the question.
>
> Hugo.
>
> 
> Hugo Mills  Anyone who says their system is completely secure
> hugo@... carfax.org.uk  understands neither systems nor security.
> http://carfax.org.uk/ 
> PGP: E2AB1DE4  Bruce Schneier

/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On 09/01/2017 05:13 AM, Raimo Niskanen
wrote:
On Fri, Sep 01, 2017 at 04:00:59AM 0700, Michael Truog wrote:
On 09/01/2017 01:54 AM, Raimo Niskanen wrote:
On Thu, Aug 31, 2017 at 10:29:34PM 0700, Michael Truog wrote:
:
I have some examples that can make this desire a bit clearer:
https://github.com/CloudI/cloudi_core/blob/a1c10a02245f0f4284d701a2ee5f07aad17f6e51/src/cloudi_core_i_runtime_testing.erl#L139L149
% use BoxMuller transformation to generate Gaussian noise
% (G. E. P. Box and Mervin E. Muller,
% A Note on the Generation of Random Normal Deviates,
% The Annals of Mathematical Statistics (1958),
% Vol. 29, No. 2 pp. 610–611)
X1 = random(),
X2 = PI2 * random(),
K = StdDev * math:sqrt(2.0 * math:log(X1)),
Result1 = erlang:max(erlang:round(Mean + K * math:cos(X2)), 1),
Result2 = erlang:max(erlang:round(Mean + K * math:sin(X2)), 1),
sleep(Result2),
Why not use rand:normal/3?
It uses the Ziggurat Method and is supposed to be much faster and
numerically more stable than the basic BoxMuller method.
The BoxMuller is simpler and producing 2 results instead of 1 . I believe I looked at the source code for rand:normal/3 and expected the BoxMuller to be faster only because it creates 2 results, though I should check that. I will have to investigate it more.
Simpler  yes.
The basic benchmark in rand_SUITE indicates that rand:normal() is only
about 50% slower than rand:uniform(1 bsl 58) (internal word size),
which I think is a very good number.
The BoxMuller transform method needs 4 calls to the 'math' module for
nontrivial floating point functions i.e log(), sqrt(), cos() and sin(),
which is why I think that "must" be slower.
But I have also not measured... :/
Looking forward to hear your results!
I have some interesting results.
These results use https://github.com/okeuday/erlbench which includes
a copy of the source code at https://github.com/okeuday/quickrand :
TEST pseudo_randomness
N == 10000 (10 runs)
18_bxor_abs get: 1612.7 us ( 1.3)
18_erlang:system_tim get: 1254.1 us ( 1.0)
18_monotonic get: 1372.5 us ( 1.1)
18_os:system_time/1 get: 1221.7 us ( 1.0)
19_os:perf_counter/1 get: 3752.2 us ( 3.1)
20_rand:normal get: 6832.0 us ( 5.6)
20_rand_exrop get: 3949.3 us ( 3.2)
20_rand_exs1024s get: 12073.3 us ( 9.9)
20_rand_exsp get: 3390.4 us ( 2.8)
os:timestamp/0 get: 1392.3 us ( 1.1)
os_time:perf_counter get: 4109.4 us ( 3.4)
quickrand_c:floatR/0 get: 5776.0 us ( 4.7)
quickrand_c:floatR/1 get: 5704.3 us ( 4.7)
quickrand_c:uni/1 get: 4015.2 us ( 3.3)
quickrand_c:uni/2 get: 3960.7 us ( 3.2)
quickrand_c_normal/2 get: 9329.5 us ( 7.6)
quickrand_c_normal/3 get: 8917.7 us ( 7.3)
random_wh06_int:unif get: 10777.5 us ( 8.8)
random_wh82:uniform/ get: 4750.0 us ( 3.9)
random_wh82_int:unif get: 4866.4 us ( 4.0)
The function names that are relevant for a normal distribution
are:
20_rand:normal > rand:normal/0 (when
using rand:seed(exsp, _))
20_rand_exsp > rand:uniform/1 (when using rand:seed(exsp, _))
quickrand_c:floatR/0 > quickrand_cache:floatR/0
quickrand_c:floatR/1 > quickrand_cache:floatR/1
quickrand_c_normal/2 > quickrand_cache_normal:box_muller/2
quickrand_c_normal/3 >
quickrand_cache_normal:box_muller/3
The rand module exsp algorithm was used here because it is the
fastest pseudorandom number generator in the rand module.
A rough look at the latency associated with the normal
distribution method, ignoring the latency for random number source
is:
rand:normal/0
3441.6 us = 6832.0 us  (rand:uniform/1 3390.4
us)
quickrand_cache_normal:box_muller/2
3553.5 us = 9329.5 us  (quickrand_cache:floatR/0
5776.0 us)
quickrand_cache_normal:box_muller/3
3213.4 us = 8917.7
us  (quickrand_cache:floatR/1 5704.3
us)
So, this helps to
show that the latency with both methods is very similar if you
ignore the random number generation. However, it likely requires
some explanation: The quickrand_cache module is what I am using
here for random number generation, which stores cached data from crypto: strong_rand_bytes/1
with a default size of 64KB for the cache. The difference between
the functions quickrand_cache_normal:box_muller/2
and quickrand_cache_normal:box_muller/3
is that the first uses the process dictionary while the second
uses a state variable. Using the large amount of cached random
data, the latency associated with individual calls to crypto:strong_rand_bytes/1 is avoided at the cost of the
extra memory consumption, and the use of the cache makes the
speed of random number generation similar to the speed of
pseudorandom number generation that occurs in the rand module.
In
CloudI, I instead use quickrand_normal:box_muller/2 to avoid the
use of cached data to keep the memory use minimal (the usecase
there doesn't require avoiding the latency associated with crypto:strong_rand_bytes/1 because it is adding latency
for testing (at
https://github.com/CloudI/cloudi_core/blob/299df02e6d22103415c8ba14379e90ca8c3d3b82/src/cloudi_core_i_runtime_testing.erl#L138)
and it is best using a cryptographic random source to keep the
functionality widely applicable). However, the same function
calls occur in the quickrand BoxMuller transformation source
code, so the overhead is the same.
I used Erlang/OTP 20.0 (without HiPE) using the hardware below:
Core i7 2670QM 2.2GHz 1 cpu, 4 cores/cpu, 2 hts/core
L2:4×256KB L3:6MB RAM:8GB:DDR31333MHz
Sandy BridgeHE4 (Socket G2)
Best Regards,
Michael
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On 1/09/17 8:49 PM, Raimo Niskanen wrote:
>> By the way, given that a common way to make random floats is to
>> generate a bitvector, consider
>> (0 to: 15) collect: [:each  ((each / 15) * 256) truncated].
>> You will notice that the spacing between the values is *almost*
>> uniform, but not at the end.
>
> That sounds interesting but I do not understand. Is that Elixir code?
Nope, Smalltalk. I wanted to use rational arithmetic. In fact I did
not need to. Here it is in Haskell:
> [(x * 256) `div` 15  x < [0..15]]
[0,17,34,51,68,85,102,119,136,153,170,187,204,221,238,256]
Let's push that a bit further. Let's generate all possible 10bit
integers and map them to the range [0..63]. We find again that
the gap sizes are not all the same. They can't be. If you
consider all vectors of N bits and map them to the range
[0..2**M] they *cannot* be uniformly distributed no matter what
method you use because (2**M+1) does not divide 2**N. You can
fix this by rejecting some of the bit vectors, but that would
be asking everyone to pay extra for a result they don't have any
particular need for.
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On 2/09/17 12:53 AM, Raimo Niskanen wrote:
>> P = random(),
>> if
>> P =< 0.3 > ...;
>> P =< 0.7 > ...;
>> P > 0.7 > ...
>> end
> Especially since decimal numbers below 1.0 have no exact representation as
> IEEE floating point numbers.
There is in fact an IEEE standard for *decimal* floating point numbers.
If memory serves me correctly, z/Series and POWER support it. There is
a supplement to the C standard for it. Software emulation is available.
See decimal32, decimal64, decimal128 in
https://en.wikipedia.org/wiki/IEEE_754or the actual IEEE 7542008 if you have a copy (which I sadly don't).
In any case,
P = 10*rand(),
if P < 3 > ...
; P < 7 > ...
; P >= 7 > ...
end
evades that issue. (With the numbers from drand48() the multiplication
is exact; with 53bit random numbers it is not.)
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On Sat, Sep 02, 2017 at 04:35:39PM 0700, Michael Truog wrote:
> On 09/01/2017 05:13 AM, Raimo Niskanen wrote:
> > On Fri, Sep 01, 2017 at 04:00:59AM 0700, Michael Truog wrote:
> >> On 09/01/2017 01:54 AM, Raimo Niskanen wrote:
> >>> On Thu, Aug 31, 2017 at 10:29:34PM 0700, Michael Truog wrote:
> >>> :
> >>>> I have some examples that can make this desire a bit clearer:
> >>>>
> >>>> https://github.com/CloudI/cloudi_core/blob/a1c10a02245f0f4284d701a2ee5f07aad17f6e51/src/cloudi_core_i_runtime_testing.erl#L139L149> >>>>
> >>>> % use BoxMuller transformation to generate Gaussian noise
> >>>> % (G. E. P. Box and Mervin E. Muller,
> >>>> % A Note on the Generation of Random Normal Deviates,
> >>>> % The Annals of Mathematical Statistics (1958),
> >>>> % Vol. 29, No. 2 pp. 610–611)
> >>>> X1 = random(),
> >>>> X2 = PI2 * random(),
> >>>> K = StdDev * math:sqrt(2.0 * math:log(X1)),
> >>>> Result1 = erlang:max(erlang:round(Mean + K * math:cos(X2)), 1),
> >>>> Result2 = erlang:max(erlang:round(Mean + K * math:sin(X2)), 1),
> >>>> sleep(Result2),
> >>> Why not use rand:normal/3?
> >>>
> >>> It uses the Ziggurat Method and is supposed to be much faster and
> >>> numerically more stable than the basic BoxMuller method.
> >>>
> >> The BoxMuller is simpler and producing 2 results instead of 1 . I believe I looked at the source code for rand:normal/3 and expected the BoxMuller to be faster only because it creates 2 results, though I should check that. I will have to investigate it more.
> > Simpler  yes.
> >
> > The basic benchmark in rand_SUITE indicates that rand:normal() is only
> > about 50% slower than rand:uniform(1 bsl 58) (internal word size),
> > which I think is a very good number.
> >
> > The BoxMuller transform method needs 4 calls to the 'math' module for
> > nontrivial floating point functions i.e log(), sqrt(), cos() and sin(),
> > which is why I think that "must" be slower.
> >
> > But I have also not measured... :/
> >
> > Looking forward to hear your results!
> I have some interesting results.
>
> These results use https://github.com/okeuday/erlbench which includes a copy of the source code at https://github.com/okeuday/quickrand :
>
> TEST pseudo_randomness
> N == 10000 (10 runs)
> 18_bxor_abs get: 1612.7 us ( 1.3)
> 18_erlang:system_tim get: 1254.1 us ( 1.0)
> 18_monotonic get: 1372.5 us ( 1.1)
> 18_os:system_time/1 get: 1221.7 us ( 1.0)
> 19_os:perf_counter/1 get: 3752.2 us ( 3.1)
> 20_rand:normal get: 6832.0 us ( 5.6)
> 20_rand_exrop get: 3949.3 us ( 3.2)
> 20_rand_exs1024s get: 12073.3 us ( 9.9)
> 20_rand_exsp get: 3390.4 us ( 2.8)
> os:timestamp/0 get: 1392.3 us ( 1.1)
> os_time:perf_counter get: 4109.4 us ( 3.4)
> quickrand_c:floatR/0 get: 5776.0 us ( 4.7)
> quickrand_c:floatR/1 get: 5704.3 us ( 4.7)
> quickrand_c:uni/1 get: 4015.2 us ( 3.3)
> quickrand_c:uni/2 get: 3960.7 us ( 3.2)
> quickrand_c_normal/2 get: 9329.5 us ( 7.6)
> quickrand_c_normal/3 get: 8917.7 us ( 7.3)
> random_wh06_int:unif get: 10777.5 us ( 8.8)
> random_wh82:uniform/ get: 4750.0 us ( 3.9)
> random_wh82_int:unif get: 4866.4 us ( 4.0)
>
> The function names that are relevant for a normal distribution are:
> 20_rand:normal > rand:normal/0 (when using rand:seed(exsp, _))
> 20_rand_exsp > rand:uniform/1 (when using rand:seed(exsp, _))
> quickrand_c:floatR/0 > quickrand_cache:floatR/0
> quickrand_c:floatR/1 > quickrand_cache:floatR/1
> quickrand_c_normal/2 > quickrand_cache_normal:box_muller/2
> quickrand_c_normal/3 > quickrand_cache_normal:box_muller/3
>
> The rand module exsp algorithm was used here because it is the fastest pseudorandom number generator in the rand module.
Thank you for sharing these numbers!
>
> A rough look at the latency associated with the normal distribution method, ignoring the latency for random number source is:
> rand:normal/0
> 3441.6 us = 6832.0 us  (rand:uniform/1 3390.4 us)
Should not the base value come from rand:uniform/0 instead. I know the
difference is not big  rand_SUITE:measure/1 suggests 3%, but it also
suggests that rand:normal/0 is about 50% slower than rand:uniform/0 while
your numbers suggests 100% slower. Slightly strange...
> quickrand_cache_normal:box_muller/2
> 3553.5 us = 9329.5 us  (quickrand_cache:floatR/0 5776.0 us)
> quickrand_cache_normal:box_muller/3
> 3213.4 us = 8917.7 us  (quickrand_cache:floatR/1 5704.3us)
It is really interesting to see that the calls to the 'math' module
does not slow that algorithm down very much (hardly noticable)!
>
> So, this helps to show that the latency with both methods is very similar if you ignore the random number generation. However, it likely requires some explanation: The quickrand_cache module is what I am using here for random number generation, which stores cached data from crypto:strong_rand_bytes/1 with a default size of 64KB for the cache. The difference between the functions quickrand_cache_normal:box_muller/2 and quickrand_cache_normal:box_muller/3 is that the first uses the process dictionary while the second uses a state variable. Using the large amount of cached random data, the latency associated with individual calls to crypto:strong_rand_bytes/1 is avoided at the cost of the extra memory consumption, and the use of the cache makes the speed of random number generation similar to the speed of pseudorandom number generation that occurs in the rand module.
We should add a 'rand' plugin to the 'crypto' module that does this
buffered crypto:strong_random_bytes/1 trick. There is something like that
in rand_SUITE, but we should really have an official one.
I also wonder where the sweet spot is? 64 KB seems like a lot of buffer.
>
> In CloudI, I instead use quickrand_normal:box_muller/2 to avoid the use of cached data to keep the memory use minimal (the usecase there doesn't require avoiding the latency associated with crypto:strong_rand_bytes/1 because it is adding latency for testing (at https://github.com/CloudI/cloudi_core/blob/299df02e6d22103415c8ba14379e90ca8c3d3b82/src/cloudi_core_i_runtime_testing.erl#L138) and it is best using a cryptographic random source to keep the functionality widely applicable). However, the same function calls occur in the quickrand BoxMuller transformation source code, so the overhead is the same.
>
> I used Erlang/OTP 20.0 (without HiPE) using the hardware below:
> Core i7 2670QM 2.2GHz 1 cpu, 4 cores/cpu, 2 hts/core
> L2:4×256KB L3:6MB RAM:8GB:DDR31333MHz
> Sandy BridgeHE4 (Socket G2)
>
> Best Regards,
> Michael
> 
Best regards

/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


In reply to this post by Richard A. O'Keefe2
On Mon, Sep 04, 2017 at 12:37:50PM +1200, Richard A. O'Keefe wrote:
>
>
> On 1/09/17 8:49 PM, Raimo Niskanen wrote:
> >> By the way, given that a common way to make random floats is to
> >> generate a bitvector, consider
> >> (0 to: 15) collect: [:each  ((each / 15) * 256) truncated].
> >> You will notice that the spacing between the values is *almost*
> >> uniform, but not at the end.
> >
> > That sounds interesting but I do not understand. Is that Elixir code?
>
> Nope, Smalltalk. I wanted to use rational arithmetic. In fact I did
> not need to. Here it is in Haskell:
> > [(x * 256) `div` 15  x < [0..15]]
> [0,17,34,51,68,85,102,119,136,153,170,187,204,221,238,256]
>
> Let's push that a bit further. Let's generate all possible 10bit
> integers and map them to the range [0..63]. We find again that
> the gap sizes are not all the same. They can't be. If you
> consider all vectors of N bits and map them to the range
> [0..2**M] they *cannot* be uniformly distributed no matter what
> method you use because (2**M+1) does not divide 2**N. You can
> fix this by rejecting some of the bit vectors, but that would
> be asking everyone to pay extra for a result they don't have any
> particular need for.
>
I see, but do not quite understand what you are getting at.
The current leftclosed float generator starts with 58 random bits and
puts 53 of these into the mantissa in an IEEE 754 double binary float,
so that would not be it.
I guess it is a generator for the closed interval [0.0,1.0] or the open
(0.0,1.0) you talk about. If so:
This oneliner generates over [0.0,1.0]:
(rand:uniform((1 bsl 53) + 1) 1) * math:pow(2, 53)
and it uses an integer range R = (2^53 + 1), which is not dividable by 2.
The implementation for that range will generate a 58 bit number and then
check if the number is 0 =< X < R*31 and if so return the number mod R
(31 repetitions of the range is all that fits completely in 58 bits).
If the generated number is R*31 =< X that is in the top truncated interval
then we start over with a new number.
The implementation may in theory retry forever before finding a good
number, but the odds for retry is about 1/32 for each round so the
accumulated time is about 32/31 times one round. And only once in a million
it will take 4 attempts or more.
I discussed a different implementation with Prof. Vigna that is to always
generate one word too much and then use mod R on that which would bury the
skew in one word of random bits hence the difference in probability between
generated numbers should be about (2^58  1)/2^58, which would take quite
some effort to measure with statistical significance. But he considered
that as a bad idea since it would get half the speed for most ranges.
So this is an already solved problem, as I see it.
We *can* write efficient and good generators for open, closed and
halfclosed intervals, if we want.
So far I have only seen the need for a (0.0,1.0] generator, which can be
implemented with:
1.0  rand:uniform()
but in some applications such as 1.0/X and math:log(X) I think that the form
N * 2^53 might be less than optimal, so I have a new suggestion:
https://github.com/erlang/otp/compare/OTP20.0...RaimoNiskanen:raimo/stdlib/randuniformNZThis variant never returns exactly 0.0 and have better precision for low
values. Comments? Especially about the name.
And so far I have not seen any actual need for (0.0,1.0) nor [0.0,1.0].

/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


Raimo and all: I got late to follow the thread. I think the new function name should be rand:uniform_non_zero/{1,2} because I've rarely seen somethingUPPERCASE names in Erlang functions. (I might be wrong.) On ranges: I'm concerned about how OTP pre20 (i.e., <= OTP 19.3.x) rand:uniform/1 code might crash or cause bugs when running on the OTP 20 implementation. At least how to write the OTP pre20equivalent code should be documented. I have no firm idea on what should be the default behavior on ranges and whether the borders should be inclusive/exclusive to the limit values. In fact, the behavior differs between languages and implementations. Some examples I've found follow: JavaScript math.random(): [0.0, 1.0) https://developer.mozilla.org/enUS/docs/Web/JavaScript/Reference/Global_Objects/Math/randomPython random.uniform(a, b): [a, b] (when a <= b) https://stackoverflow.com/questions/6088077/howtogetarandomnumberbetweenafloatrangeC++ std::uniform_real_distribution<> gen(a, b): [a, b) http://en.cppreference.com/w/cpp/numeric/random/uniform_real_distributionRuby 2.4.1 Random class: rand(max): [0.0, max) https://rubydoc.org/core2.4.1/Random.html"When max is a Float, rand returns a random floating point number unless max = min or maxmin is small compared to min, and in particular not for the default arguments."
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On Wed, Sep 06, 2017 at 09:30:17PM +0900, Kenji Rikitake wrote:
> Raimo and all:
>
> I got late to follow the thread.
>
> I think the new function name should be
> rand:uniform_non_zero/{1,2}
> because I've rarely seen somethingUPPERCASE
> names in Erlang functions.
> (I might be wrong.)
At least it would be the first (only) in the rand module, so you surely
have a point.
Another possible name would be rand:uniform_right_closed/{1,2}, or
rand:uniform_left_open/{1,2} but the mathematical reference might be
a bit obscure.
>
> On ranges:
>
> I'm concerned about how OTP pre20 (i.e., <= OTP 19.3.x) rand:uniform/1 code
> might crash or cause bugs when running on the OTP 20 implementation.
> At least how to write the OTP pre20equivalent code should be documented.
If you mean pre20 code using rand:uniform/1 that later runs on 20.0 or
later, then I am not aware of any drastic changes, at least not crashes.
Code calling rand:seed/{1,2} will select the explicit algorithm from pre20
and behave exactly as before. Code not calling rand:seed/{1,2} will as
before get a random number sequence it has never seen before.
The same applies for code using rand:uniform/0  it could pre20 get
exactly 0.0 and the same on 20.0 and later.
You can, of course, only use the explicit algorithms that existed pre20
when executing pre20, and the same algorithms still exists but are not
documented on OTP 20.0 and should behave exactly the same.
So since there should be nothing more to writing pre20 code than to read
the documentation for the release in question, and that such code should
work just as before on 20.0 and later, I do not think there should be any
need for documenting that. This should be exactly as expected.
There is a note that the upper bound may depending on rounding not be
returned. Really nasty.
That is a really sloppy definition!
>
> MySQL 5.7 RAND(): [0.0, 1.0)
> https://dev.mysql.com/doc/refman/5.7/en/mathematicalfunctions.html#function_rand>
> PostgreSQL 10 random(): [0.0, 1.0)
> https://www.postgresql.org/docs/10/static/functionsmath.html#functionsmathrandomtable>
> MS SQL Server: (0.0, 1.0)
> https://docs.microsoft.com/enus/sql/tsql/functions/randtransactsql> "Returns a pseudorandom float value from 0 through 1, exclusive."
>
> dSFMT: "[1, 2), [0, 1), (0, 1] and (0, 1)"
> http://www.math.sci.hiroshimau.ac.jp/~mmat/MT/SFMT/#dSFMT> (dSFMT: SFMT PRNG implementation for doubleprecision floating point only)
>
> It took me an hour to investigate this.
> Lesson learned: don't take the range definitions for granted.
You serched better than me and did find some that excludes 0.0.
Good to know!
It seems [0.0, 1.0) dominates but does not rule.
Best regards
/ Raimo
>
> Regards,
> Kenji Rikitake
>
> On Mon, Sep 4, 2017 at 6:49 PM, Raimo Niskanen <
> [hidden email]> wrote:
>
> > On Mon, Sep 04, 2017 at 12:37:50PM +1200, Richard A. O'Keefe wrote:
> > >
> > >
> > > On 1/09/17 8:49 PM, Raimo Niskanen wrote:
> > > >> By the way, given that a common way to make random floats is to
> > > >> generate a bitvector, consider
> > > >> (0 to: 15) collect: [:each  ((each / 15) * 256) truncated].
> > > >> You will notice that the spacing between the values is *almost*
> > > >> uniform, but not at the end.
> > > >
> > > > That sounds interesting but I do not understand. Is that Elixir code?
> > >
> > > Nope, Smalltalk. I wanted to use rational arithmetic. In fact I did
> > > not need to. Here it is in Haskell:
> > > > [(x * 256) `div` 15  x < [0..15]]
> > > [0,17,34,51,68,85,102,119,136,153,170,187,204,221,238,256]
> > >
> > > Let's push that a bit further. Let's generate all possible 10bit
> > > integers and map them to the range [0..63]. We find again that
> > > the gap sizes are not all the same. They can't be. If you
> > > consider all vectors of N bits and map them to the range
> > > [0..2**M] they *cannot* be uniformly distributed no matter what
> > > method you use because (2**M+1) does not divide 2**N. You can
> > > fix this by rejecting some of the bit vectors, but that would
> > > be asking everyone to pay extra for a result they don't have any
> > > particular need for.
> > >
> >
> > I see, but do not quite understand what you are getting at.
> >
> > The current leftclosed float generator starts with 58 random bits and
> > puts 53 of these into the mantissa in an IEEE 754 double binary float,
> > so that would not be it.
> >
> > I guess it is a generator for the closed interval [0.0,1.0] or the open
> > (0.0,1.0) you talk about. If so:
> >
> > This oneliner generates over [0.0,1.0]:
> > (rand:uniform((1 bsl 53) + 1) 1) * math:pow(2, 53)
> > and it uses an integer range R = (2^53 + 1), which is not dividable by 2.
> >
> > The implementation for that range will generate a 58 bit number and then
> > check if the number is 0 =< X < R*31 and if so return the number mod R
> > (31 repetitions of the range is all that fits completely in 58 bits).
> >
> > If the generated number is R*31 =< X that is in the top truncated interval
> > then we start over with a new number.
> >
> > The implementation may in theory retry forever before finding a good
> > number, but the odds for retry is about 1/32 for each round so the
> > accumulated time is about 32/31 times one round. And only once in a
> > million
> > it will take 4 attempts or more.
> >
> > I discussed a different implementation with Prof. Vigna that is to always
> > generate one word too much and then use mod R on that which would bury the
> > skew in one word of random bits hence the difference in probability between
> > generated numbers should be about (2^58  1)/2^58, which would take quite
> > some effort to measure with statistical significance. But he considered
> > that as a bad idea since it would get half the speed for most ranges.
> >
> > So this is an already solved problem, as I see it.
> >
> > We *can* write efficient and good generators for open, closed and
> > halfclosed intervals, if we want.
> >
> > So far I have only seen the need for a (0.0,1.0] generator, which can be
> > implemented with:
> > 1.0  rand:uniform()
> > but in some applications such as 1.0/X and math:log(X) I think that the
> > form
> > N * 2^53 might be less than optimal, so I have a new suggestion:
> >
> > https://github.com/erlang/otp/compare/OTP20.0...
> > RaimoNiskanen:raimo/stdlib/randuniformNZ
> >
> > This variant never returns exactly 0.0 and have better precision for low
> > values. Comments? Especially about the name.
> >
> > And so far I have not seen any actual need for (0.0,1.0) nor [0.0,1.0].
> >
> > 
> >
> > / Raimo Niskanen, Erlang/OTP, Ericsson AB
> > _______________________________________________
> > erlangquestions mailing list
> > [hidden email]
> > http://erlang.org/mailman/listinfo/erlangquestions> >
> _______________________________________________
> erlangquestions mailing list
> [hidden email]
> http://erlang.org/mailman/listinfo/erlangquestions
/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On 7/09/17 12:30 AM, Kenji Rikitake wrote:
> Raimo and all:
>
> I got late to follow the thread.
>
> I think the new function name should be
> rand:uniform_non_zero/{1,2}
> because I've rarely seen somethingUPPERCASE
> names in Erlang functions.
> (I might be wrong.)
"nonzero" is one word, so rand:uniform_nonzero/{1,2}
would be better still.
>
> JavaScript math.random(): [0.0, 1.0)
> https://developer.mozilla.org/enUS/docs/Web/JavaScript/Reference/Global_Objects/Math/random>
> Python random.uniform(a, b): [a, b] (when a <= b)
> https://stackoverflow.com/questions/6088077/howtogetarandomnumberbetweenafloatrangeAd fontes! Check the official documentation.
https://docs.python.org/3/library/random.htmlsays quite explicitly
"Almost all module functions depend on the basic function random(),
which generates a random float uniformly in the semiopen range
[0.0, 1.0)."
The actual source code for random.uniform(a, b) is
def uniform(self, a, b):
"Get a random number in the range
[a, b) or [a, b] depending on rounding."
return a + (ba) * self.random()
Now of course in exact real arithmetic, if 0 <= u < 1
and a < b then a <= (ba)*u + a < b. So they are using
an algorithm that *would* exclude b except for roundoff.
And they are leaving it to users to deal with the
consequences of the numerical error, and weaselling out
of it by blaming the computer.
In any case, Python is a [0,1) example.
> C++ std::uniform_real_distribution<> gen(a, b): [a, b)
> http://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution>
> Ruby 2.4.1 Random class: rand(max): [0.0, max)
> https://rubydoc.org/core2.4.1/Random.html> "When max is a Float, rand returns a random floating point number
> between 0.0 and max, including 0.0 and excluding max."
>
> R runif(min=0.0, max=1.0): [0.0, 1.0] (See Note)
> https://stat.ethz.ch/Rmanual/Rdevel/library/stats/html/Uniform.html> Note: "runif will not generate either of the extreme values
> unless max = min or maxmin is small compared to min,
> and in particular not for the default arguments."
>
> MySQL 5.7 RAND(): [0.0, 1.0)
> https://dev.mysql.com/doc/refman/5.7/en/mathematicalfunctions.html#function_rand>
> PostgreSQL 10 random(): [0.0, 1.0)
> https://www.postgresql.org/docs/10/static/functionsmath.html#functionsmathrandomtable>
> MS SQL Server: (0.0, 1.0)
> https://docs.microsoft.com/enus/sql/tsql/functions/randtransactsql> "Returns a pseudorandom float value from 0 through 1, exclusive."
All of the systems you have mentioned to this point use [0,1)
as the building block (or in the case of R, (0,1).
>
> dSFMT: "[1, 2), [0, 1), (0, 1] and (0, 1)"
> http://www.math.sci.hiroshimau.ac.jp/~mmat/MT/SFMT/#dSFMT> (dSFMT: SFMT PRNG implementation for doubleprecision floating point only)
That is to say, it has these functions in the interface:
dsfmt_genrand_close1_open2() => [1,2)  their primitive
dsfmt_genrand_close_open() => [0,1)
dsfmt_genrand_open_close() => (0,1]
dsfmt_genrand_open_open() => (0,1)
You *can't* take the range for granted because you have to choose.
The one range it does not offer is [0,1].
The key lesson is that none of these systems offers [0,1] as a
basic building block, and that to the extent that it is possible
to tell, the ones that offer [a,b] as a derived version  R and
Python  do so as a sloppy implementation of [a,b).
For what it's worth, my Smalltalk library now has
aRandom
next [0,1)
nextNonzero (0,1]
nextBetween: a andExclusive: b [a,b)
nextBetweenExclusive: a and: b (a,b]
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


To conclude
===========
It would be convenient to have functions in the rand module that
generates on the interval (0.0, 1.0] e.g rand:uniform_nonzero/0
and rand:uniform_nonzero_s/1 or maybe rand:uniform_nz/0
and rand:uniform_nz_s/1.
But since it is very easy to use e.g (1.0  rand:uniform()) I see
little value in adding them. Adding hints in the documentation could
suffice.
However, I think we could add functions with the same names and
interval that also has got a different uniformity i.e still uniform,
but not equidistant, instead increasing precison towards 0. This would
have a greater value. Such a uniformity would work better for some
suggested algorithms such as BoxMuller.
Ironically, the implementation by Kenji that I changed was a few bits in
that direction i.e used the generators' extra bits over 53 (5 or 11)
for increased precision, but I want to have at least 53 extra bits which I
hope is close enough to infinity.
I have implemented such functions in my own GitHub repo, but ran into
problems with the number distribution that I suspect is caused by
rounding errors in the current implementation of erlang:float/1.
So I will try to find the time to investigate that further...
/ Raimo
On Thu, Sep 07, 2017 at 12:02:36PM +1200, Richard A. O'Keefe wrote:
>
>
> On 7/09/17 12:30 AM, Kenji Rikitake wrote:
> > Raimo and all:
> >
> > I got late to follow the thread.
> >
> > I think the new function name should be
> > rand:uniform_non_zero/{1,2}
> > because I've rarely seen somethingUPPERCASE
> > names in Erlang functions.
> > (I might be wrong.)
>
> "nonzero" is one word, so rand:uniform_nonzero/{1,2}
> would be better still.
> >
> > JavaScript math.random(): [0.0, 1.0)
> > https://developer.mozilla.org/enUS/docs/Web/JavaScript/Reference/Global_Objects/Math/random> >
> > Python random.uniform(a, b): [a, b] (when a <= b)
> > https://stackoverflow.com/questions/6088077/howtogetarandomnumberbetweenafloatrange>
> Ad fontes! Check the official documentation.
> https://docs.python.org/3/library/random.html> says quite explicitly
> "Almost all module functions depend on the basic function random(),
> which generates a random float uniformly in the semiopen range
> [0.0, 1.0)."
>
> The actual source code for random.uniform(a, b) is
>
> def uniform(self, a, b):
> "Get a random number in the range
> [a, b) or [a, b] depending on rounding."
> return a + (ba) * self.random()
>
> Now of course in exact real arithmetic, if 0 <= u < 1
> and a < b then a <= (ba)*u + a < b. So they are using
> an algorithm that *would* exclude b except for roundoff.
> And they are leaving it to users to deal with the
> consequences of the numerical error, and weaselling out
> of it by blaming the computer.
>
> In any case, Python is a [0,1) example.
>
> > C++ std::uniform_real_distribution<> gen(a, b): [a, b)
> > http://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution> >
> > Ruby 2.4.1 Random class: rand(max): [0.0, max)
> > https://rubydoc.org/core2.4.1/Random.html> > "When max is a Float, rand returns a random floating point number
> > between 0.0 and max, including 0.0 and excluding max."
> >
> > R runif(min=0.0, max=1.0): [0.0, 1.0] (See Note)
> > https://stat.ethz.ch/Rmanual/Rdevel/library/stats/html/Uniform.html> > Note: "runif will not generate either of the extreme values
> > unless max = min or maxmin is small compared to min,
> > and in particular not for the default arguments."
> >
> > MySQL 5.7 RAND(): [0.0, 1.0)
> > https://dev.mysql.com/doc/refman/5.7/en/mathematicalfunctions.html#function_rand> >
> > PostgreSQL 10 random(): [0.0, 1.0)
> > https://www.postgresql.org/docs/10/static/functionsmath.html#functionsmathrandomtable> >
> > MS SQL Server: (0.0, 1.0)
> > https://docs.microsoft.com/enus/sql/tsql/functions/randtransactsql> > "Returns a pseudorandom float value from 0 through 1, exclusive."
>
> All of the systems you have mentioned to this point use [0,1)
> as the building block (or in the case of R, (0,1).
> >
> > dSFMT: "[1, 2), [0, 1), (0, 1] and (0, 1)"
> > http://www.math.sci.hiroshimau.ac.jp/~mmat/MT/SFMT/#dSFMT> > (dSFMT: SFMT PRNG implementation for doubleprecision floating point only)
>
> That is to say, it has these functions in the interface:
> dsfmt_genrand_close1_open2() => [1,2)  their primitive
> dsfmt_genrand_close_open() => [0,1)
> dsfmt_genrand_open_close() => (0,1]
> dsfmt_genrand_open_open() => (0,1)
>
> You *can't* take the range for granted because you have to choose.
> The one range it does not offer is [0,1].
>
> The key lesson is that none of these systems offers [0,1] as a
> basic building block, and that to the extent that it is possible
> to tell, the ones that offer [a,b] as a derived version  R and
> Python  do so as a sloppy implementation of [a,b).
>
> For what it's worth, my Smalltalk library now has
> aRandom
> next [0,1)
> nextNonzero (0,1]
> nextBetween: a andExclusive: b [a,b)
> nextBetweenExclusive: a and: b (a,b]
> _______________________________________________
> erlangquestions mailing list
> [hidden email]
> http://erlang.org/mailman/listinfo/erlangquestions
/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On Mon, Sep 04, 2017 at 10:48:15AM +0200, Raimo Niskanen wrote:
:
> We should add a 'rand' plugin to the 'crypto' module that does this
> buffered crypto:strong_random_bytes/1 trick. There is something like that
> in rand_SUITE, but we should really have an official one.
>
> I also wonder where the sweet spot is? 64 KB seems like a lot of buffer.
>
I have created a pull request for such an extension to the crypto module:
https://github.com/erlang/otp/pull/1573

/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions


On Fri, Sep 08, 2017 at 02:26:28PM +0200, Raimo Niskanen wrote:
> To conclude
> ===========
>
> It would be convenient to have functions in the rand module that
> generates on the interval (0.0, 1.0] e.g rand:uniform_nonzero/0
> and rand:uniform_nonzero_s/1 or maybe rand:uniform_nz/0
> and rand:uniform_nz_s/1.
>
> But since it is very easy to use e.g (1.0  rand:uniform()) I see
> little value in adding them. Adding hints in the documentation could
> suffice.
>
> However, I think we could add functions with the same names and
> interval that also has got a different uniformity i.e still uniform,
> but not equidistant, instead increasing precison towards 0. This would
> have a greater value. Such a uniformity would work better for some
> suggested algorithms such as BoxMuller.
>
> Ironically, the implementation by Kenji that I changed was a few bits in
> that direction i.e used the generators' extra bits over 53 (5 or 11)
> for increased precision, but I want to have at least 53 extra bits which I
> hope is close enough to infinity.
>
> I have implemented such functions in my own GitHub repo, but ran into
> problems with the number distribution that I suspect is caused by
> rounding errors in the current implementation of erlang:float/1.
>
> So I will try to find the time to investigate that further...
I have created a pull request with such an extension to the rand module:
https://github.com/erlang/otp/pull/1574
/ Raimo Niskanen, Erlang/OTP, Ericsson AB
_______________________________________________
erlangquestions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlangquestions

123
