# Is it really possible for rand to throw a zero?

**URL:** <https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505>\
**Category:** Statistics\
**Tags:** question\
**Created:** [September 7, 2019, 4:29pm UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505 "2019-09-07T16:29:09Z")\
**Posts on this page:** 20\
**Page:** 3

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 9, 2019, 7:51pm UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/41 "2019-09-09T19:51:35Z")

</div>

> [@rfourquet](#):
>
> One pretty efficient solution to generate in `(0, 1)` is to set the least significant bit to `1` (which is always set to `0` by default)

I think that this would also make the expected value exactly equal to 0.5, rather than 0.49999999999999988897…

It should be possible to do in zero clock cycles by modifying the last step of the radom-number generator to subtract `prevfloat(1.0)` instead of subtracting `1.0`.

---

<div class="post-metadata">

**Author:** ![aerdely](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerdely/32/37506_2.png) [@aerdely](https://discourse.julialang.org/u/aerdely)\
**Post date:** [September 9, 2019, 10:37pm UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/42 "2019-09-09T22:37:13Z")

</div>

Could you please provide examples of the advantage with the half open interval?

---

<div class="post-metadata">

**Author:** ![PeterB](https://avatars.discourse-cdn.com/v4/letter/p/e8c25b/32.png) [@PeterB](https://discourse.julialang.org/u/PeterB)\
**Post date:** [September 11, 2019, 6:22am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/43 "2019-09-11T06:22:28Z")

</div>

(Answering @aerdely) If you want an event to happen with probability p:

if rand() \< p  
do\_something  
end

That way do\_something never happens if p==0 and always happens if p==1. It certainly makes my code easier!

Historically, I suspect the reason goes back (at least) to the K&R C rand() function, which (on a 32 bit system) returns an integer on [0, 2^32 - 1] (for hardware reasons), and then dividing by 2^32 is the simple way to get a floating point number on [0, 1).

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 11, 2019, 9:21am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/44 "2019-09-11T09:21:58Z")

</div>

The more I’ve thought about this, the more I’ve come to believe that the current Julia (and C, C++, Fortran, Python, Perl and Ruby) implementations are wrong. Here’s why:

The `rand` function does not actually draw from the infinite set of points on the (open or closed) interval. It draws from the finite set `0:eps(T):1-eps(T)`. (Which has the convenient property that the elements can be represented in type `T`.)

We might want to think about `rand` as drawing from the infinite set (the open or closed interval) and then projecting (rounding) onto a finite set of floating-point numbers. But for that metaphor to work with the above-mentioned set, the interval must be `(-eps(T)/2, 1-eps(T)/2)`. (The open or the closed interval. Mathematically, there’s no difference. The probability of hitting an endpoint in a finite number of draws is zero.)

Setting the least significant bit to 1, as suggested by @rfourquet would be equiivalent to drawing from `eps(T)/2:eps(T):1-eps(T)/2`, which in turn would be equivalent to drawing from the (open or closed) continuous interval (0,1) and then correctly rounding to the nearest member of that set. This would not cost any extra time, and is the correct fix, i.m.o.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 11, 2019, 9:32am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/45 "2019-09-11T09:32:52Z")

</div>

> [@Per](#):
>
> We might want to think about `rand` as drawing from the infinite set (the open or closed interval) and then projecting (rounding) onto a finite set of floating-point numbers.

I am not sure we want that, as that would be very tricky to implement (note that floating point numbers get denser near 0).

> [@Per](#):
>
> the current Julia (and C, C++, Fortran, Python, Perl and Ruby) implementations are wrong

These implementations just use a reasonable compromise between speed and “correctness”, approximating the uniform CDF with a step function. It is of course understood that there exist _floating point_ numbers between 0 and 1 which can _never_ be drawn by `rand`, eg `eps()^4`.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 11, 2019, 9:38am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/46 "2019-09-11T09:38:27Z")

</div>

Also, it is not even `eps(T)`. If you want to see what numbers can be generated, try something like

```julia
julia> z(n) = reinterpret(Float64, 0x3ff0000000000000 | UInt64(n)) - 1
z (generic function with 1 method)

julia> z(0)
0.0

julia> z(2^52-1)
0.9999999999999998

```

_edit_ Thinking of this, actually correcting `rand(Float64)` by `2^(-53)` would be a reasonable idea: it would give numbers in (0, 1), and remove the bias.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 11, 2019, 9:46am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/47 "2019-09-11T09:46:59Z")

</div>

> [@Tamas\_Papp](#):
>
> that would be very tricky to implement

Sorry if I was unclear. By “a set of floating point numbers” I meant a set like `0:eps(T):1-eps(T)` as previously mentioned. I did not mean “the set of all floating point numbers of a given type”.

> [@Tamas\_Papp](#):
>
> it is not even `eps(T)`

```julia
julia> z(1) - z(0) === eps(Float64)
true

```

Yes it is.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 11, 2019, 9:50am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/48 "2019-09-11T09:50:20Z")

</div>

Good point, thanks for the correction.

I am wondering whether to open an issue about just adding `eps(T)/2` to `rand`. It would remove the bias, and solve the issue in this topic.

_edit_ now with your clarification I see this is the same thing you proposed originally.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 11, 2019, 10:02am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/49 "2019-09-11T10:02:45Z")

</div>

Please do! From the point of view of Monte Carlo simulations etc, this should only make things better, and since it can be accomplished in zero clock cycles, there’s no performance argument.

I think the main argument against would be that some people might have tests which rely on random numbers being predictable, like:

```julia
@test my_function(rand(100)) == 1.87554603778140

```

The solution is then to fix those tests.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 11, 2019, 10:17am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/50 "2019-09-11T10:17:04Z")

</div>

If this gets implemented, then another reason to prefer Julia over Python will be that “Julia’s random numbers are bigger!”

🙂

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 11, 2019, 10:24am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/51 "2019-09-11T10:24:31Z")

</div>

[https://github.com/JuliaLang/julia/issues/33222](https://github.com/JuliaLang/julia/issues/33222)

Incidentally, this would make all random floats interior, so the confusion in this topic would never arise.

---

<div class="post-metadata">

**Author:** ![aerdely](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aerdely/32/37506_2.png) [@aerdely](https://discourse.julialang.org/u/aerdely)\
**Post date:** [September 13, 2019, 4:49am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/52 "2019-09-13T04:49:17Z")

</div>

In the comments of  
[https://github.com/JuliaStdlibs/Random.jl/blob/master/src/normal.jl](https://github.com/JuliaStdlibs/Random.jl/blob/master/src/normal.jl)  
it is mentioned that

```julia
# The Ziggurat Method for generating random variables - Marsaglia and Tsang
# Paper and reference code: http://www.jstatsoft.org/v05/i08/

```

is used for `randn` and `randexp` but such paper requires the use of a uniform generator in the unit open interval (0,1) and Julia’s code uses `rand` which generates from [0,1) therefore, strictly speaking, the assumption of the Ziggurat Method is being violated, and so it is that Julia’s code has to define a function `randexp_unlikely` to fix it and avoid throwing a zero. That is why, from a probabilistic point of view, it is better to have uniform random generators in the open unit interval (0,1).

Since Julia’s code for `randexp` avoids throwing a zero, it is straightforward to obtain random variates from the open interval (0,1) because the distribution of \exp(-X) is Uniform in the open interval (0,1) whenever X has exponential distribution with scale parameter equal to 1.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 13, 2019, 7:16am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/53 "2019-09-13T07:16:05Z")

</div>

The reason for `randexp_unlikely` is not that `rand(Float64)` might return an exact zero. According to the comments, `randexp_unlilkely` gets called about 1.2 % of the time, while an exact zero only happens about 0.0000000000002 % of the time.

It seems that in extremely rare cases (less than once per 10^17 draws) `randexp` will return `Inf`, which is a bug. (A 1000-hour Monte-Carlo simulation would have a 1 % chance of returning Inf or NaN instead of the desired result.)

[Edit: I think I may have overestimated the probability. Maybe it’s more like once in 10^19 draws, or a 0.01 % chance per 1000 CPU-hours of simulation. I don’t have access to a compute cluster big enough that I can test…]

With the fix to `rand` proposed above, `randexp` would instead return ~~36.7368005696771~~ 44.43391803980815 in those cases.

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 13, 2019, 10:51am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/54 "2019-09-13T10:51:47Z")

</div>

For what it’s worth, I investigated what some other languages are doing:

In Python and Ruby (at least the builds that happened to be installed on my computer) there seems to be one more bit of randomness, so the bias is -eps/4.

Perl seems to do exactly the same thing as Julia.

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [September 15, 2019, 10:29pm UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/55 "2019-09-15T22:29:49Z")

</div>

A highly redacted and rewritten answer to this question exists on [stackoverflow](https://stackoverflow.com/questions/57948193/julia-is-it-really-possible-for-rand-to-throw-a-zero) as an “asked and answered” entry.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 16, 2019, 6:07am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/56 "2019-09-16T06:07:11Z")

</div>

I made a PR for this, currently waiting for review:

[https://github.com/JuliaLang/julia/pull/33251](https://github.com/JuliaLang/julia/pull/33251)

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 16, 2019, 7:01am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/57 "2019-09-16T07:01:49Z")

</div>

I made a T-shirt, in case the PR gets accepted. 😉

 ![t-shirt](https://global.discourse-cdn.com/julialang/original/3X/9/a/9a0169e1866cb99f2a0aca7e9775fe790510ea13.jpeg)

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [September 16, 2019, 7:40am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/58 "2019-09-16T07:40:19Z")

</div>

> [@Per](#):
>
> The more I’ve thought about this, the more I’ve come to believe that the current Julia (and C, C++, Fortran, Python, Perl and Ruby) implementations are wrong.

To summarize the issue again: Julia and most frameworks sample a uniform fixed-point number in (0,1] and convert it to floating point (no rounding). One can argue that the proper way of generating random floats would be to generate a uniform infinite-precision number and round it to float. That is very different because floats are more precise than fixed-point numbers close to zero.

The “improper” way has probability of `nextfloat(T(1))-T(1)` of emitting zero, while the “proper” way has probability `nextfloat(T(0))`. The proper way has two more bit of entropy.

To put numbers on it: `1.2f-7` vs `1.0f-45` in `Float32` and `2e-16` vs `5e-324` for `Float64`. We see that “proper” draws should give zero or subnormals approximately never, and Float32 draws a lot of arguably spurious zeros.

For reasonably fast generation of “almost proper” random floats, see [Conversion to Float · Issue #8 · JuliaRandom/RandomNumbers.jl · GitHub](https://github.com/sunoru/RandomNumbers.jl/issues/8#issuecomment-494339113)

However, that would be slower as long as we use a weird RNG like Mersenne twister (make 54 random bits in correct position) instead of “random bitstring”-style RNGs (make UInt32 / UInt64 / UInt128 / UInt256 only, and convert later).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 16, 2019, 8:10am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/59 "2019-09-16T08:10:02Z")

</div>

> [@foobar\_lv2](#):
>
> One can argue that the proper way of generating random floats would be to generate a uniform infinite-precision number and round it to float.

I don’t think @Per is arguing for that (I also misunderstood this originally).

That said, my inner geek would find it fascinating to come up with a scheme for this (eg first draw and exponent with a geometric distribution, then fill in the bits with a 52-bit uniform, there must be some wrinkles one would need to think of), but in practice if one needs “proper” randomness around 0 then something like a transformed value with `randexp` is the way to go (which, in turn, is using a very clever algorithm along these lines).

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [September 16, 2019, 9:57am UTC](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505/60 "2019-09-16T09:57:37Z")

</div>

Yes, to be clear: I think 51 bits of randomness is absolutely fine as the default option. I was only arguing for removing the bias, as in the PR.

As a side note, and unrelated to the bias issue, I wouldn’t mind if there were an option to get what’s referred to above as “proper” draws. This can easily be accomplished as `rand(T) + eps(T)*rand(T) + eps(T)^2*rand(T) + ` … where the series can be truncated after two terms in almost all cases, so the average runtime is approximately twice that of `rand(T)`.

One place where one might want to use such draws is when testing the precision of floating-point code. Otherwise one might miss bugs such as using `log(1+x)` instead of `log1p(x)` that only manifest themselves when the input is not a multiple of `eps(T)˜.

[Previous page](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505.md?page=2)

[Next page](https://discourse.julialang.org/t/is-it-really-possible-for-rand-to-throw-a-zero/28505.md?page=4)