Picking A, C and M for Linear congruential generat

问题:

I am looking to implement a simple pseudorandom number generator (PRNG) that has a specified period and guaranteed no collisions for the duration of that period. After doing some research I came across the very famous LCG which is perfect. The problem is, I am having trouble understanding how to properly configure it. Here is my current implementation:

    function LCG (state)
    {
        var a = ?;
        var c = ?;
        var m = ?;

        return (a * state + c) % m;
    }

It says that in order to have a full period for all seed values the following conditions must be met:

1. c and m are relatively prime
2. a-1 is divisible by all prime factors of m
3. a-1 is a multiple of 4 if m is a multiple of 4

1 and 3 are simple to understand and test for. However what about 2, I don't quite understand what that means or how to check for it. And what about C, can it be zero? what if it's non-zero?

Overall I need to select A, C and M in such a way that I have a period of 48^5 - 1. M is equal to the period, I am not sure about A and C.

回答1:

From Wikipedia:

Provided that c is nonzero, the LCG will have a full period for all seed values if and only if:

1. c and m are relatively prime,
2. a-1 is divisible by all prime factors of m,
3. a-1 is a multiple of 4 if m is a multiple of 4.

You said you want a period of 48⁵-1, so you must choose m≥48⁵-1. Let's try choosing m=48⁵-1 and see where that takes us. The conditions from the Wikipedia article prohibit you from choosing c=0 if you want the period to be m.

Note that 11, 47, 541, and 911 are the prime factors of 48⁵-1, since they're all prime and 11*47*541*911 = 48⁵-1.

Let's go through each of those conditions:

1. For c and m to be relatively prime, c and m must have no common prime factors. So, pick any prime numbers other than 11, 47, 541, and 911, then multiply them together to choose your c.
2. You'll need to choose a such that a-1 is divisible by all the prime factors of m, i.e., a = x*11*47*541*911 + 1 for any x of your choosing.
3. Your m is not a multiple of 4, so you can ignore the third condition.

In summary:

• m = 48⁵-1,
• c = any product of primes other than 11, 47, 541, and 911 (also, c must be less than m),
• a = x*11*47*541*911 + 1, for any nonnegative x of your choice (also, a must be less than m).

---

Here's a smaller test case (in Python) using a period of 48²-1 (which has prime factors 7 and 47):

def lcg(state):
    x = 1
    a = x*7*47 + 1
    c = 100
    m = 48**2 - 1
    return (a * state + c) % m

expected_period = 48**2 - 1
seeds = [5]
for i in range(expected_period):
    seeds.append(lcg(seeds[-1]))
print(len(set(seeds)) == expected_period)

It outputs True, as it should. (If you have any trouble reading Python, let me know and I can translate it to JavaScript.)

回答2:

Based on Snowball's answer and the comments I've created a complete example. You can use the set == list comparison for smaller numbers. I could not fit 48^5-1 into memory.

To circumvent the a < m problem, I'm incrementing the target a few times to find a number where a is able to be < m (where m has duplicated prime factors). Surprisingly +2 is enough for a lot of numbers. The few extra numbers are later skipped while iterating.

import random

def __prime_factors(n):
    """
    https://stackoverflow.com/a/412942/6078370
    Returns all the prime factors of a positive integer
    """
    factors = []
    d = 2
    while n > 1:
        while n % d == 0:
            factors.append(d)
            n //= d
        d += 1
        if d * d > n:
            if n > 1: factors.append(n)
            break
    return factors

def __multiply_numbers(numbers):
    """multiply all numbers in array"""
    result = 1
    for n in numbers:
        result *= n
    return result

def __next_good_number(start):
    """
    https://en.wikipedia.org/wiki/Linear_congruential_generator#c%E2%89%A00
    some conditions apply for good/easy rotation
    """
    number = start
    factors = __prime_factors(number)
    while len(set(factors)) == len(factors) or number % 4 == 0:
        number += 1
        factors = __prime_factors(number)
    return number, set(factors)

# primes < 100 for coprime calculation. add more if your target is large
PRIMES = set([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
        43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97])

def create_new_seed(target):
    """be aware, m might become > target"""
    m, factors = __next_good_number(target)
    a = __multiply_numbers(factors) + 1
    # https://en.wikipedia.org/wiki/Coprime_integers
    otherPrimes = [p for p in PRIMES if p not in factors]
    # the actual random part to get differnt results
    random.shuffle(otherPrimes)
    # I just used arbitary 3 of the other primes
    c = __multiply_numbers(otherPrimes[:3])
    # first number
    state = random.randint(0, target-1)
    return state, m, a, c

def next_number(state, m, a ,c, limit):
    newState = (a * state + c) % m
    # skip out of range (__next_good_number increases original target)
    while newState >= limit:
        newState = (a * newState + c) % m

    return newState

if __name__ == "__main__":
    target = 48**5-1
    state, m, a, c = create_new_seed(target)
    print(state, m, a, c, 'target', target)

    # list and set can't fit into 16GB of memory
    checkSum = sum(range(target))
    randomSum = 0
    for i in range(target):
        state = newState = next_number(state, m, a ,c, target)
        randomSum += newState
    print(checkSum == randomSum) # true