More pseudo random number generators for calculators

[Namir]

Hi All,

I posted a new study about "more random number generators for calculators." The study takes the efficient algorithm used by HP in the HP-67 Stat Pac and attempts to tweak it to find better variants. The study looks at 55 algorithms and ranks them using Matlab code. I use Matlab because it can handle very large arrays which I needed to get better results.

To download the pdf article click here.

To download the Matlab code click here.

Enjoy!

Namir

[Erik-M]

Hello Namir,

For me this is a completely different way to look at PRNG's.
I am active in telecom and know the shift register based maximum length generators using a polynome.
Thanks a lot for this study, I have something to read and play with this weekend.

Erik

[ttw]

A short search on the HP50g gives the following parameters that should work there.

X(j+1)=X(j)*227716313950 mod 549755813507
X(0) not 0

This should generate all integers less than the modulus except zero. The numbers should be small enough not to screw up the division on the HP50g.

(All this assumes that I copied the numbers from the display into the computer correctly.)

[Albert Chan]

Hi, Namir

Some algorithms you tested might get stuck, even with non-zero fractional part guess

Example, say with 12 digits calculator:

r = 0.999 999 999
r = frac(997/r)     → r = frac(997.000000 997000000 997000000 ...) = .000 000 997
r = frac(997/r)     → r = 0
r = frac(997/r)     → r = NaN

To avoid rounding errrors like above, multiply based versions is safer.
With 15 digits internal precision, and 3 digits factor k, r = frac(k*r) is guaranteed exact.

This may be why baseline use r = frac(997*r)

[ttw]

Good points about getting stuck. The theory may show cyclic (or at least unending wandering) but rounding can cause fixed points to appear, usually 0.0.

For my big integer version above. I set the modulus with MODSTO just use MULTMOD. Then I can divide by the modulus to scale (I kept the modulus so it should be representable.) No problem with low order bits. (All assuming I copied things correctly.)

[Namir]

(07-09-2020 11:19 AM)Albert Chan Wrote:  Hi, Namir

Some algorithms you tested might get stuck, even with non-zero fractional part guess

Example, say with 12 digits calculator:

r = 0.999 999 999
r = frac(997/r)     → r = frac(997.000000 997000000 997000000 ...) = .000 000 997
r = frac(997/r)     → r = 0
r = frac(997/r)     → r = NaN

To avoid rounding errrors like above, multiply based versions is safer.
With 15 digits internal precision, and 3 digits factor k, r = frac(k*r) is guaranteed exact.

This may be why baseline use r = frac(997*r)

You certainly have a point. However, in each "run" I would generate a billion random numbers and I never had a runtime error. So the likelihood of getting a random number very close to 1 is small. Of course dealing with floating point number depends on many aspects of the software being used.

[Gene]

(07-09-2020 07:35 PM)Namir Wrote:  You certainly have a point. However, in each "run" I would generate a billion random numbers and I never had a runtime error. So the likelihood of getting a random number very close to 1 is small. Of course dealing with floating point number depends on many aspects of the software being used.

Gene: Which is CERTAINLY more than any HP model other than the PRIME would ever need. :-)

[Namir]

(07-09-2020 09:15 PM)Gene Wrote:  

(07-09-2020 07:35 PM)Namir Wrote:  You certainly have a point. However, in each "run" I would generate a billion random numbers and I never had a runtime error. So the likelihood of getting a random number very close to 1 is small. Of course dealing with floating point number depends on many aspects of the software being used.

Gene: Which is CERTAINLY more than any HP model other than the PRIME would ever need. :-)

Absolutely right!

[Namir]

(07-09-2020 03:45 AM)ttw Wrote:  A short search on the HP50g gives the following parameters that should work there.

X(j+1)=X(j)*227716313950 mod 549755813507
X(0) not 0

This should generate all integers less than the modulus except zero. The numbers should be small enough not to screw up the division on the HP50g.

(All this assumes that I copied the numbers from the display into the computer correctly.)

I ran your algorithm in my extended calculations. It came out 9th among 56 algorithm (counting the algorithm itself). So it's an OK algorithm and works with calculators that can store the two big numbers and handle well the multiplication of large integers.

Cheers,

Namir

[ttw]

The generator I posted was the result of a 5-minute search. I'm going to find one for some computations but getting "better" (however measured) takes more time.

I use long LCGs for Monte Carlo (whenever my quasi-Monte Carlo stuff doesn't apply) as the precision does matter. I have a short generator on the calculator (or I did, it's easy to redesign) with three shift registers lengths 15,16,17, combined with a 16-bit GCD. The total cycle is (2^15-1)*(2^16-1)*(2^17-1)*2^16 (the cycle lengths are prime) giving just under 2^64 as a cycle. However, the precision is only 15 bits. A true 64-bit RNG would have a precision of 64 bits.

A problem arises in simulating exponential or Gaussian distributions by the Box-Muller method. With fifteen bits, the greatest exponential number is ln(2^16) or about 11 whereas with 64 bits one gets ln(2^64) or 44. With the Gaussian, the largest number generated is about Sqrt(2*Ln(2^16)) or a sigma of 4.7. Not even 6-sigma quality. With 64 bits one gets Sqrt(2*Ln(2^64)) or 9.4. The smaller precision does not explore rare events. In real life, one should use a rare-event generator.

[Namir]

(07-10-2020 11:25 PM)Namir Wrote:  

(07-09-2020 03:45 AM)ttw Wrote:  A short search on the HP50g gives the following parameters that should work there.

X(j+1)=X(j)*227716313950 mod 549755813507
X(0) not 0

This should generate all integers less than the modulus except zero. The numbers should be small enough not to screw up the division on the HP50g.

(All this assumes that I copied the numbers from the display into the computer correctly.)

I ran your algorithm in my extended calculations. It came out 9th among 56 algorithm (counting the algorithm itself). So it's an OK algorithm and works with calculators that can store the two big numbers and handle well the multiplication of large integers.

Cheers,

Namir

I also tested:

X(j+1)=(X(j)*227716313950 + 22771631395) mod 549755813507

which gave satifactory (better than r = frac997*r) results. The above algorithm handles zero values for X(j).