Hp Algorithm to create random numbers

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