find more indipendent seed value for a 64 bit LCG (MMIX (by Knuth))

Question

I'm using a 64 bit LCG (MMIX (by Knuth)). It generate a certain block of random numbers inside my code, which use them to perform some operations. My code works in single core and I would like to parallelize the work to reduce the execution time.
Before start thinking to more advanced methods in this sense I'd like to simply execute more identical codes in parallel (in fact the code repeats the same task over a certain numbers of indipendent simulation, so I can simply split the number of simulation between more identical codes and run them in parallel).
My only problem now is to find a seed for each code; in particular, to avoid the possibility of unwanted non trivial correlation between data generated in different codes, I have to be sure that the random number generated in the various codes don't overlap. To do so, starting from a certain seed in the first code I have to find a way to find a value (the next seed) very distant not in absolute value but in the pseudo-random sequence (so, such that, to go from the first to the second seed, I need a huge number of steps of LCG).
My first attempt was this:
starting from the LCG relation between 2 consecutive numbers generated in the sequence

So, in principle, I could calculate the above relation with, say, n = 2^40 and I_0 equal to the value of the first seed, and obtain a new seed distant 2^40 steps in the random CLG sequence from the first one.
The problem is that, doing so, I necessary go in overflow calculating a^n. In fact for MMIX (by Knuth) a~2^62 and i use unsigned long long int (64 bit). Note that the only problem here is the fraction in the above relation. If there only were sum and multiplication I could ignore the overflow problem due to the following modular properties (in fact I'm using 2^64 as c (64 bit generator)):

So, starting from a certain value (first seed), how can I find a second one distant a huge number of step in the LC pseudo-random sequence?

[EDIT]
r3mainer solution is perfectly suited for python codes. I'm trying now to implement it in c using unsigned __int128 variables. I have only one problem: in principle I should compute:

Say, for simplicity, I want to compute:

with n = 2^40 and c(a-1)~2^126. I proceed with a cycle.Starting with temp = a, in each iteration I compute temp = temp*temp, then I compute temp%c(a-1). The problem is in the second step (temp = temp*temp). temp in fact could be, in principle any number < c(a-1)~2^126. If temp is a big number, say > 2^64, I'll go in overflow, reaching 2^128 - 1, before the next module operation. So is there a way to avoid it? For now the only solution I see is to perform each multiplication with a loop over bit, as suggested here: c code: prevent overflow in modular operation with huge modules (modules near the overflow treshold) Is there another way to perform module operation during the multiplication?
(note that being c = 2^64, with mod(c) operation I don't have the same problem because the overflow point (for ull int variables) coincides with the module)

user_1679849 · Answer

Any LCG of the form x[n+1] = (x[n] * a + c) % m can be skipped to an arbitrary position very quickly.

Starting with a seed value of zero, the first few iterations of the LCG will give you this sequence:

x₀ = 0
x₁ = c % m
x₂ = (c(a + 1)) % m
x₃ = (c(a² + a + 1)) % m
x₄ = (c(a³ + a² + a + 1)) % m

It's pretty easy to see that each term is actually the sum of a geometric series, which can be calculated with a simple formula:

x_n = (c(a^{n-1} + a^{n-2} + ... + a + 1)) % m
    = (c * (a^n - 1) / (a - 1)) % m

The (a^n - 1) term can be calculated quickly by modular exponentiation, but dividing by (a-1) is a bit tricky because (a-1) and m are both even (i.e., not coprime), so we can't calculate the modular multiplicative inverse of (a-1) mod m directly.

Instead, calculate (a^n-1) mod m*(a-1), then perform a straightforward (non-modular) division of the result by a-1. In Python, the calculation would go something like this:

def lcg_skip(m, a, c, n):
    # Calculate nth term of LCG sequence with parameters m (modulus),
    # a (multiplier) and c (increment), assuming an initial seed of zero
    a1 = a - 1
    t = pow(a, n, m * a1) - 1
    t = (t * c // a1) % m
    return t

def test(nsteps):
    m = 2**64
    a = 6364136223846793005
    c = 1442695040888963407
    #
    print("Calculating by brute force:")
    seed = 0
    for i in range(nsteps):
        seed = (seed * a + c) % m
    print(seed)
    #
    print("Calculating by fast method:")
    # Calculate nth term by modular exponentiation
    print(lcg_skip(m, a, c, nsteps))

test(1000000)

So to create LCGs with non-overlapping output sequences, all you would need to do is use initial seed values generated by lcg_skip() with values of n that are far enough apart.

user_4044696 · Answer

Well, for LCG it is known property to jump forward and backward in O(log₂(N)) time where N is the distance between jump points, paper by F. Brown, "Random Number Generation with Arbitrary Stride," Trans. Am. Nucl. Soc. (Nov. 1994).

It means if you have LCG parameters (a, c) satisfying Hull–Dobell Theorem, then whole period would be 2⁶⁴ numbers before repeating themself, and say for Nt number pf threads you make jump distance of 2⁶⁴ / Nt, and all threads start with the same seed and just jump after initializing LCG by (2⁶⁴ / Nt)*threadId, and you would be completely safe from RNG correlations due to sequences overlap.

For simplest case of common 64 unsigned modulo math, as implemented in NumPy, code below should work fine

import numpy as np

class LCG(object):

    UZERO: np.uint64 = np.uint64(0)
    UONE : np.uint64 = np.uint64(1)

    def __init__(self, seed: np.uint64, a: np.uint64, c: np.uint64) -> None:
        self._seed: np.uint64 = np.uint64(seed)
        self._a   : np.uint64 = np.uint64(a)
        self._c   : np.uint64 = np.uint64(c)

    def next(self) -> np.uint64:
        self._seed = self._a * self._seed + self._c
        return self._seed

    def seed(self) -> np.uint64:
        return self._seed

    def set_seed(self, seed: np.uint64) -> np.uint64:
        self._seed = seed

    def skip(self, ns: np.int64) -> None:
        """
        Signed argument - skip forward as well as backward

        The algorithm here to determine the parameters used to skip ahead is
        described in the paper F. Brown, "Random Number Generation with Arbitrary Stride,"
        Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
        O(log2(N)) operations instead of O(N). It computes parameters
        A and C which can then be used to find x_N = A*x_0 + C mod 2^M.
        """

        nskip: np.uint64 = np.uint64(ns)

        a: np.uint64 = self._a
        c: np.uint64 = self._c

        a_next: np.uint64 = LCG.UONE
        c_next: np.uint64 = LCG.UZERO

        while nskip > LCG.UZERO:
            if (nskip & LCG.UONE) != LCG.UZERO:
                a_next = a_next * a
                c_next = c_next * a + c

            c = (a + LCG.UONE) * c
            a = a * a

            nskip = nskip >> LCG.UONE

        self._seed = a_next * self._seed + c_next


#%%
np.seterr(over='ignore')

seed = np.uint64(1)

rng64 = LCG(seed, np.uint64(6364136223846793005), np.uint64(1))

print(rng64.next())
print(rng64.next())
print(rng64.next())

#%%
rng64.skip(-3) # back by 3
print(rng64.next())
print(rng64.next())
print(rng64.next())

rng64.skip(-3) # back by 3
rng64.skip(2) # forward by 2
print(rng64.next())

Tested in Python 3.9.1, x64 Win 10

find more indipendent seed value for a 64 bit LCG (MMIX (by Knuth))

2 Answers