Recursive random number generator based on irrational numbers

Question

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.

Answer 1

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

Answer 2

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);