Information on calculator Random Number Generators from PPC Journal articles

[Dieter]

(11-15-2016 09:32 PM)Paul Dale Wrote:  You didn't include the \( frac ( 2^{20} x_n ) \) part. The leading bits don't pass the randomness tests but going down the mantissa a bit does.

Sorry, it was not clear to me that R_n is returned and not x_n. Adding the 20-bit shift improves the results very much, if (!) the algorithm runs in a regular binary double precision environment with 53 bit mantissa. On the other hand the method may not be the first choice on a standard 10-digit BCD calculator: here only three digits may remain in the fractional part after the multiplication by 2^20. I tried shifting by 13 bits, i.e. x_n*8192, which yields 2..3 more digits, but here other problems arise.

Dieter

[Paul Dale]

For a BCD device, I'd multiply by a power of ten instead.


Pauli

[Paul Dale]

(11-15-2016 05:12 PM)Namir Wrote:  I checked the method with Excel VBA and realized that the random calculation noise (generated when you take the fractional part of large numbers) was particular to Matlab!

Might the noise be correct and Excel is doing naughty things in the background?
See e.g. pages 7 & 8 of https://people.eecs.berkeley.edu/~wkahan/ARITH_17.pdf


Pauli

[Dieter]

(11-15-2016 10:40 PM)Paul Dale Wrote:  For a BCD device, I'd multiply by a power of ten instead.

OK, but 2^20 = 1.048.576 or 1.000.000 does not make much of a difference. ;-)

Dieter

[Paul Dale]

It might be possible to reduce the multiplier to \( 10^4 \) or \( 10^5 \)without losing the apparent entropy. The choice of \( 2^{20} \) was fairly arbitrary.

Pauli

[Dieter]

(11-15-2016 11:51 PM)Paul Dale Wrote:  It might be possible to reduce the multiplier to \( 10^4 \) or \( 10^5 \)without losing the apparent entropy. The choice of \( 2^{20} \) was fairly arbitrary.

I tried a multiplier of 1000 which leaves 5 or 6 digits on a 10-digit calculator. Which is not too much since some game programs use the first four digits to get two 2-digit random numbers. #-)

After a few runs with 10 batches of 1000 numbers each I discovered a similar phenomenon like the one you can see in the table at post #39. This time there were several 1000-number-runs where not a single random number appeared between 0,7 and 0,8. The pattern was like this, leading to Chi² = 202,9. While the first three 10.000-number-runs looked fine, all following ones had this pattern, maybe here and there ±1 in the one or other bin, but always the 8th one empty:

Code:

0,0...    117
0,1...    103
0,2...     78
0,3...    143
0,4...    182
0,5...     91
0,6...    117
0,7...      0  (!!!)
0,8...     78
0,9...     91

Very strange. I wonder if it's the RNG, Excel or my coding.
The random numbers were calculated this way, starting with a seed of x0 = 0.546703671.

Code:

    z = Exp(Round(x + x + 1, 9))
    x = z - Int(z)
    If z < 9.9999999995 Then x = Round(x, 9) Else x = Round(x, 8)
    rn = Round(1000 * x - Int(1000 * x), 6)

Any idea what would cause strange results like these?

Dieter

[Dieter]

(11-16-2016 01:53 PM)Dieter Wrote:  Any idea what would cause strange results like these?

I did another test. This time starting with x0 = 0,45642414. The implementation mentioned above generates the following sequence:

0,350905
0,42875
0,62931
0,93117
0,02565
0,42874
0,643333
...

A sample of 1000 consecutive random numbers showed a period length of 77. The 78th number is 0,350905 again, and the above sequence starts over. In fact the 78th x-value is 0,45642414 again, the original seed. Thus in 1000 numbers every value occured 13 times. OK, except one since 13x77=1001. This exeption is the seed itself, i.e. Rn=0, 42414. And in fact there is not one single Rn in [0,7; 0,8[. In the sorted list of these 1000 random numbers 0,69777 is directly followed by 0,84506.

I have also checked this with an HP41 implementation:

Code:

RCL 00
ST+ X
1
+
e^x
FRC
STO 00
1E3
*
FRC

Starting with 0,45642414 in R00 the program yielded exactly the same results, and after 77 numbers the sequence restarted, with 0,45642414 in R00 again.

For the record, here is the complete list of the first 80 Rn:

Code:

 1.  0,350905
 2.  0,42875
 3.  0,62931
 4.  0,93117
 5.  0,02565
 6.  0,42874
 7.  0,643333
 8.  0,85718
 9.  0,10318
10.  0,65815
11.  0,557927
12.  0,47512
13.  0,092
14.  0,95182
15.  0,41278
16.  0,07307
17.  0,505942
18.  0,303103
19.  0,382599
20.  0,094476
21.  0,436495
22.  0,03523
23.  0,02344
24.  0,285095
25.  0,85653
26.  0,675471
27.  0,37482
28.  0,54911
29.  0,95747
30.  0,57554
31.  0,214103
32.  0,849191
33.  0,309908
34.  0,46272
35.  0,4397
36.  0,12143
37.  0,428283
38.  0,038915
39.  0,116376
40.  0,383202
41.  0,29357
42.  0,25395
43.  0,516251
44.  0,94027
45.  0,63507
46.  0,687494
47.  0,342444
48.  0,06445
49.  0,875943
50.  0,521511
51.  0,24125
52.  0,17458
53.  0,424301
54.  0,141801
55.  0,332297
56.  0,213162
57.  0,144471
58.  0,10338
59.  0,070719
60.  0,872262
61.  0,370813
62.  0,99494
63.  0,482137
64.  0,69777
65.  0,11986
66.  0,960758
67.  0,3351
68.  0,43555
69.  0,69399
70.  0,43742
71.  0,391572
72.  0,84506
73.  0,596742
74.  0,41389
75.  0,913563
76.  0,60827
77.  0,42414
78.  0,350905   ...and back to start
79.  0,42875
80.  0,62931
     ...

Hm, hmm, hmmm...

Dieter

[Paul Dale]

Working in higher precision, the \( e^{2x+1} \) generator seems to work okay with a scale factor of \( 10^5 \) but no so well with \( 10^4 \). However, I've not tried limiting calculations to 10 digits as you are doing.

Pauli

[Namir]

In the previous messages in this thread, the following similar PRNGs were mentioned:

r=frac((pi+r)^3)
r=frac((pi+r)^5)

Their general form is:

r=frac((a+r)^b)

I chose a = pi + cos(k * c2 / pi) and b = 4 + sin(k * c1 / pi)

Where c1=127, c2=113, and k is a counter that is incremented every time we invoke the PRNG function.

I tested 5 runs of 10,000 random numbers. The maximum absolute auto-correlation was less than 0.03. The Chi-Square was in (4, 10.3) range. The mean and sdev of all 50,000 random numbers were 0.500016257 and 0.288179175, respectively.

The above PRNG should be adequate for calculation apps, simulations, and games.

Enjoy!

Namir

PS: You can also use a = pi + sin(k * c2 / pi) and b = 4 +cos(k * c1 / pi) and get similar results.