HP41C random number generator help

[ttw]

Here's the pseudo code for a uniform random number generator. It's not the smallest that I've got but it is extremely fast. There's some theory showing that one gets a uniform distribution in theory.

There is a constant: c=fractional part of fractional irrational. If 4*k+1 is prime, Frac(Sqrt(4*k+1)) for any k is good as is Frac((Sqrt(4k+1)-1)/2). Frac(Sqrt(2)) is good too.
Sqrt(2)
(Sqrt(5)-1)/2
(Sqrt(13)-1)/2
etc.

To avoid some rounding error, one can also use the auxiliary constant c2= Frac(2*c).

Two sequences are generated, the "random" sequence U and an auxiliary sequence A using the following rules.

U(0)=0
A(0)=0

Then compute:
U(n+1)=Frac( u(n)+A(n)+c )
A(n+1)=Frac( A(n)+c2 )
0r A(n+1)=Frac(A(n)+2*c)
This is a finite difference sequence to computing Frac(n^2*c).

There is a rational alternative using c as BigInteger1/BigInteger2 which also works with no rounding error but needs long integer arithmetic.

There are two easy extensions to the sequence. One can generate two independent streams by generating separate sequence with different c's. For example Frac(Sqrt(2)) and Frac(!Sqrt(5)-1)/2). One can also start the sequence with U=Frac(X) where X is a number. I've used X as some type of combination of clock readings. I can run a simulation using different values of c for independent things within the simulation and using different starting values for U to get different samples.

If one computes (slightly different in rounding error)
U(0)=0
A(0)=0
U(n+1)=Frac(U(n)+2*A(n)+c)
A(n+1)=A(n)+c
Then sequence A(n) is a low-discrepancy sequence; uniformly distributed by approaching the uniform distribution much faster than a random sequence. The related sequence U(n) behaves more like a random sequence (when sorted, the differences between members of the sequence have a Poisson distribution.)