Generating Uniform Random Numbers on [0,1]

Generating Uniform Random Numbers on [0,1]

Chris Lomont, March 2017

Common advice for generating uniform random numbers in [0,1] in many languages looks like this [1,2,3]:

This however does not generate many possible floating-point numbers in $[0,1]$, leaving gaps – there are many floating-point numbers that should occur but cannot due to this method. To understand why, you need to understand how floating-point numbers are (usually) stored on computers.

Throughout we’ll assume the random number generators produce uniformly random integers in a known range.

Efficient divisibility testing

Efficient divisibility testing

The example

Here is a low computational cost method to test if an integer $x$ is divisible by another integer $d$, useful in C-like languages (assembly included). It replaces a possibly costly modulus % operation with a usually lower cost multiplication *.

Here is an example illustrating the method for divisibility by 5:

What trickery is this? How does it work? Where do the constants 3435973837U and 858993459U come from?

Note the magic numbers depend on the type of $x$. The above example requires $x$ is an unsigned 32 bit integer, and the language has multiplication overflow computed modulus 32. Many common languages (C,C++,C#,Java, ..) do this.

Making Three NOT Gates from Two

Making Three NOT Gates from Two

Here is an old and interesting logic puzzle.


Recall an AND gate takes in two inputs A and B, and has one output X. The truth table looks like

A AND B 0 1
0 0 0
1 0 1

Similarly, OR looks like

A OR B 0 1
0 0 1
1 1 1

And NOT simply inverts its input:

0 1
1 0

LZCL: An embedded compression library

LZCL: An embedded compression library

Copyright Chris Lomont, Sept 2016
Version 1.0


See my Compression Notes article for compression background on useful for this article. The README at the GitHub page has usage details .

For a recent project, I needed to design a low resource compression scheme. Over the years I’ve designed many, and this time I decided to unify some code, clean it up, and make some tools to ease development of such algorithms (mostly for my own use). The result this particular time is several very low memory decoders, necessarily slow, including Huffman, Arithmetic, LZ77, and my own new invention, LZCL, which outperforms the others quite well. All are designed for minimal RAM usage.

How Many Prisoners?

How many prisoners?

This riddle comes from the XKCD forums, at and is as follows:

You, together with a finite number n-1 of other ideal mathematicians, have been arrested on a whim by a generic evil dictator and are about to be locked up in a prison. The prison is circular, with n identical windowless cells arranged in a ring around a central court. There are some problems with the lighting system – the light switch in each cell controls the light in the next cell clockwise around the ring. Even worse, electric power is only provided to the lights for one tenth of a second each night, just after midnight.

Pick a Number

Pick a Number

Here is a logic puzzle I find interesting. Unfortunately I don’t recall the source.

Alice picks two distinct real numbers, and puts one in each of two envelopes. Bob wants to pick the envelope containing the larger of the two numbers. Bob picks an envelope, and look at the number inside it.

After seeing the number, Bob can elect to switch envelopes. If he switches he has to take the second. If he doesn’t switch, he has to keep the first.

Is there a strategy that allows Bob to select the larger value with greater than a 50% chance?

Average Segment Length

Average segment length

Here is an elementary derivation of the average length of a segment chosen randomly on the unit square. By random I mean $x$ and $y$ coordinates are chosen uniformly in $[0,1]$. Let the endpoints be $(x_1,y_1)$ and $(x_2,y_2)$. Then the average is given by
$$\int_0^1 \int_0^1 \int_0^1 \int_0^1 \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2} dx_1 dx_2 dy_1 dy_2 $$
which unfortunately is difficult to evaluate directly, so I’ll use the deltas $\Delta x = x_1-x_2$ and $\Delta y = y_1-y_2$ instead. Now I need the probability distributions on the $\Delta$s to evaluate the resulting integral.
The distribution of the difference of two uniform random variables is the triangle distribution. Here is a quick derivation:
Let $X_1$ and $X_2$ be uniform random variables on $[0,1]$. They have probability distribution functions (PDFs) of $f_{X_i}(x)=1$ and cumulative probability distribution functions (CDFs) of $F_{X_i}(x)=x$. The joint PDF of $X_1$ and $X_2$ is
$$f_{X_1,X_2}(x_1,x_2)=1 \qquad 0\le X_1,X_2 \le 1$$
Let $Z=X_1-X_2$ be the random variable of the difference. Using the method of cumulative probability functions,
F_Z(z) &= P(Z\le z) \\
&= P(X_1-X_2\le z) \\
&= \begin {cases}
\int_0^{1+z} \int_{x_1-z}^1 f_{X_1,X_2}(x_1,x_2) dx_2 dx_1, & -1\le z < 0 \\
1-\int_{z}^{1} \int_0^{x_1-z} f_{X_1,X_2}(x_1,x_2) dx_2 dx_1, & 0\le z \le 1
\end{cases} \\
&= \begin {cases}
z^2/2 + z + 1/2, & -1\le z < 0 \\
-z^2/2 + z + 1/2, & 0\le z \le 1
Differentiating w.r.t $z$ then gives the PDF
f_Z(z) =
\begin {cases}
1+z, & -1\le z < 0 \\
1-z, & 0\le z \le 1
which simplifies to $f_Z(z)=1-|z|$ for $-1\le z \le 1$.
To evaluate the average length, let $x=x_1-x_2$ and $y=y_1-y_2$, then the average length of a random segment is given by
$$\int_{-1}^1\int_{-1}^1 \sqrt{x^2+y^2}(1-|x|)(1-|y|) dx dy$$
Over the integration region of the square $[-1,1]\times[-1,1]$ this is symmetric over the four quadrants, so integrating over $[0,1]\times[0,1]$ obtains the value $I$ as
$$I=4\int_0^1\int_0^1 \sqrt{x^2+y^2}(1-x)(1-y) dx dy$$
Switch to polar coordinates, $x=r \cos\theta$, $y=r\sin\theta$, $\sqrt{x^2+y^2}=r$, using $xy$ symmetry, integrate over the region $0\le\theta\le\pi/4$ and $0\le r\le \sec\theta$. Note the right side of the square, $x=1$, gives the $r$ bound $x=r\cos\theta=1$, i.e., $r=1/\cos\theta=\sec\theta$. Remembering the Jacobian for the change of variables gives $dx dy=r dr d\theta$. The symmetry gives another factor of 2 in front, resulting in
I &= 2\times 4 \int_0^{\pi\over 4} \int_0^{\sec\theta} r (1-r\cos\theta)(1-r\sin\theta) r dr d\theta \\
&= 8 \int_0^{\pi\over 4} \int_0^{\sec\theta} r^2(1-r\cos\theta-r\sin\theta+r^2\sin\theta\cos\theta) dr d\theta \\
&= 8 \int_0^{\pi\over 4} \left. \frac{r^3}{3}-\frac{r^4}{4}\cos\theta-\frac{r^4}{4}\sin\theta+\frac{r^5}{5}\sin\theta\cos\theta \right|_0^{\frac{1}{\cos\theta}} d\theta \\
&= 8 \int \frac{1}{12\cos^3\theta} – \frac{\sin\theta}{20\cos^4\theta}d\theta \\
&= \frac{2}{3}\int \frac{1}{\cos^3\theta} d\theta + \frac{2}{5}\int \frac{-\sin\theta}{\cos^4\theta}d\theta \\
&= I_1 + I_2
Do the easier integral $I_2$ first via the substitution $u=\cos\theta$, $du=-\sin\theta d\theta$.
I_2 &= \frac{2}{5}\int_1^{\frac{1}{\sqrt{2}}} \frac{du}{u^4} \\
&= \frac{2}{5}\left.\frac{u^{-3}}{-3}\right|_1^{\frac{1}{\sqrt{2}}} \\
&= \frac{2}{15}(1-2\sqrt{2})
Evaluate $I_1$ using the substitution $v=\sin\theta$, $dv=\cos\theta d\theta$, and the identity $\cos^2\theta=1-\sin^2\theta$:
I_1 &= \frac{2}{3}\int_0^{\pi\over 4} \frac{d\theta}{\cos^3\theta} \\
&= \frac{2}{3}\int_0^{\pi\over 4} \frac{\cos\theta d\theta}{\cos^4\theta} \\
&= \frac{2}{3}\int_0^{1\over \sqrt{2}} \frac{dv}{(1-v^2)^2}
Split the integrand using $1-v^2=(1+v)(1-v)$ and the method of partial fractions
\frac{1}{(1-v^2)^2} &= \frac{1}{(1-v)^2(1+v)^2} \\
&= \frac{A}{1-v}+\frac{B}{(1-v)^2}+\frac{C}{1+v}+\frac{D}{(1+v)^2}
for some to be determined constants $A$, $B$, $C$, and $D$. Clearing fractions yields the identity
This can be expanded and powers of $v$ set equal on each side, and the resulting linear system solved. A useful trick here is to plug in values for $v$ to get enough relations.
v &=+1 &\implies& 1= 0A + 4B + 0C + 0D \\
v &=-1 &\implies& 1= 0A + 0B + 0C + 4D \\
v &=+2 &\implies& 1=-9A + 9B + 3C + D \\
v &=-2 &\implies& 1= 3A + B – 9C + 9D
Lines 1 and 2 give $B=D=1/4$. Plugging into lines 3 and 4 gives
0 &= \frac{3}{2} – 9A + 3C \\
0 &= \frac{3}{2} + 3A – 9C
Subtracting line 1 from line 2 gives $0=12A-12C$ so $A=C$. Plugging into line 1 and solving gives $A=C=1/4$.
Plugging into $I_1$ and pulling the $1/4$ out front gives
I_1 &= \frac{1}{6}\int_0^{1\over\sqrt{2}} (1+v)^{-2}+(1-v)^{-2}+(1+v)^{-1}+(1-v)^{-1} dv \\
&= \frac{1}{6} \left(\left. \frac{(1+v)^{-1}}{-1} – \frac{(1-v)^{-1}}{-1} + \ln|1+v| – \ln|1-v| \right|_0^{1\over\sqrt{2}}\right) \\
&= \frac{1}{6}\left(\frac{1}{1-{1\over\sqrt{2}}}-\frac{1}{1+{1\over\sqrt{2}}} -1 + 1 + \ln\frac{1+{1\over\sqrt{2}}}{1-{1\over\sqrt{2}}}-\ln\frac{1}{1} \right) \\
&= \frac{\sqrt{2}}{3} + \frac{1}{6} \ln\left(3+2\sqrt{2}\right)
Combining $I_1$ and $I_2$ gives the final answer
I &= I_1+I_2\\
&= \frac{\sqrt{2}}{3} + \frac{1}{6} \ln\left(3+2\sqrt{2}\right) + \frac{2}{15}(1-2\sqrt{2}) \\
&= \frac{2+\sqrt{2}}{15} +\frac{\ln(3+2\sqrt{2})}{6}
As a sanity check simulate it. Note $$\frac{2+\sqrt{2}}{15} +\frac{\ln(3+2\sqrt{2})}{6} = 0.5214….$$
Here is a simulation in F#:

returns 0.520409, agreeing with the symbolic result. (It can be tested online, for example, at

Initial post – testing features

Initial testing post
Some inline math $E=mc^2$.

Markdown testing




This italic word bold word. Here is some of both.

horiz rule is three or more —
Here is strikethrough

  1. List
  2. More
    • unordered
    • more unordered
  3. Final

Inline link for google
Inline code backticked command code.

Tables Are Cool
left center right

block quotes go here after initial >