I ran across the Cliff random number generator yesterday. Given a starting value x0 in the open interval (0, 1), the generator proceeds by
xn+1 = | 100 log(xn) mod 1 |
for n > 0. The article linked to above says that this generator passes a test of randomness based on generating points on a sphere.
Real numbers
While the long term distribution of the generator may be good, it has a problem with its sequential behavior, namely that it has an infinite number of fixed points. If the generator ever reaches one of these points, it gets stuck forever.
Here’s a proof. Since x is between 0 and 1, log(x) is always negative. So we can replace the absolute value above with a negative. A fixed point of the generator is a solution to
x = -100 log(x) – k
for some integer k. Define
f(x, k) = -100 log(x) – x – k.
For each non-negative k, the limit of f(x, k) as x goes to 0 is ∞ and the limit as x goes to 1 is negative, so somewhere in between it is zero.
Floating point numbers
So in exact operations over the reals, there is a fixed point for each non-negative integer k. However, when implemented in finite precision arithmetic, it manages to get itself unstuck as the following Python code shows with k arbitrarily chosen to be 20.
from numpy import log
from scipy.optimize import bisect
r = bisect(lambda x: -100*log(x) - x - 20, 0.4, 0.999)
print(r)
for i in range(10):
r = abs(-100*log(r)) % 1
print(r)
This produces
0.8121086949937715
0.8121086949252678
0.8121087033605505
0.8121076646716787
0.8122355649800639
0.7964876237426814
0.7543688054472710
0.1873898667258977
0.4563983607272064
0.4389252746489802
0.3426097596058071
In infinite precision, r above would be a fixed point. In floating point precision, r is not. But it does start out nearly fixed. The first iteration only changes r in the 11th decimal place, and it doesn’t move far for the next few iterations.
Update: See the next post for how the random number generator does on the DIEHARDER test suite.
Plotting fixed points
The kth fixed point is the solution to f(x, k) = 0. The following code plots the fixed points as a function of k.
t = arange(300)
y = [bisect(
lambda x: -100*log(x)-x-k,
0.001,
0.999)
for k in t]
plt.plot(t, y)
plt.xlabel("k")
plt.ylabel("fixed point")
The fixed points cluster toward zero, or they would in infinite precision arithmetic. As we showed above, the Cliff random number generator performs better in practice than in theory. Maybe the generator does well after wandering close to zero, but I wouldn’t be surprised if it has a bias toward the low end of the interval.