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#L139-L149 > >> > >> % use Box-Muller 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 Box-Muller method. > > > The Box-Muller is simpler and producing 2 results instead of 1 . I believe I looked at the source code for rand:normal/3 and expected the Box-Muller 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 Box-Muller transform method needs 4 calls to the 'math' module for non-trivial 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 _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Hugo Mills-2
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 floating-point 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 non-zero, 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 > 8-bit 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 fully-closed, fully-open and half-open > intervals. IMO, the fully-closed and half-open are probably the most > useful (and, modulo any floating-point issues which I'm not qualified > to talk about, [0,1) can be turned into (0,1] with > 1-random_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 fully-closed interval, it should be possible to write > helpers for generating the other three by simply calling > random_closed() again if you get an undesirable end-point. You can't > easily extend the range of the half-open 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 as-needed 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 _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Raimo Niskanen-2
On Thu, Aug 31, 2017 at 02:30:08PM +0200, Raimo Niskanen wrote:
: > Should I make a pull request of this? > > https://github.com/erlang/otp/compare/OTP-20.0...RaimoNiskanen:raimo/stdlib/rand-uniformNZ > > Is the name uniformNZ good enough? > Are uniform floats complete enough with this addition? I have updated the suggestion above. Comments? -- / Raimo Niskanen, Erlang/OTP, Ericsson AB _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Raimo Niskanen-2
On 09/01/2017 05:13 AM, Raimo Niskanen
wrote:
I have some interesting results.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#L139-L149 % use Box-Muller 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 Box-Muller method.The Box-Muller is simpler and producing 2 results instead of 1 . I believe I looked at the source code for rand:normal/3 and expected the Box-Muller 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 Box-Muller transform method needs 4 calls to the 'math' module for non-trivial 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! 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 pseudo-random 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 pseudo-random 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 use-case 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 Box-Muller 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
_______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Raimo Niskanen-2
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 10-bit 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. _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Raimo Niskanen-2
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_754 or the actual IEEE 754-2008 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 53-bit random numbers it is not.) _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Michael Truog
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#L139-L149 > >>>> > >>>> % use Box-Muller 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 Box-Muller method. > >>> > >> The Box-Muller is simpler and producing 2 results instead of 1 . I believe I looked at the source code for rand:normal/3 and expected the Box-Muller 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 Box-Muller transform method needs 4 calls to the 'math' module for > > non-trivial 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 pseudo-random 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 pseudo-random 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 use-case 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 Box-Muller 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:DDR3-1333MHz > Sandy Bridge-HE-4 (Socket G2) > > Best Regards, > Michael > | Best regards -- / Raimo Niskanen, Erlang/OTP, Ericsson AB _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Richard A. O'Keefe-2
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 10-bit > 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 left-closed 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 one-liner 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 half-closed 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/OTP-20.0...RaimoNiskanen:raimo/stdlib/rand-uniformNZ 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 _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
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 pre-20 (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 pre-20-equivalent 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/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/random Python random.uniform(a, b): [a, b] (when a <= b) https://stackoverflow.com/questions/6088077/how-to-get-a-random-number-between-a-float-range 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://ruby-doc.org/core-2.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/R-manual/R-devel/library/stats/html/Uniform.html Note: "runif will not generate either of the extreme values unless max = min or max-min 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/mathematical-functions.html#function_rand PostgreSQL 10 random(): [0.0, 1.0) https://www.postgresql.org/docs/10/static/functions-math.html#functions-math-random-table MS SQL Server: (0.0, 1.0) https://docs.microsoft.com/en-us/sql/t-sql/functions/rand-transact-sql "Returns a pseudo-random float value from 0 through 1, exclusive." dSFMT: "[1, 2), [0, 1), (0, 1] and (0, 1)" http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT (dSFMT: SFMT PRNG implementation for double-precision floating point only) It took me an hour to investigate this. Lesson learned: don't take the range definitions for granted. 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: _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
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 pre-20 (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 pre-20-equivalent code should be documented. If you mean pre-20 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 pre-20 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 pre-20 get exactly 0.0 and the same on 20.0 and later. You can, of course, only use the explicit algorithms that existed pre-20 when executing pre-20, 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 pre-20 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. > > 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/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/random > > Python random.uniform(a, b): [a, b] (when a <= b) > https://stackoverflow.com/questions/6088077/how-to-get-a-random-number-between-a-float-range There is a note that the upper bound may depending on rounding not be returned. Really nasty. > > 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://ruby-doc.org/core-2.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/R-manual/R-devel/library/stats/html/Uniform.html > Note: "runif will not generate either of the extreme values > unless max = min or max-min is small compared to min, > and in particular not for the default arguments." That is a really sloppy definition! > > MySQL 5.7 RAND(): [0.0, 1.0) > https://dev.mysql.com/doc/refman/5.7/en/mathematical-functions.html#function_rand > > PostgreSQL 10 random(): [0.0, 1.0) > https://www.postgresql.org/docs/10/static/functions-math.html#functions-math-random-table > > MS SQL Server: (0.0, 1.0) > https://docs.microsoft.com/en-us/sql/t-sql/functions/rand-transact-sql > "Returns a pseudo-random float value from 0 through 1, exclusive." > > dSFMT: "[1, 2), [0, 1), (0, 1] and (0, 1)" > http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT > (dSFMT: SFMT PRNG implementation for double-precision 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 10-bit > > > 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 left-closed 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 one-liner 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 > > half-closed 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/OTP-20.0... > > RaimoNiskanen:raimo/stdlib/rand-uniformNZ > > > > 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 > > _______________________________________________ > > erlang-questions mailing list > > [hidden email] > > http://erlang.org/mailman/listinfo/erlang-questions > > > _______________________________________________ > erlang-questions mailing list > [hidden email] > http://erlang.org/mailman/listinfo/erlang-questions -- / Raimo Niskanen, Erlang/OTP, Ericsson AB _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Kenji Rikitake-4
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/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/random > > Python random.uniform(a, b): [a, b] (when a <= b) > https://stackoverflow.com/questions/6088077/how-to-get-a-random-number-between-a-float-range 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 semi-open 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 + (b-a) * self.random() Now of course in exact real arithmetic, if 0 <= u < 1 and a < b then a <= (b-a)*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://ruby-doc.org/core-2.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/R-manual/R-devel/library/stats/html/Uniform.html > Note: "runif will not generate either of the extreme values > unless max = min or max-min 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/mathematical-functions.html#function_rand > > PostgreSQL 10 random(): [0.0, 1.0) > https://www.postgresql.org/docs/10/static/functions-math.html#functions-math-random-table > > MS SQL Server: (0.0, 1.0) > https://docs.microsoft.com/en-us/sql/t-sql/functions/rand-transact-sql > "Returns a pseudo-random 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.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT > (dSFMT: SFMT PRNG implementation for double-precision 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] _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
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 Box-Muller. 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/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/random > > > > Python random.uniform(a, b): [a, b] (when a <= b) > > https://stackoverflow.com/questions/6088077/how-to-get-a-random-number-between-a-float-range > > 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 semi-open 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 + (b-a) * self.random() > > Now of course in exact real arithmetic, if 0 <= u < 1 > and a < b then a <= (b-a)*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://ruby-doc.org/core-2.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/R-manual/R-devel/library/stats/html/Uniform.html > > Note: "runif will not generate either of the extreme values > > unless max = min or max-min 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/mathematical-functions.html#function_rand > > > > PostgreSQL 10 random(): [0.0, 1.0) > > https://www.postgresql.org/docs/10/static/functions-math.html#functions-math-random-table > > > > MS SQL Server: (0.0, 1.0) > > https://docs.microsoft.com/en-us/sql/t-sql/functions/rand-transact-sql > > "Returns a pseudo-random 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.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT > > (dSFMT: SFMT PRNG implementation for double-precision 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] > _______________________________________________ > erlang-questions mailing list > [hidden email] > http://erlang.org/mailman/listinfo/erlang-questions -- / Raimo Niskanen, Erlang/OTP, Ericsson AB _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Raimo Niskanen-2
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 _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
In reply to this post by Raimo Niskanen-2
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 Box-Muller. > > 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 _______________________________________________ erlang-questions mailing list [hidden email] http://erlang.org/mailman/listinfo/erlang-questions |
Free forum by Nabble | Edit this page |