I recently browsed through someone else's code and found a section where a pseudo random number generator is implemented. I know that random number generation is not an easy task, some even regard as art, so my question is: Is this a known generator? Can I trust it? What are good/bad seeds?
The generator
The first thing I noticed was that the algorithm operates on floating point numbers instead on integers. It takes three doubles
as a seed, which I will call $(s_1, s_2, s_3)$ in the following. These are converted to the internal states $(a_n, b_n, c_n)$ by this map:
$$a_0 = \sqrt{s_1} - 1 \\ b_0 = \sqrt{s_2} - 1 \\ c_0 = \sqrt{s_3} - 2$$
Then the algorithm proceeds as follows:
$$\begin{pmatrix} a_{n+1} \\ b_{n+1} \\ c_{n+1} \end{pmatrix} = \operatorname{mod} \left[ \begin{pmatrix} 1 & 1 & 1 \\ 1 & 2 & 2 \\ 2 & 3 & 4 \end{pmatrix} \begin{pmatrix} a_n \\ b_n \\ c_n \end{pmatrix} , 1 \right] $$
The first random number then is $a_1$ and $a_2, a_3, \ldots$ are the next ones. $\operatorname{mod}(x,1)$ means that we take the decimal places, e.g. $\operatorname{mod}(4.345,1) = 0.345$. This makes sure that the generated numbers are smaller than one.
If you take a look at the seeding procedure again, then you can see that seeding with $(1,1,4)$ is obviously a bad choice - this generates only zeros. Lets assume all seeds are larger than this tuple, then it apparently generates numbers in [0,1).
What I have found out so far
In order to classify this generator I think it is helpful to first rewrite it into an algorithm that operates on integers instead of floating point numbers: I can assume that $a_0 = a_0' / q_1$, $b_0 = b_0' / q_2$, and $c_0 = c_0' / q_3$ where the nominators and denominators are integers (because double
has finite precision). Then the first line reads
$$a_1 = \operatorname{mod} \left( \frac{q_2 q_3 \cdot a_0' + q_1 q_3 \cdot b_0' + q_1 q_2 \cdot c_0'}{m}, 1 \right)$$
where $m = q_1 q_2 q_3$ is a common denominator. Pulling out the factor $m$, all three lines can be rewritten to
$$\begin{pmatrix} a_1 m \\ b_1 m \\ c_1 m \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 2 & 2 \\ 2 & 3 & 4 \end{pmatrix} \begin{pmatrix} q_2 q_3 \cdot a_0' \\ q_1 q_3 \cdot b_0' \\ q_1 q_2 \cdot c_0' \end{pmatrix} \mod m$$
with $x \mod y$ having the usual meaning.
This looks like some kind of coupled system of linear congruential generators (LCG). Maybe one can use the chinese remainder theorem and reduce the sum of the three LCGs to one, similar to the answer in this question: Flaw or no flaw in MS Excel's RNG? I am not sure how to deal with the matrix, maybe this introduces some kind of lag?