Picking A, C And M For Linear Congruential Generator
Solution 1:
From Wikipedia:
Provided that c is nonzero, the LCG will have a full period for all seed values if and only if:
- c and m are relatively prime,
- a-1 is divisible by all prime factors of m,
- 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:
- 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.
- 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.
- 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):
deflcg(state):
x = 1
a = x*7*47 + 1
c = 100
m = 48**2 - 1return (a * state + c) % m
expected_period = 48**2 - 1
seeds = [5]
for i inrange(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.)
Solution 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 = 2while n > 1:
while n % d == 0:
factors.append(d)
n //= d
d += 1if d * d > n:
if n > 1: factors.append(n)
breakreturn factors
def__multiply_numbers(numbers):
"""multiply all numbers in array"""
result = 1for 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)
whilelen(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])
defcreate_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 notin 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
defnext_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 = 0for i inrange(target):
state = newState = next_number(state, m, a ,c, target)
randomSum += newState
print(checkSum == randomSum) # true
LCG is quite fascinating and usable in things like games. You can iterate a giant list of things in a deterministic random order. Shuffeling and saving the whole list is not required:
defriter(alist):
""" iterate list using LCG """
target = len(alist)
state, m, a, c = create_new_seed(target)
for _ inrange(target):
yield alist[state]
state = next_number(state, m, a ,c, target)
It is easy to save the state in between iteration steps:
savedState = '18:19:25:6:12047269:20'print('loading:', savedState)
i, state, m, a, c, target = (int(i) for i in savedState.split(':'))
state = next_number(state, m, a, c, target)
i += 1print('i:', i, 'is random number:', state, 'list done:', i+1 == target)
print('saving:', '{}:{}:{}:{}:{}:{}'.format(i, state, m, a, c, target))
Post a Comment for "Picking A, C And M For Linear Congruential Generator"