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.

How Many Prisoners?

How many prisoners?

This riddle comes from the XKCD forums, at http://forums.xkcd.com/viewtopic.php?f=3&t=70558 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,
$$\begin{align}
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
\end{cases}
\end{align}
$$
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
\end{cases}
$$
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
$$\begin{align}
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
\end{align}
$$
Do the easier integral $I_2$ first via the substitution $u=\cos\theta$, $du=-\sin\theta d\theta$.
$$\begin{align}
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})
\end{align}
$$
Evaluate $I_1$ using the substitution $v=\sin\theta$, $dv=\cos\theta d\theta$, and the identity $\cos^2\theta=1-\sin^2\theta$:
$$\begin{align}
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}
\end{align}
$$
Split the integrand using $1-v^2=(1+v)(1-v)$ and the method of partial fractions
$$\begin{align}
\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}
\end{align}
$$
for some to be determined constants $A$, $B$, $C$, and $D$. Clearing fractions yields the identity
$$1=A(1-v)(1+v)^2+B(1+v)^2+C(1-v)^2(1+v)+D(1-v)^2$$
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.
$$\begin{align}
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
\end{align}
$$
Lines 1 and 2 give $B=D=1/4$. Plugging into lines 3 and 4 gives
$$\begin{align}
0 &= \frac{3}{2} – 9A + 3C \\
0 &= \frac{3}{2} + 3A – 9C
\end{align}
$$
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
$$\begin{align}
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)
\end{align}
$$
Combining $I_1$ and $I_2$ gives the final answer
$$\begin{align}
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}
\end{align}
$$
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 http://tryfs.net/).