New PRNG for calculators

[Namir]

I was just tinkering with an HP-21 calculator and I tried the following (I had to do the INT part by hand):

r1 = Frac(exp(exp(exp(r0))))

Where r0 is the old random number, and r1 is the new random number. I then wrote a Matlab program to generate a million random number. Here is a sample set of results:

min = 1.646694727241993e-07
max = 0.999999882210702
mean = 0.500182006746859
sdev = 0.288245734345533

For the auto correlation between 1000 elements:

min = -0.003854280151534
max = 0.003319062287133
mean = 1.567008799132003e-05
sdev = 0.001003357884267

The method looks good for random numbers for calculators.

Namir

[KeithB]

Don't you mean Frac() to get a number between 0 and 1?

[Namir]

(06-06-2024 08:30 PM)KeithB Wrote:  Don't you mean Frac() to get a number between 0 and 1?

Yes and I corrected it! Thanks!

[ttw]

Does the loss of significance (Exp(Exp(Exp(1))) is 3814279.1...) give any trouble? On an HP50g, the smallest result before the MOD 1 reduction is 15.1542622415 and the largest is 3814279.10484. The next inputs would differ in the last 4 digits. This may have no significance as the output wraps around the unit interval 15 or more times.

Interval transformations are not preserved.
Transformation of (.10001,.10002) yields: (.486654,.48734).
Transformation of (.20001, .20002) yields: (.72485, .72609).
Transformation of (.65001, .65002) yields ( .57199, .68771).
Transformation of (.69001, .69002) yields (.20483, .43117).
Transformation of(.69003, .69004) yields (.65744, .88376).

The successors of numbers falling to a small interval on the 0 side of the unit interval have a very different distribution from those falling to a small interval on the 1 side.

[Namir]

(06-07-2024 01:35 AM)ttw Wrote:  Does the loss of significance (Exp(Exp(Exp(1))) is 3814279.1...) give any trouble? On an HP50g, the smallest result before the MOD 1 reduction is 15.1542622415 and the largest is 3814279.10484. The next inputs would differ in the last 4 digits. This may have no significance as the output wraps around the unit interval 15 or more times.

Interval transformations are not preserved.
Transformation of (.10001,.10002) yields: (.486654,.48734).
Transformation of (.20001, .20002) yields: (.72485, .72609).
Transformation of (.65001, .65002) yields ( .57199, .68771).
Transformation of (.69001, .69002) yields (.20483, .43117).
Transformation of(.69003, .69004) yields (.65744, .88376).

The successors of numbers falling to a small interval on the 0 side of the unit interval have a very different distribution from those falling to a small interval on the 1 side.

I guess MATLAB uses more significant figures than typical calculators. How many random numbers did you generate?

Namir

[Namir]

(06-07-2024 02:05 AM)Namir Wrote:  [quote='ttw' pid='187994' dateline='1717724102']
Does the loss of significance (Exp(Exp(Exp(1))) is 3814279.1...) give any trouble?

I feed the random-number generating functions with numbers < 1.

Namir

[Matt Agajanian]

Hi. Excellent challenge. I just created an HP-25 version.

I also prefaced the code with a routine that scales down seeds to values 0<=seed<=1

Here it is:

Code:


To use:

Enter seed
Press CLEAR: PRGM
Press R/S repeatedly to generate successive results.

01 - 31          ENTER
02 - 15 03       ABS
03 - 01          1
04 - 51          +
05 - 14 01       INT
06 - 14 08       LOG
07 - 01          1
08 - 51          +
09 - 14 01       INT
10 - 15 08       10ˣ
11 - 71          ÷
12 - 23 01       STO 1
13 - 15 07       eˣ
14 - 15 07       eˣ
15 - 15 07       eˣ
16 - 15 01       FRAC
17 - 74          R/S
18 - 13 12       GTO 12

[ttw]

(06-07-2024 02:05 AM)Namir Wrote:  I guess MATLAB uses more significant figures than typical calculators. How many random numbers did you generate?

Namir

None. I just looked at small intervals. The image of sufficiently small intervals with the same length transforms to another image that doesn't exceed 1.0 in length, but it does differ with the position of the interval. Knowing the interval allows one to guess the output of the next number (at least its range) with pretty good probability. For those that wrap around quite a bit, (near 1.0), one knows little about the next number; .90001's successor is .555669 (to HP50g precision)' .90002's successor is .276095. These results come from the pre-reduction interval 120626.555669 and 120661.27695. Neither of these results helps in figuring out .900015's successor (which is .914477). The successor of .10001 is 20.4866596716 and of .10002 is 20.487343122; the successor of .10015 is 20.487001537. This lies between the successors of .10001 and .10002 both before and after the MOD 1 mapping.

Exp(x) expands differently for each x. For this reason, mappings within finite fields (like shift registers) or rings (like linear congruential generators) are often used. Points within a small interval, no matter where in the range of the generator, will get spread out approximately uniformly over the whole space.

Without doing any analysis, I think you could make the exp(exp(exp(x))) work by mapping the output to binary ,then XORing with the output of a linear congruential generator.

[Namir]

Here are more results for TEN MILLION random numbers and TEN THOUSAND lags. These results are more than enough for using random numbers in calculator simulations. I also rounded the random numbers to 12 decimals to give results similar to popular programmable calculators.

Random numbers rounded to 12 decimals
With random init seed
Number of points = 1.000000e+07
min(x) = 0.0000018142
max(x) = 0.9999992789
mean(x) = 0.5002791084
std(x) = 0.2878110140
Max lags = 10000
min(acf(2:end) = -9.59039166e-03
max(acf(2:end)) = 9.51402836e-03
mean(acf(2:end)) = 2.32021982e-05
std(acf(2:end)) = 2.60088239e-03

In the above results the mean is close to 0.5 and the standard deviation is close to 0.288. The minimum and maximum random numbers cover a good part of the range [0, 1].

The auto-correlation values are taken of 10,000 lags and show a weak correlation between the random numbers with values in the range [-0.0095, 0.0095].

These two sets of results give a preliminary assessment that the algorithm is good and needs no further modifications.

I rand two batches for the initial seeds of 0 and 1. The results are similar to the above. In practical use, one would provide an initial seed other than 0 or 1, and with many decimal places.

With init seed of 0
Number of points = 1.000000e+07
min(x) = 0.0000065295
max(x) = 0.9999998871
mean(x) = 0.4995039774
std(x) = 0.2882581865
Max lags = 10000
min(acf(2:end) = -4.78069967e-03
max(acf(2:end)) = 4.49636143e-03
mean(acf(2:end)) = 1.24859064e-05
std(acf(2:end)) = 1.17159700e-03

With init seed of 1
Number of points = 1.000000e+07
min(x) = 0.0000000721
max(x) = 0.9999992789
mean(x) = 0.5003003014
std(x) = 0.2878265888
Max lags = 10000
min(acf(2:end) = -9.23315691e-03
max(acf(2:end)) = 9.30401849e-03
mean(acf(2:end)) = 2.19537931e-05
std(acf(2:end)) = 2.53885074e-03

[Albert Chan]

Hi, Namir

I just tested e^e^e random generator ... so far, numbers look good.
Chi square numbers is a bit weak though.

ent.lua generate 1e6 random byte, to be feed to ent.exe for analysis.

Code:

-- usage: run ent.lua [seed] | ent.exe
local seed  = tonumber((...) or os.time())

math.randomseed(seed)
local exp, floor = math.exp, math.floor
local s = math.random()
local r = function(n) s = exp(exp(exp(s)))%1; return 1+floor(n*s) end

io.stderr:write('seed = ', seed, '\n')
local write, c = io.write, string.char

for i=1,1e5 do
    write(
        c(r(256)-1), c(r(256)-1), c(r(256)-1), c(r(256)-1), c(r(256)-1),
        c(r(256)-1), c(r(256)-1), c(r(256)-1), c(r(256)-1), c(r(256)-1)
    )
end

Note: seed = 0 --> s = 0.8438413544573831

c:\random> lua ent.lua 0 | ent
seed = 0
Entropy = 7.999764 bits per byte.

Optimum compression would reduce the size
of this 1000000 byte file by 0 percent.

Chi square distribution for 1000000 samples is 327.60, and randomly
would exceed this value 0.14 percent of the times.

Arithmetic mean value of data bytes is 127.3877 (127.5 = random).
Monte Carlo value for Pi is 3.144708579 (error 0.10 percent).
Serial correlation coefficient is 0.002659 (totally uncorrelated = 0.0).

---

For comparison, this is what happen if I use default, r = math.random
Random generator used https://prng.di.unimi.it/xoroshiro128plus.c

c:\random> lua ent.lua 0 | ent
seed = 0
Entropy = 7.999804 bits per byte.

Optimum compression would reduce the size
of this 1000000 byte file by 0 percent.

Chi square distribution for 1000000 samples is 271.73, and randomly
would exceed this value 22.52 percent of the times.

Arithmetic mean value of data bytes is 127.6241 (127.5 = random).
Monte Carlo value for Pi is 3.145572582 (error 0.13 percent).
Serial correlation coefficient is 0.000624 (totally uncorrelated = 0.0).

[ttw]

The function as written Frac(exp(exp(exp(x))) is uniformly distributed mod 1. So the OP's pRNG should work, except for rounding and discretization errors (All pRNGs have this problem.).

Frac(Exp(x)) is also uniformly distributed, Frac(Ln(x)) isn't. I found a paper proving this result and giving criteria for the uniform distribution of some types of functions.

https://dwc.knaw.nl/DL/publications/PU00018885.pdf

Daniel Shanks: "ln( ln( ln ( x ) ) ) approaches infinity with great dignity."

(It's a comment about the inverse function. I've only seen ln(ln(ln(ln(x)))) in practice.)

[Albert Chan]

(06-08-2024 02:51 AM)ttw Wrote:  The function as written Frac(exp(exp(exp(x))) is uniformly distributed

How to we know this?

Quote:Frac(Exp(x)) is also uniformly distributed

I don't think it is.

< local r = function(n) s = exp(exp(exp(s)))%1; return 1+floor(n*s) end
> local r = function(n) s = exp(s)%1; return 1+floor(n*s) end

c:\random> lua ent2.lua 0 | ent
seed = 0
Entropy = 2.686859 bits per byte.

Optimum compression would reduce the size
of this 1000000 byte file by 66 percent.

Chi square distribution for 1000000 samples is 126198648.37, and randomly
would exceed this value less than 0.01 percent of the times.

Arithmetic mean value of data bytes is 10.5856 (127.5 = random).
Monte Carlo value for Pi is 3.971823887 (error 26.43 percent).
Serial correlation coefficient is 0.871967 (totally uncorrelated = 0.0).

[ttw]

After more reading and a night's sleep, it's the entire line from -inf to +inf for which frac(x) is uniform. From (0-1) it may not be. Also as an iteration, a single x0's orbit may not be uniform. My pRNGs tend to be seedless for that reason. (I have stolon the idea from St Augustin grass.)

[Albert Chan]

(06-08-2024 03:44 PM)ttw Wrote:  After more reading and a night's sleep, it's the entire line from -inf to +inf for which frac(x) is uniform.

I would have guessed FP(exp(x)) would skew to tiny values, not uniform.

Negative side, x<ln(0.1), FP(exp(x)) = exp(x) < 0.1
Positive side, numerically FP(exp(x)) = 0, if x goes big. (assume it did not overflow)

The other two exp() is needed to keep x in reasonable domain.

e^e^0 .. e^e^1 = e .. e^e ≈ 2.71828 .. 15.1543

[SlideRule]

Stumbled on Example 4: HP-25 Pocket Calculator, Computational Statistics, 25^th Conference, page 131

"The pseudorandom numbers of the HP-25 pocket calculator are generated
by the algorithm
                               u_i ≡ ( u_i-1+3)⁵    (mod 1)

With a seed u₀, 0 < u₀ < 1, the floating-point arithmetic creates sequences
which behave as demonstrated in Example 1. After a more or less ran-
dom behaviour they enter a short subsequence repeating itself. The HP-25
calculator may enter such a cycle of only 29 numbers.

This effect can be attributed to the rounding errors of the floating-point
arithmetic. The generators of the BASIC interpreters on the Commodore
and Apple microcomputers (CBM PET 2001 Series and Apple II europlus )
are further examples for generators having this defect. They are described
and analysed in Afflerbach (1985) and Ripley ( 1987, p.18 ). In Afferbach
(1985) it is shown that they may reach short cycles of period 202 and 703,
respectively. Examples 3 and 4 elucidate why congruential generators de-
fined by integer calculations ( and division by the modulus ) should be pre-
ferred to congruential generators mod 1. More about the implementation of
pseudorandom number generators can be found in Gentle (1990) and Ripley
( 1990 )."

BEST!9SlideRule
argh - yet another senior moment!

[KeithB]

A reference from 1900!

How can u0 be less than u0?

[EdS2]

I'm guessing it's Ripley from 1990
Thoughts on pseudorandom number generators B.D. RIPLEY
(pdf here)

[Namir]

I did a preliminray test for the mean, standard deviation and autocorrelation values for random numbers generated with the above algorithm and the results are good.

I am thinking about doing a study for the general algorithm r = frac((r+ n)^m) for a parameters matrix for various integer values of n and m (and perhaps n+pi and m). The study would use my penalty factor stat that I have presented in various HHC meetings. The penalty factor calculations use a number of tests (beyond just the three mentioned above) to investigate the quality of randomness generated. Of course the results wouldbe presented in an upcoming HHC meeting.

Namir

[KeithB]

Wow, you type all these by hand?

[Namir]

(08-08-2024 05:48 PM)SlideRule Wrote:  Stumbled on Example 4: HP-25 Pocket Calculator, Computational Statistics, 25^th Conference, page 131

"The pseudorandom numbers of the HP-25 pocket calculator are generated
by the algorithm
                               u_i ≡ ( u_i-1+3)⁵    (mod 1)

With a seed u₀, 0 < u₀ < 1, the floating-point arithmetic creates sequences
which behave as demonstrated in Example 1. After a more or less ran-
dom behaviour they enter a short subsequence repeating itself. The HP-25
calculator may enter such a cycle of only 29 numbers.

This effect can be attributed to the rounding errors of the floating-point
arithmetic. The generators of the BASIC interpreters on the Commodore
and Apple microcomputers (CBM PET 2001 Series and Apple II europlus )
are further examples for generators having this defect. They are described
and analysed in Afflerbach (1985) and Ripley ( 1987, p.18 ). In Afferbach
(1985) it is shown that they may reach short cycles of period 202 and 703,
respectively. Examples 3 and 4 elucidate why congruential generators de-
fined by integer calculations ( and division by the modulus ) should be pre-
ferred to congruential generators mod 1. More about the implementation of
pseudorandom number generators can be found in Gentle (1990) and Ripley
( 1990 )."

BEST!9SlideRule
argh - yet another senior moment!

I found my old Matlab code that I had used to test PRNGs for calculators. I adapted the code for testing r=frac((3+r)^5) (as the first of many variants). I also made Matlab round the random numbers to 10 decimals to make then comparable to vintage HP calculators. The results were impressive and better than my previous favorite r=frac(r*997).

So now I will study the variants of that algorithms by changing the integer shift and power values.

Namir