# Recursive random number generator based on irrational numbers

Here $$\{\cdot\}$$ and $$\lfloor \cdot\rfloor$$ denote the fractional part and floor functions respectively. For a negative, non-integer number $$x$$, we use the following definition: $$\{x\}=1-\{-x\}$$. If $$x$$ is a negative integer, $$\{x\} =0$$. We are dealing with the following recurrence:

$$X_{k+2}=\{b_2 X_{k+1}+b_1 X_k\}$$

where $$X_1$$ is a uniform random variable on $$[0,1]$$ and $$X_0\in [0,1]$$ is a constant. Thus all the $$X_k$$'s are in $$[0,1]$$. Also, $$b_1, b_2$$ are integers, called bases; they represent bases in a numeration system.

The simple case: $$b_1=0$$

I extensively studied the case $$b_1=0, b_2 > 1$$ corresponding to a first-order recurrence, see here. The main results are:

• The sequence $$\lfloor b_2X_k \rfloor$$ corresponds to the digits of $$X_1$$ in base $$b_2$$. These digits behave as independently and identically distributed discrete uniform variables on $$\{0, 1,\cdots,b-1\}$$.

• The sequence $$X_k$$ behaves as identically distributed continuous uniform variables on $$[0, 1]$$. The correlation between $$X_k$$ and $$X_{k+m}$$ is equal to $$b_2^{-m}$$.

• For a specific value of $$X_1$$, say $$X_1=c$$ with $$c$$ a normal number (say $$c=\log 2$$), the empirical process of observed $$X_k$$'s (corresponding to a specific realization of the theoretical stochastic process) satisfies the same properties for the empirical statistics: convergence of the empirical distribution to uniform on $$[0, 1]$$, convergence of the empirical auto-correlations to the theoretical values mentioned above, etc.

This happens because the sequence is ergodic. Note that almost all numbers are normal, though no one knows if any of $$e,\pi,\sqrt{2},\log 2$$ is normal. They are believed to be normal.

The general case, and my question

The general case is when both $$b_1$$ and $$b_2$$ are non zero. For simplicity, we can focus on the following specific case, which seems to behave very nicely: $$X_0=\frac{1}{2}, b_1=-3, b_2 = -5$$. More specifically, it now looks like the $$X_k$$'s are not only uniformly distributed on $$[0, 1]$$, but also asymptotically independently distributed. Thus we can use that sequence as a random number generator, with $$X_1$$ being the seed. This is a big contrast with the simple case discussed in the first section.

For instance (this is an illustration of what I mean by asymptotic independence), if $$X_1=\frac{\sqrt{2}}{2}$$, the empirical probabilities satisfy

$$\hat{P}\Big[\bigcap_{i=0}^m (X_{k+i}<\alpha_i)\Big]\rightarrow \prod_{i=0}^m \hat{P}\Big[X_{k+i}<\alpha_i\Big]\rightarrow\prod_{i=0}^m \alpha_i$$

regardless of $$m$$ and $$0\leq \alpha_0,\cdots,\alpha_m\leq 1$$, when more and more terms (that is more and more $$k$$'s) are used to estimate these probabilities. I thus assume (maybe erroneously) that it must also be true for the theoretical probabilities. This is illustrated further in the Appendix (last section).

My question is whether my conjecture (independence of the $$X_k$$'s) is true. It was verified empirically when $$X_0=\frac{1}{2}, X_1=\frac{\sqrt{2}}{2}, b_1=-3, b_2=-5$$, as well as for many other parameter sets. The generated deviates seem to approach randomness better than those generated using Excel, based on various statistical tests. Note that not any parameter set works; there are plenty of exceptions, and identifying these exceptions would be a bonus.

Computational considerations

No need to read this section, only if you are interested, but it is not directly related to my question.

When you compute the successive $$X_k$$'s, you lose one bit of precision at each iteration. This is not an issue thanks to ergodicity, it's like re-starting the sequence with new seeds every 45 or so iterations. It's only an issue if you look e.g. at long-range auto-correlations.

Also, it is possible to carry the computations very efficiently. You start with getting and storing several billions of binary digits of $$X_1=\frac{\sqrt{2}}{2}$$. See here how you can get these digits. Then you only need to perform simple additions and bit shifting with a big number library. For instance, $$5x = 4x + x$$, and computing $$4x$$ is just a bit shifting operation (no multiplication involved). Likewise with $$3x=2x + x$$. In my Perl code, if I use this little trick, it runs 10 times faster than doing an actual multiplication.

Appendix

I estimated the probability $$P(X_k<\alpha_0, X_{k+1}<\alpha_1, X_{k+2}<\alpha_2)$$ for a thousand randomly selected triplets $$(\alpha_0,\alpha_1,\alpha_2)$$ in $$[0, 1]^3$$ and 100,000 $$(X_k,X_{k+1},X_{k+2})$$'s. Assuming uniform distribution and independence between $$X_k, X_{k+1}$$ and $$X_{k+2}$$, the theoretical value is always $$\alpha_0\cdot \alpha_1\cdot \alpha_2$$. The data and source code is available in an Excel spreadsheet, here. It is very easy to replicate my results. The observed and theoretical values are extremely close, supporting the conjecture of stochastic independence and uniformity. Below is a scatter plot where each point corresponds to one of the $$(\alpha_0, \alpha_1, \alpha_2)$$'s, with the X-axis being the observed (estimated) probability, and the Y-axis being the theoretical probability (the product of $$\alpha_0,\alpha_1,\alpha2$$).

By contrast, here is an example where the independence assumption is violated, corresponding to $$b_1=0, b_2=3$$ with $$\mbox{Correl}(X_k,X_{k+1}) = \frac{1}{3}$$ and $$\mbox{Correl}(X_k,X_{k+2}) = \frac{1}{9}$$:

In short, in the example with independence (first chart) you seem to have, assuming $$X_1$$ is uniformly distributed on $$[0,1]$$: $$P(X_{k+2}\in A_2 | X_{k+1}\in A_1, X_k\in A_0)=\mu(A_2)$$ where $$A_0, A_1, A_2$$ are Borel subsets of $$[0, 1]$$ and $$\mu(\cdot)$$ is the Lebesgue measure. Yet it is obviously true that $$X_{k+2}-b_2X_{k+1}-b_1X_k$$ can only take on $$|b_1|+|b_2|$$ distinct integer values. Note that if you consider the sequence $$Y_k=X_{3k}$$ instead, then the triplets $$(Y_k,Y_{k+1},Y_{k+2})$$ don't lie in such a small number of planes, unlike the $$X_k$$'s.

Rephrased differently, my question is whether or not the uniform distribution with independence is the main fixed point (also called equilibrium distribution or attractor) of the stochastic / chaotic system in question. Usually finding the attractor requires solving a stochastic integral equation, yet here if we suspect Uniform/Independent might be the solution, you just need to plug that presumed solution in the integral equation and see if it solves it.

Update on 9/25/2020

Based purely on pattern recognition techniques, I've found this:

$$X_k=\{A(k)X_1\}, \mbox{ with } A(k)=b_2A(k-1)+b_1A(k-2)$$

with $$A(0)=\frac{X_0}{X_1}$$ and $$A(1)=1$$. I don't have a proof, but this looks like something very easy to prove. In addition, it helps prove whether or not the stochastic uniform/independence solution is correct or not. More about this next week.

• Version 14 of this question. Sep 24 '20 at 5:24
• @ Gerry: I am only editing my answer, not sure if this counts in the 14 versions. Of course I want to minimize edits, but still have to make a lot of effort in that direction. For instance now I have found the exact number of planes for the sequence $Y_k$ in my answer, as a function of $b_1,b_2$.. Sep 24 '20 at 5:53
• No, it's 14 versions of the question. Each edit brings the question back to the front page, bumping off some other question. Each edit makes the question more of a moving target, more difficult to answer, more work for the reader to see what has really changed. Each edit may cause some readers to wonder why the author couldn't get it right the first time, or the second, or ... or the 13th. Sep 24 '20 at 6:00
• Is there a way to make the question invisible / hidden for (say) 2 weeks so you can make many edits without bothering anyone until it's finalized? If I delete it, can I still edit it in private, and then restore it 2 weeks later? I suppose MO works differently, but on my blog you can save as draft (nobody but you see it) until you decide to publish. And you can un-publish / re-publish. Wondering if there is such a feature here. Sep 24 '20 at 6:17
• The workaround I have is to edit the question, save the content in Notepad on my desktop, and not publish the edit on MO. Do that 15 times, bring the notepad version back to MO for a "private edit" without saving it on MO, and instead save the edits in Notepad each time, until after a while it's good enough that I can actually publish the final version on MO. Not ideal, but that way nobody but me will see the intermediate edits. Sep 24 '20 at 6:42

Of course the $$X_k$$ are not independent as random variables. So I assume you are referring to some notion of asymptotic independence, and it would help if you state your conjecture more precisely. One natural guess is equidistribution (see [1]) of r-tuples $$(X_k,\ldots,X_{k+r-1})$$. However, The triples $$(X_k, X_{k+1},X_{k+2})$$ will lie on a bounded number of planes in $$[0,1]^3$$; there will be at most 9 such planes in your example. This can be seen if one graphs these triples in 3D. A similar problem arose in the classical RANDU random number generator, see [3].

[1] Kuipers, L.; Niederreiter, H. (2006) [1974]. Uniform Distribution of Sequences. Dover Publications. [2] https://en.wikipedia.org/wiki/Equidistributed_sequence# [3] https://en.wikipedia.org/wiki/RANDU

• Thank you, I clarified my conjecture according to your comment. Sep 20 '20 at 0:03
• @ Yuval: I will look into that in more details. I noticed that for the simple case ($b_1=0$) you have the same problem in 2D, with couples $(X_k, X_{k+1})$ fitting on exactly $b_2$ parallel lines. If you consider the sequence $X_{2k+1}$ instead of $X_k$, the problem is gone. I would expect in the general case (3D issue) the number of planes should depend on $b_1, b_2$. The larger $b_1, b_2$, the more planes, and thus the better. But I haven't checked this yet, just a wild guess at this point. Sep 20 '20 at 2:16
• @ Yuval: how did you come up with number of at must 8 planes? Is it by visual inspection or based on some logic? That number 8 may just be $|b_1|+|b_2|$. Maybe just a coincidence, maybe not. Maybe working with the sequence $X_{3k+1}$ would fix this at least in 3D. Sep 20 '20 at 2:18
• Notice that $X_{k+1} - b_2 X_k - b_1 X_{k-1}$ is by definition an integer, and you can check that the expression can acheive at most $|b_1| + |b_2|$ different values. Sep 20 '20 at 5:37
• Suppose $b_1$ and $b_2$ are positive. As Random wrote, by definition $X_{k+2}-b_2 X_{k+1}−b_1X_k$ is an integer between $0$ and $- b_1-b_2$. So that actually gives up to 9 planes in your example. No triple can lie outside these planes. Sep 21 '20 at 17:06

Here $$b_1, b_2>0$$ are integers. I investigated the sequence $$Y_k=X_{3k}$$, which has far more communal planes, and thus more useful to build a random generator. Of course, choosing large values for $$b_1,b_2$$ will further drastically improve the generator by adding a lot more planes. I suggest choosing values larger than (say) $$2^{30}$$ for $$b_1,b_2$$.

There are $$M=b_2^3+3b_1b_2+b_1^3$$ communal planes and they all have an equation of the form

$$b_1^3\cdot Y_k+b_2(b_2^2+3b_1) \cdot Y_{k+1}-Y_{k+2} =d.$$

The possible values for $$d$$ are $$0,1,\cdots,M-1.$$ Each plane (identified by $$d$$) contains a different proportion of triplets $$(Y_k,Y_{k+1},Y_{k+2})$$. The empirical distribution for these proportions is featured in the histogram below (corresponding to $$b_1=5,b_2=3$$), where the X-axis represents $$d$$, and the Y-axis the proportion of triplets lying in plan $$d$$.

Of course it is easy by looking at this chart to guess what the exact theoretical distribution is. To identify these planes, I used the program below and some experimental math.

# Compute equations of planes containing 3 random vectors
#   P(k) = (x[k], x[k+1], x[k+2])
#   P(l) = (x[l], x[l+1], x[l+2])
#   P(m) = (x[m], x[m+1], x[m+2])
# (k, l, m) are randomly selected (M triplets)
#
# Equation of planes is x + s*y + t*z = intercept
# For each (k,l,m) output the coefficients s, t, intercept
#
# Goal: Find communal planes absorbing many (P(k), P(l), P(m))
# Once the planes are computed, sort them by s, t, intercept

$n=100000; $$b1=5;$$b2=3; # xx[] is the original sequence $$xx[0]=0.5;$$xx[1]=sqrt(2)/2; for ($$k=2;$$k<$$n;$$k++) { $$xx[$$k]=$$b2*$$xx[$$k-1]+$$b1*$$xx[$$k-2]-int($$b2*$$xx[$$k-1]+$$b1*$$xx[$$k-2]); if ($$xx[$$k]<0) { $$xx[$$k]=1+$$xx[$$k]; } } # we actually use 1 out of 3 consecutive terms from original sequence xx[] # to see if it the new sequence x[] also has a small number of communal planes for ($$k=0;$$k< $$n/3;$$k++) { $$x[$$k]=$$xx[3*$$k]; }$M=10000; # must be < n/3
open(OUT,">coplanes2.txt");

for ($$iter=0;$$iter<$$M;$$iter++) {

$$k=int($$M*rand());
$$l=int($$M*rand());
$$m=int($$M*rand());

# in case k=l or k=m or l=m, an ERROR message is reported

$$a=$$x[$$k];$$b=$$x[$$k+1]; $$c=$$x[$$k+2];$$d=$$x[$$l]; $$e=$$x[$$l+1];$$f=$$x[$$l+2];
$$p=$$x[$$m];$$q=$$x[$$m+1]; $$r=$$x[$$m+2];$$u=($$e-$$b)*($$r-$$c)-($$f-$$c)*($$q-$$b);
$$v=-($$d-$$a)*($$r-$$c)+($$f-$$c)*($$p-$$a);$$w=($$d-$$a)*($$q-$$b)-($$e-$$b)*($$p-$$a);

if ($$u != 0) { s=v/u; t=w/u; intercept=(u*a + v*b + w*c)/$$u;

print OUT "$$k\t$$l\t$$m\t"; print OUT "$$s\t$$t\t$$intercept\n";

} else {
print OUT "$$k\t$$l\t\$m\tERROR (u=0)\n";
}
}
close(OUT);