Hp Algorithm to create random numbers

[Gil]

When pressing on HP50G
NumberX RDZ... RAND, we get always the same random number, whatever... is.
That's OK and that's what we want.

Question.
Does somebody know what is the algorithm (the internal special function) used by HP calculators to generate pseudo random numbers ?

Thanks for your answer.

[Thomas Klemm]

(08-12-2024 11:42 AM)Gil Wrote:  Does somebody know what is the algorithm (the internal special function) used by HP calculators to generate pseudo random numbers ?

• HP-48: How RAND Works
• HP-15C: A Python program that simulates the RAN# function
• 50g: the random generator source code
• HP's original annotated source code for the HP 48's RAND and RANDOMIZE functions
  

[Gil]

Thanks for the reference, Thomas.

Regards,
Gil

[Gil]

With the first reference
I press, on my HP50G, the command RAND
and get R(n) =.821118792263 (1)
and build R(n)_int =821118792263 (an integer) (2)

Then I multiply 2851130928467 (3)
by R(n) and get 2.34111718457E12 (4)

Then I multiply also 2851130928467 (3)
by R(n)_int and get 2341117184566508886050821 (5)

Pressing a second time RAND gives
R(n+1)=.994325264552 (6)

I thought, with the explanations of the first reference,
that the digits of (6), R(n+1), should be found somewhere in (5) or (6).

But it's clearly not the case.

Could somebody explain my error or give a detailed example?

Thanks in advance.

Regards,
Gil

[Albert Chan]

Hi, Gil

It is a difference of random state vs random output

Quote:R(n) =.821118792263 (1)

Output is 12 digits, but random state is really 15 (red numbers cracked by brute force)

821118792263521 * 2851130928467 % 10^15 = 994325264552307

We drop next state last 3 digits for next random output.

Quote:Pressing a second time RAND gives
R(n+1)=.994325264552 (6)

[DavidM]

As Albert also indicated, the internal seed is actually 15 decimal digits. Here's another series of examples that might help to show what's happening.

First, I'll set the seed to a nice number with RDZ so that this can be easily replicated:

1E-13 RDZ

The internal seed now has a value of 000000000000001.

Executing RAND at this stage causes the following:
new_seed = (000000000000001 x 2851130928467) MOD 10^15 = 002851130928467
RAND result is normalized to 2.85113092846E-3 *

Next RAND:
new_seed = (002851130928467 x 2851130928467) MOD 10^15 = 261097470970089
RAND result is normalized to 0.26109747097

Next RAND:
new_seed = (261097470970089 x 2851130928467) MOD 10^15 = 335429755623563
RAND result is normalized to 0.335429755623

Next RAND:
new_seed = (335429755623563 x 2851130928467) MOD 10^15 = 468090732667921
RAND result is normalized to 0.468090732667

* The normalization step appears to shift out any leading zeros before truncating the last three digits.

[Gil]

Thanks to Albert for his cracking effort
and to David for his detailed examples and explanations.

Gil

[Gil]

And what is the algorithm/function to get the first RANDom number when entering x RDZ RAND (x a real number)?

Thanks for partaking this information.

Regards,
Gil

[DavidM]

(08-19-2024 10:23 PM)Gil Wrote:  And what is the algorithm/function to get the first RANDom number when entering x RDZ RAND (x a real number)?

RAND is no different in this scenario than previously described, so what you're really asking is "how does RDZ compute a seed?", which turns out to be difficult to explain without some Saturn assembly knowledge.

The Saturn code for RDZ is concise and efficient, but it doesn't map well to higher-level mathematical translation because it is heavily dependent on the internal representation of reals and "bit-twiddling" steps that are much easier to do in Saturn code than they are with traditional numeric calculations.

That said, I've attempted to create a User RPL version of RDZ that returns an exact integer seed for the given basis, leaving it on the stack instead of storing it internally. It treats an input of 0 in the same way that RDZ does, ie. using the system clock for the seed basis. This code requires a 49g/49g+/50g due to its use of exact integers. The use of reals and integers is important here, so the program should be transferred while the destination calculator is in exact mode.

Code:

RDZU
\<<
  @ if given 0, use the low-order 20 bits of the system clock instead
  IF
    DUP NOT
  THEN
    TICKS NIP
    RCWS 20 STWS
    SWAP B\->R SWAP
    STWS
  END

  @ we need a non-negative real number for the seed basis
  I\->R ABS

  @ separate the mantissa and exponent fields as integers
  DUP MANT 1E14 * R\->I
  SWAP XPON 1000. + 1000. MOD R\->I

  @ denormalize
  DO
    1 +
    CASE
      OVER NOT THEN 1. END
      DUP 2 * 1000 IDIV2 SWAP NOT THEN DROP 1. END
      NOT THEN 1. END
      SWAP 10 IQUOT SWAP
      0.
    END
  UNTIL
  END

  @ replace exponent field after left-shift/increment
  10 * 1 + 1000 MOD
  SWAP 1000 IQUOT 1000 *
  +
\>>

Note that the resulting seed should be interpreted as having leading 0s if its length is less than 15 digits. So a result of 999001 should be interpreted as 000000000999001 when stored internally.

Some examples

1E-13 RDZU => 1
12345. RDZU => 123450000000051
9.87654321098E295 RDZU => 987654321098961
0.000000000999 RDZU => 999001

Full RAND check

12345 RDZU        => 123450000000051
2851130928467 *   => 351972113119396557677351817
10 15 ^ MOD       => 396557677351817
...so RAND should be 0.396557677351

12345 RDZ RAND    => 0.396557677351



9.99E-10 RDZU     => 999001
2851130928467 *   => 2848282648669461467
10 15 ^ MOD       => 282648669461467
...so RAND should be 0.282648669461

9.99E-10 RDZ RAND => 0.282648669461

[Gil]

A great thank to you, David.

I will study your programs later with calm.

Regards,
Gil

[Albert Chan]

(08-21-2024 05:43 PM)DavidM Wrote:  Some examples

1E-13 RDZU => 1
12345. RDZU => 123450000000051
9.87654321098E295 RDZU => 987654321098961
0.000000000999 RDZU => 999001

RDZU seems to use 12 digits mantissa + 2 digits exponents (by MOD 100) + final 1 to make 15 digits

Code:

.000000 000000 1 E000   --> 000000 000000 001
.123450 000000   E005   --> 123450 000000 051
.987654 321098   E296   --> 987654 321098 961
.000000 000999   E000   --> 000000 000999 001

Mantissa first 12 digits is zero only when seed is small.
I would expect smaller seed (< 1E-13) would adjust exponent field. (to get state ??1)
This was my guess.

1E-13 --> 001
1E-14 --> 991
1E-15 --> 981
1E-16 --> 971
...

But when I try RDZ, above 4 cases all give random state 001
Only when seed get smaller, does this state pattern activated, why the gap?

1E-13 --> 001
1E-14 --> 001
1E-15 --> 001
1E-16 --> 001
1E-17 --> 991
1E-18 --> 981
1E-19 --> 971
...

[DavidM]

(08-23-2024 01:42 PM)Albert Chan Wrote:  But when I try RDZ, above 4 cases all give random state 001
Only when seed get smaller, does this state pattern activated, why the gap?

1E-13 --> 001
1E-14 --> 001
1E-15 --> 001
1E-16 --> 001
1E-17 --> 991
1E-18 --> 981
1E-19 --> 971

There's more to RDX's algorithm than your table implies.

RDX separates the mantissa and exponent fields of the real number given as input, and then treats the mantissa as a 15-digit integer while combining the two into a seed. The exponent stays as a 3-digit number, but it may not be what you expect due to the special way negative exponents for reals are encoded in the Saturn architecture. Specifically, negative exponents are encoded as 1000-ABS(exponent), but RDX simply treats that result as a 3-digit positive integer.

The denormalization section of the code performs right-shifts on the mantissa based on the exponent value.

Note that there are 3 independent exit conditions in RDX's denormalization loop:
- the mantissa becomes 0 as a result of repeated right-shifts, OR
- (2*exponent) <= 998 (ie. a carry does NOT occur as a result of 2*exponent), OR
- (2*exponent) MOD 1000 is 0

So here's the critical values for some specific inputs that will hopefully show what's happening.

Input 1E-12
Initial mantissa: 100000000000000
Initial exponent:             988
Mantissa after denormalization: 000000000001000
Exponent after denormalization:             000
Exponent after shift-left:                  000
Exponent after increment:                   001
Seed:                           000000000001001

Input 1E-13
Initial mantissa: 100000000000000
Initial exponent:             987
Mantissa after denormalization: 000000000000100
Exponent after denormalization:             000
Exponent after shift-left:                  000
Exponent after increment:                   001
Seed:                           000000000000001

Input 1E-14
Initial mantissa: 100000000000000
Initial exponent:             986
Mantissa after denormalization: 000000000000010
Exponent after denormalization:             000
Exponent after shift-left:                  000
Exponent after increment:                   001
Seed:                           000000000000001

Input 1E-15
Initial mantissa: 100000000000000
Initial exponent:             985
Mantissa after denormalization: 000000000000001
Exponent after denormalization:             000
Exponent after shift-left:                  000
Exponent after increment:                   001
Seed:                           000000000000001

Input 1E-16
Initial mantissa: 100000000000000
Initial exponent:             984
Mantissa after denormalization: 000000000000000
Exponent after denormalization:             000
Exponent after shift-left:                  000
Exponent after increment:                   001
Seed:                           000000000000001

Input 1E-17
Initial mantissa: 100000000000000
Initial exponent:             983
Mantissa after denormalization: 000000000000000
Exponent after denormalization:             999
Exponent after shift-left:                  990
Exponent after increment:                   991
Seed:                           000000000000991

Input 123456.123456
Initial mantissa: 123456123456000
Initial exponent:             005
Mantissa after denormalization: 123456123456000
Exponent after denormalization:             006
Exponent after shift-left:                  060
Exponent after increment:                   061
Seed:                           123456123456061

Note how the last 3 digits of the mantissa get obliterated by the altered exponent in each example.

[Albert Chan]

(08-24-2024 12:05 AM)DavidM Wrote:  RDX separates the mantissa and exponent fields of the real number given as input,
and then treats the mantissa as a 15-digit integer while combining the two into a seed.

Thanks, 15 digits mantissa make sense.
However, 16 digits mantissa is simpler to understand RDZ logic.

Code:

1E-12 = .0000 0000 0001 0000 E000   --> 0000 0000 0001 00 1
1E-13 = .0000 0000 0000 1000 E000   --> 0000 0000 0000 00 1
1E-14 = .0000 0000 0000 0100 E000   --> 0000 0000 0000 00 1
1E-15 = .0000 0000 0000 0010 E000   --> 0000 0000 0000 00 1
1E-16 = .0000 0000 0000 0001 E000   --> 0000 0000 0000 00 1
1E-17 = .0000 0000 0000 0001 E-01   --> 0000 0000 0000 99 1    ; -01 % 100 = 99

[Gil]

Despite all your detailed explanations, David, would it be possible to work out a detailed example, again step by step, to get the RAND om number .510262440213
after entering 2.567 RDZ on the HP50G?

And perhaps also, if I may ask, for
RANDom number .163041130928
after entering 0.002567 RDZ?

Thanks in advance for your painstaking, David.

Regards,
Gil

[DavidM]

(08-25-2024 10:52 AM)Gil Wrote:  Despite all your detailed explanations, David, would it be possible to work out a detailed example, again step by step, to get the RAND om number .510262440213
after entering 2.567 RDZ on the HP50G?

And perhaps also, if I may ask, for
RANDom number .163041130928
after entering 0.002567 RDZ?

Thanks in advance for your painstaking, David.

Regards,
Gil

In your first example, entering 2.567 RDZ produces an internal seed value of 256700000000011.
When executing RAND immediately following that, we have the following:
new_seed_value = (256700000000011 * 2851130928467) MOD 10^15 = 510262440213137

The normalization process for that integer yields 0.510262440213 for the final RAND result.

In your second example, it appears that you may have inadvertently left out the "6" when executing "0.002567 RDZ":
0.00257 RDZ RAND => 0.163041130928
0.002567 RDZ RAND => 0.377640130928

The only reason I figured that out was that I did the same exact thing when testing this after reading your post -- it had me confused for a minute or so until I subsequently typed it in correctly on a second try.

So here's the same process for that input.

Entering 0.002567 RDZ produces an internal seed value of 2567000000001.
When executing RAND immediately following that, we have the following:
new_seed_value = (2567000000001 * 2851130928467) MOD 10^15 = 377640130928467

The normalization process for that integer yields 0.377640130928 for the final RAND result.

Does this help?

[Albert Chan]

(08-25-2024 12:02 PM)DavidM Wrote:  In your second example, it appears that you may have inadvertently left out the "6" when executing "0.002567 RDZ":
0.00257 RDZ RAND => 0.163041130928
0.002567 RDZ RAND => 0.377640130928

The only reason I figured that out was that I did the same exact thing when testing this after reading your post -- it had me confused for a minute or so until I subsequently typed it in correctly on a second try.

Assuming 1E-16 ≤ seed < 1, we know original state must end in 001

--> state * 2851130928467 ≡ 163041130928467 (mod 10^15)

1/2851130928467 ≡ -46007610876197 ≡ 953992389123803 (mod 10^15)      (*)

--> state ≡ 163041130928467 * 953992389123803 ≡ 002570000000001 (mod 10^15)
--> seed = 0.00257

---

(*) Calculations based on Euclidean GCD calculations, then scale-up for mod inverse, see here

1E15            --> -floor(446607260/9707247377 * 1E15) = -46007610876197
2851130928467 +
2104175036550 -
746955891917 +
610263252716 -
136692639201 +
63492695912 -
9707247377 --> +floor(13843/300885 * 9707247377) = 446607260, 9.71E9*6.35E10 ≈ 6.2E20
5249211650 -
4458035727 +
791175923 -
502156112 +
289019811 -
213136301 +
75883510 -
61369281 +
14514229 -
3312365 +
1264769 -
782827 +
481942 -
300885 --> -floor(-19/413 * 300885) = 13843, 3.01E5*4.82E5 ≈ 1.4E10
181057 +
119828 -
61229 +
58599 -
2630 +
739 -
413 --> -floor(1/21 * 413) = -19, 413*739 = 305207
326 +
87 -
65 +
22 -
21 --> 1/1 ≡ 1 (mod 21), 21*22 = 462
1 = gcd

[Gil]

Thanks for your answer, David.

The seed is OK and RAND calculation also.

But could you detail the internal_seed with 2.567.
—> 257....001?

Thanks again in advance, David.

Gil

[DavidM]

(08-25-2024 05:51 PM)Gil Wrote:  But could you detail the internal_seed with 2.567.
—> 257....001?

The computation of RAND is a straightforward chain of mathematical operations that's easy to describe. The conversion of a real number to a random seed (RDZ) isn't as easy to explain, mainly due to the variability of what happens during the denormalization process. Your best bet would be to step through the code with DBUG so you could watch the stack contents change with various inputs.

At a high level, RDZ/RDZU works as follows (using 2.567 as the example here):

1) The input is altered to conform to standard (non-zero positive real): 2.567E0

2) That value is separated into two intermediate pieces. mantissa (256700000000000) and exponent (0)

3) Denormalization occurs (the actual steps here vary with the input). For this example, the end result is:

Code:

2:      256700000000000       (mantissa)
1:                  001       (exponent)

4) exponent is then shifted left 1 digit:

Code:

2:      256700000000000       (mantissa)
1:                  010       (exponent)

5) exponent is incremented by 1:

Code:

2:      256700000000000       (mantissa)
1:                  011       (exponent)

6) Finally, the last three digits of mantissa are replaced by the full 3-digit form of exponent:

Code:

1:      256700000000011       (final seed value)

[Albert Chan]

(08-25-2024 02:50 PM)Albert Chan Wrote:  Assuming 1E-16 ≤ seed < 1, we know original state must end in 001

--> state * 2851130928467 ≡ 163041130928467 (mod 10^15)

It may be better to reduce modulus from huge 10^15 down to size, to simplify mod inverse calculations.

Let state = 1000 * x + 001, where x is 12 digits integer

x * 2851130928467 ≡ (163041130928467 - 2851130928467) / 1000 ≡ 16019 * 10^7 (mod 10^12)

Last 7 digits of x must be 0. Let x = 10^7 * y, where y is 5 digits integer

y * 28467 ≡ 16019 (mod 10^5)

Again, Euclidean GCD, follow by scaling up to get modulo inverse.

10^5 --> +floor(3475/14599 * 10^5) = 23803
28467 -
14599 --> +floor(5/21 * 14599) = 3475, 14599*28467 > 10^5
13868 -
731 +
710 -
21 --> +floor(1/4 * 21) = 5, 21*710 = 14910
17 -
4 --> 1/1 = 1 (mod 4), 4*17 = 68
1 = gcd

y ≡ 16019 / 28467 ≡ 16019 * 23803 ≡ 00257 (mod 10^5)

--> state = 00257 0000000 001
--> seed = 0.00257

---

Another way, by checking last digits

y * 28467 ≡ 16019 (mod 10^5)      // multiply both side by 3, to make 7 to 1
y * 85401 ≡ 48057 (mod 10^5)      // y must end in 57

Let y = 100*z + 57

z * 85401 ≡ (48057 - 57*85401) / 100 ≡ -48198 (mod 10^3)
z * 401 ≡ 802 (mod 10^3)             // z must end in 02
z ≡ 802/401 ≡ 002 (mod 10^3)

--> y ≡ 002 57 (mod 10^5)

[Albert Chan]

(08-25-2024 02:50 PM)Albert Chan Wrote:  state * 2851130928467 ≡ 163041130928467 (mod 10^15)

Faster way to solve, scale to reduce multiplier to 1, 2nd column is state.

Cas> m := 10^15
Cas> [2851130928467, 163041130928467] * 3; /* force multiplier last digit = 1 */

[8553392785401,489123392785401]

Cas> (Ans*(2-Ans[1])) MOD m

[667229546840001,469799546840001]
[246014400000001,248584400000001]
[1,2570000000001]

--> state = 002570000000001

This trick work for other modulus too. Also, MOD can be signed mod

Example: solve 7x ≡ 9 (mod 3^5)
Base 3, last digit 7%3 = 1, we keep scaling to zeroed out other digits.

Cas> [7, 9] %% 3^5
Cas> Ans*(2-Ans[1])
[-35 %% 243,-45 %% 243]
[-80 %% 243,36 %% 243]
[1 %% 243,36 %% 243]

--> x ≡ 9/7 ≡ 36 (mod 3^5)