HP41C random number generator help

[Trond]

Hi folks!
I just got my HP 41CV and I am a bit rusty in HP programming in general (never mind this being an older one than I ever had before)

I just looked at the different random number generator programs and I am puzzled. How is this even supposed to work?
Take this simple one:

01 LBL "RNG3"
02 LN
03 E2
04 *
05 1
06 MOD
07 END

So this program takes the natural logarithm of whatever is in you x stack register, and I am a bit uncertain what the E2 does here? And should I just type it as it stands, or do they mean e^2 (but that would be typed 2 E^X)? I tried one of the other programs too, but it did not return numbers between 0 and 1.

[Paul Dale]

E2 is EEX 2 which puts 100 on the stack.

The formula boils down to: x_n+1 = FRAC(100 * LN(x_n)).
I suspect this isn't a great generator but it is very short.

Pauli

[Trond]

Thanks! Man, I should have gotten that.

Anyway, the 1 MOD, just removes the whole integers doesn't it? Would it not be easier to type FRC?

[Paul Dale]

How does FRC deal with negative numbers? LN(x) for 0 < x < 1 is negative.


Pauli

[Sylvain Cote]

HP used this random number generator in lots of places over the years.

FRAC((x * 9821) + 0.211327)

From the GAMES pac:

LBL "RNDM"
RCL 00
9821
*
.211327
+
FRC
STO 00
RTN

Sylvain

[Trond]

Right.

Hm, still does not work. Let me know if it works for anybody else, or if you would recommend another random number generator.

EDIT: OK thanks Sylvain. I will try that one.

[Trond]

(12-05-2017 06:08 AM)Sylvain Cote Wrote:  HP used this random number generator in lots of places over the years.

FRAC((x * 9821) + 0.211327)

From the GAMES pac:

LBL "RNDM"
RCL 00
9821
*
.211327
+
FRC
STO 00
RTN

Sylvain

Looks like it works. Thanks!

[Dieter]

(12-05-2017 05:45 AM)Trond Wrote:  So this program takes the natural logarithm of whatever is in you x stack register,

It it supposed to take the ln of the previous random number. This is how pseudo random number generators work: they calculate the next number from the last one.

(12-05-2017 05:45 AM)Trond Wrote:  and I am a bit uncertain what the E2 does here? And should I just type it as it stands, or do they mean e^2 (but that would be typed 2 E^X)?

As already noted, this is a shortcut for 1E2 or simply 100. With a bit of synthetic programming the leading "1" can be removed, saving one byte. If you happen to see a line with a single "E" that's the same thing, just the nerd way of writing a simple "1". ;-)

(12-05-2017 05:45 AM)Trond Wrote:  I tried one of the other programs too, but it did not return numbers between 0 and 1.

?!? What "other programs"?

The last step of most PRNGs is a FRC command. This removes the integer part and leaves only the decimals, so you automatically get something between 0 and (less than) 1.

Random number generators have been discussed here several times, just do a search and you'll find some suggestions. You may also take a look at PPC Journal V4N8 and V8N3.

The following examples assume the current random number (the seed) in R00.

This one is fast and good enough for games and similar purposes. I use it in most of my programs.

Code:

RCL 00
Pi
+
x²
FRC
STO 00

This one is supposed to have a somewhat better number distribution:

Code:

RCL 00
Pi
+
ENTER
x²
*
FRC
STO 00

Edit: replaced x⁵ with x³.

The following generator was used in early HP programs, but best results are only obtained if the seed has certain properties. And it obviously does not work with a seed of zero.

Code:

RCL 00
997
*
FRC
STO 00

And finally this one, already mentioned by others, has good statistical properties and produces six-digit random numbers. It runs faster if you store the two constants in data registers.

Code:

RCL 00
9821
*
,211327
+
FRC
STO 00

Dieter

[Dieter]

(12-05-2017 06:06 AM)Paul Dale Wrote:  How does FRC deal with negative numbers? LN(x) for 0 < x < 1 is negative.

FRC of a negative number is negative as well, while 1 MOD is positive as it adds 1:
FRC(–3,14)=–0,14  while  –3,14 MOD 1 returns 0,86. At least on the HP41.

Maybe x := FRC(–100*ln x) should be used here. Or even better x := FRC(ABS(100*ln x)).

Dieter

[Namir]

A good and simple algorithm (used by HP since the HP-65 Stat Pacs) is:

r = frac(997 * r)

Code:

LBL "RAND"
RCL 00
997
*
FRC
STO 00
RTN

In my HHC2017 presentation about PRNGs, I pointed out that the above legacy algorithm is recommended for its speed, simplicity, and relatively good results for calculators.

Namir

[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.)

[Trond]

(12-05-2017 08:23 AM)Dieter Wrote:  ?!? What "other programs"?

http://www.hpmuseum.org/software/41/41ranjm.htm

(12-05-2017 12:30 PM)Namir Wrote:  A good and simple algorithm (used by HP since the HP-65 Stat Pacs) is:

r = frac(997 * r)

Code:

LBL "RAND"
RCL 00
997
*
FRC
STO 00
RTN

In my HHC2017 presentation about PRNGs, I pointed out that the above legacy algorithm is recommended for its speed, simplicity, and relatively good results for calculators.

Namir

Wow, that is very simple. But it seems to assume that you start up with a number (seed) that is already non-0 and has a fraction part. Some of the other programs seem to assume similar things. Sorry if I am missing some assumptions that everyone is making here.

EDIT: ah never mind. I see Dieter already covered the same.

I see that Dieter already suggested that we can start off with pi. This is a technique I remember someone mentioning years ago, another time, another place. Maybe I could start Namir's program with adding pi?

[Trond]

Just out of curiosity, does anyone know which one is used in the RAN function of HP 42S ?

[Thomas Okken]

(12-05-2017 03:53 PM)Trond Wrote:  Just out of curiosity, does anyone know which one is used in the RAN function of HP 42S ?

The RPL random number generator is described in full detail in this thread: https://groups.google.com/d/msg/comp.sys...JTaWW0OeUJ

The only difference between it and the HP-42S RAN function is the initial seed after hard reset: in the RPL version, it is 999500333083533, while in the 42S version, it is 2787.

N.B. Free42 >= 2.0.7 uses this RNG, so if you're looking for working code, you can find it there.

[Trond]

Thanks again folks!

Wow, this calculator is really kickstarting my brain into action again. I am loving this.

I have programmed some good user keys; one for random numbers, with the algorithm Sylvain suggested. And one that is a "D6", simulating a normal die with six sides (yes I am a gamer, at least sometimes) I also assigned two keys for SDEV and MEAN, just to easily check on how the random generators are working.

[Dieter]

(12-05-2017 03:24 PM)Trond Wrote:  Wow, that is very simple. But it seems to assume that you start up with a number (seed) that is already non-0 and has a fraction part.

The frac(997*x) PRNG indeed is considered a generator with good properties. But the results heavily rely on the seed. If I remember correctly, this seed has to provide 7 decimals after the multiplication by 997, and the last (7th) digit must be a 1, 3, 7 or 9. Choosing a suitable seed is so important that HP even recommended a special value, which was 0,5284163. The mentioned PPC article suggests 0,7742713 as an alternative.

Always starting with the same seed of course means that always the same sequence of "random" numbers is generated. Thirty years ago the suggested fix was an initialization routine with a loop that permamently generates new numbers, and let this routine run for a few seconds.

All this shows that a good PRNG with good statistical properties (which includes much more than mean and standard devation) is not trivial at all. In this regard the frac(9821x+0,211327) generator seems to perform exceptionally well.

Dieter

[Garth Wilson]

One way which has been addressed here only indirectly, through Trond's link, is using the time module if you have it. It is unknown ahead of time what the hundredths of seconds, maybe the seconds too, will be when you reach for a random number, so bringing the time into the method will be quite effective.

[Trond]

(12-05-2017 06:54 PM)Dieter Wrote:  

(12-05-2017 03:24 PM)Trond Wrote:  Wow, that is very simple. But it seems to assume that you start up with a number (seed) that is already non-0 and has a fraction part.

The frac(997*x) PRNG indeed is considered a generator with good properties. But the results heavily rely on the seed. If I remember correctly, this seed has to provide 7 decimals after the multiplication by 997, and the last (7th) digit must be a 1, 3, 7 or 9. Choosing a suitable seed is so important that HP even recommended a special value, which was 0,5284163. The mentioned PPC article suggests 0,7742713 as an alternative.

Always starting with the same seed of course means that always the same sequence of "random" numbers is generated. Thirty years ago the suggested fix was an initialization routine with a loop that permamently generates new numbers, and let this routine run for a few seconds.

All this shows that a good PRNG with good statistical properties (which includes much more than mean and standard devation) is not trivial at all. In this regard the frac(9821x+0,211327) generator seems to perform exceptionally well.

Dieter

Thanks for the info, Dieter. Yes, I was assuming that it takes more than the mean and SD to really check on a random number generator (let me know if you know further ways to check, that does not involve heavy computing power).
I guess adding 0,211327 each time would not work, as it would mess up the distribution. I suppose I could program around it if I want to try this random generator, at least avoiding that the calc starts from 0.

[Dieter]

(12-05-2017 07:50 PM)Trond Wrote:  Thanks for the info, Dieter. Yes, I was assuming that it takes more than the mean and SD to really check on a random number generator (let me know if you know further ways to check, that does not involve heavy computing power).

Take a look at the PPC journals I mentioned. In the V4N8 article about a dozen tests are mentioned.

(12-05-2017 07:50 PM)Trond Wrote:  I guess adding 0,211327 each time would not work, as it would mess up the distribution.

No, of course it doesn't! It's exactly this method (r_n+1 = 9821 · r_n+0,211327) that makes this a high-quality generator that is hard to beat.
Of course you can start with essentially any seed between 0 and 1.

Dieter

[Trond]

(12-05-2017 06:08 AM)Sylvain Cote Wrote:  HP used this random number generator in lots of places over the years.

FRAC((x * 9821) + 0.211327)

From the GAMES pac:

LBL "RNDM"
RCL 00
9821
*
.211327
+
FRC
STO 00
RTN

Sylvain

Thanks again for this one! I added one more line to avoid the program messing up in case the variable 00 happens to be negative:

LBL "RNDM"
RCL 00
ABS
9821
*
.211327
+
FRC
STO 00
RTN