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