zIFBoards - Free Forum Hosting
Free Forums with no limits on posts or members.

Learn More · Register Now
Welcome to Dozensonline. We hope you enjoy your visit.
You're currently viewing our forum as a guest. This means you are limited to certain areas of the board and there are some features you can't use. If you join our community, you'll be able to access member-only sections, and use many member-only features such as customizing your profile, and sending personal messages. Registration is simple, fast, and completely free. (You will be asked to confirm your email address before we sign you on.)
Join our community!
If you're already a member please log in to your account to access all of our features:

Name:   Password:


 

 Optimizing The Linear Congruential Generator, for pseudo-random numbers
Dan
Posted: Nov 7 2014, 04:41 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



This topic is highly related to Optimal base for random guessing by humans, but I decided to create a new thread because I wanted to (1) focus on one particular algorithm, and (2) discuss "computer-scale" parameters rather than just "human-scale" ones.


The Linear Congruential Generator (LCG) is one of the most common algorithms for deterministic generation of "random" numbers. It is defined by (in Python syntax):

CODE
def random_sequence(multiplier, increment, modulus, seed):
   value = seed
   for index in range(modulus):
       value = (multiplier * value + increment) % modulus
       yield value


(The sequence is periodic, so technically the "for" statement should be "while True", but the sequence will be easier to analyze by treating it as finite.)

Each value of the sequence is an integer between 0 (inclusive) and the modulus (exclusive).

For dozenal counting, an obvious choice for the modulus is 12 (*10). This yields a sequence of single-digit numbers. The hard part is deciding on the multiplier and increment.

Some choices are quite obviously bad. For example, the multiplier a=3 and increment c=3 (starting with a seed of 0) gives the sequence [3, 0, 3, 0, 3, 0, ...]. This is not only way too obvious of a pattern, but only uses two of the twelve possible values.

The most basic desirable property of an LCG is to have a full period of m values, each used once per period. This ensures a uniform distribution of values.

With m=12, there are only four possible combinations of parameters meeting this criterion:
  • a=1, c=1: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, ...
  • a=1, c=5: 5, 10, 3, 8, 1, 6, 11, 4, 9, 2, 7, 0, ...
  • a=1, c=7: 7, 2, 9, 4, 11, 6, 1, 8, 3, 10, 5, 0, ...
  • a=1, c=11: 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, ...

But not of these sequences is "random enough" (a concept I'll more formally define in a later post). So, m=12 just isn't a good choice.

So, let's try two dozenal digits. With m=144 (*100), there are 576 (*400) possible combinations of a and c to chose from. How do we choose?

To be continued...
Top
Dan
Posted: Nov 17 2014, 12:19 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



It's hard to find an example of a good random number generator, but it's easy to find an example of a bad one.

Let Xi = The day of the week of January 1, 2000+i (Sunday=0, ..., Saturday=6)

This isn't an LCG, but it is periodic, with a period of n=400. Computing this "random" sequence over the full period gives:

ValueFreq.Rel. Freq.
0580.1450
1560.1400
2580.1450
3570.1425
4570.1425
5580.1450
6560.1400


The relative frequencies are close to the naively-expected 1/7 = 0.142 857 142 857..., deviating only because of a mathematical quirk of the Gregorian leap-day adjustment. (On the Julian calendar, they'd be exactly 1/7.)

To quantify this, define the disuniformity of a periodic PRNG as the sum of (actual - expected relative frequency)^2 over a whole number of full periods. (Note that this is similar to the chi-squared test statistic, but is independent of the number of periods taken.) So, for this particular "random" number generator, the disuniformity is 2*(0.14-1/7)^2 + 2*(0.1425 - 1/7)^2 + 3*(0.145-1/7)^2 = 17/560000 = 0.0000303 571428 571428...

Again, this is close to uniform, and would be a passable PRNG if one-dimensional uniformity were the only criterion. But random numbers aren't always used one at a time; they're often used in pairs. For example, the roll of a pair of dice in a board game, or the generation of random points in a plane.

In theory, taking two random numbers from the range 0-6 gives 49 possible results. But this particular PRNG makes only 14 of them possible:

ValueFreq.Rel. Freq.
(0, 1)430.1075
(0, 2)150.0375
(1, 2)430.1075
(1, 3)130.0325
(2, 3)440.1100
(2, 4)140.0350
(3, 4)430.1075
(3, 5)140.0350
(4, 5)440.1100
(4, 6)130.0325
(5, 0)150.0375
(5, 6)430.1075
(6, 0)430.1075
(6, 1)130.0325


This is because for all 7 possible days a year can start on, there are only 2 possibilities for what day the next year can start on, depending on whether the year is a leap year or not.

The disuniformity of this PRNG applied to pairs of values works out to 12646433/192080000 = 0.06583940545605997.
Top
Dan
Posted: Mar 17 2016, 02:56 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



For a general-purpose random-number generator, we can't assume that the range will be the same as the modulus. Sometimes, we'll want to flip a (2-sided) coin. Sometimes, we'll want to roll a 6-sided die. Sometimes, we'll want to draw a card from a deck of 52.

In computer programming, there are two main ways of generating a range(k) random number from a range(m) generator.

One is the C-style way, rand() % k. But this is a bad approach, especially in combination with the LCG method.

The other is the BASIC way, INT(RND() * k), where RND is a function that returns a floating-point number between 0 (inclusive) and 1 (exclusive). We can convert an integer in range(m) to such a number simply by dividing by m. So, given a generator that returns a random number r in range(m), we can scale it to range(k) by doing INT(k * r / m). If we use a floored integer division operator (like Python's //), then we can simply write k * r // m and avoid troublesome floating-point arithmetic altogether.

This still gives "subtly non-uniform" results, but at least it avoids the systematic bias towards smaller numbers that rand() % k has when k is not a divisor of m.

Top
Dan
Posted: Mar 17 2016, 03:41 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



As illustrated by my second post in this thread, it's not enough for individual "random" numbers to be uniformly distributed. We also need sequences of random numbers to be uniformly distributed as much as possible.

But which sequences? We can not reasonably expect, for example, pairs of numbers from range(144) to be uniformly distributed for the 20 736 possible values, if our generator is to have only 144 possible states.

But we can expect pairs of numbers from range(12) to be uniformly-distributed. Or for any factor of 12. For other numbers less than 12, we can't ensure exact uniformity, but we can strive to minimize the disuniformity (as defined above).

By similar reasoning, we should have the 125 possible triples of numbers from range(5) be as uniformly distributed as possible. (5 being the floor of the cube root of 144). And for quadruples from range(3) and 5-to-7-tuples from range(2).
Top
Dan
Posted: Mar 20 2016, 03:11 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



Another common use case of random numbers is shuffling a deck of cards. The standard algorithm for this is the Fisher–Yates shuffle, which for a deck of n cards requires n-1 random integers: the first from range(n), the second from range(n-1), etc., down to range(2). Given an ideal random number generator, all n! permutations should be possible.

Because 6!=720 but our proposed PRNG has only 144 possible states, there will be unreachable permutations of 6 or more cards. So let's evaluate the shuffles of 5 cards or less. We can ignore the degenerate cases of "shuffling" 1 card (which is a no-op) or 2 cards (which is equivalent to a coin toss).
Top
jrus
Posted: Mar 20 2016, 11:43 PM


Regular


Group: Members
Posts: 195
Member No.: 1,156
Joined: 23-October 15



Forget Linear Congruential Generators. Just use http://www.pcg-random.org :-)
Top
Dan
Posted: Mar 22 2016, 12:23 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



CODE
uint32_t pcg32_random_r(pcg32_random_t* rng)
{
   uint64_t oldstate = rng->state;
   // Advance internal state
   rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
   // Calculate output function (XSH RR), uses old state for max ILP
   uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
   uint32_t rot = oldstate >> 59u;
   return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}


Nice simple algorithm. Looks like an enhanced version of Xorshift. I'll have to try it the next time I do a Monte Carlo simulation.

But this thread was supposed to be about human-generated random numbers. Bitshifts and XOR may be fast for a computer, but not so much for pencil-and-paper (unless you're already using a power-of-two base).
Top
jrus
Posted: Mar 23 2016, 03:12 AM


Regular


Group: Members
Posts: 195
Member No.: 1,156
Joined: 23-October 15



Do you have some context where you need to generate random base-twelve numbers by hand, or is this just a fun curiosity?
Top
Dan
Posted: Mar 23 2016, 03:13 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



In summary, the tests I have proposed for a 144-state random number generator are:

Test# CombinationsMinimum disuniformity
Flip 2 coins40
Flip 3 coins80
Flip 4 coins160
Flip 5 coins321/2592
Flip 6 coins641/1728
Flip 7 coins1287/10368
Roll 2 d390
Roll 3 d3271/3456
Roll 4 d3817/10368
Roll 2 d4160
Roll 3 d4641/1728
Roll 2 d52519/86400
Roll 3 d51251007/1296000
Roll 2 d6360
Roll 2 d74923/169344
Roll 2 d8641/1728
Roll 2 d981 7/10368
Roll 2 d1010077/64800
Roll 2 d111211127/1254528
Roll 2 d121440
Shuffle 3 cards60
Shuffle 4 cards240
Shuffle 5 cards1201/1080


So, the question is: Is there any combination of multiplier and increment that passes all of these tests?
Top
Dan
Posted: Apr 12 2016, 02:39 AM


Dozens Disciple


Group: Members
Posts: 1,463
Member No.: 19
Joined: 8-August 05



The answer turns out to be "no". So the question is: How close can we get?

Test Minimum disuniformity Combinations of (multiplier, increment)
Flip 2 coins 1/5184 = 0.00019290123456790122 (1, 35), (1, 37), (1, 107), (1, 109), (13, 1), (13, 11), (13, 13), (13, 23), (13, 25), (13, 29), (13, 35), (13, 37), (13, 41), (13, 43), (13, 47), (13, 49), (13, 55), (13, 59), (13, 61), (13, 71), (13, 73), (13, 83), (13, 85), (13, 95), (13, 97), (13, 101), (13, 107), (13, 109), (13, 113), (13, 115), (13, 119), (13, 121), (13, 127), (13, 131), (13, 133), (13, 143), (25, 1), (25, 5), (25, 7), (25, 11), (25, 13), (25, 17), (25, 19), (25, 23), (25, 25), (25, 29), (25, 31), (25, 35), (25, 37), (25, 41), (25, 43), (25, 47), (25, 49), (25, 53), (25, 55), (25, 59), (25, 61), (25, 65), (25, 67), (25, 71), (25, 73), (25, 77), (25, 79), (25, 83), (25, 85), (25, 89), (25, 91), (25, 95), (25, 97), (25, 101), (25, 103), (25, 107), (25, 109), (25, 113), (25, 115), (25, 119), (25, 121), (25, 125), (25, 127), (25, 131), (25, 133), (25, 137), (25, 139), (25, 143), (37, 1), (37, 5), (37, 7), (37, 11), (37, 13), (37, 17), (37, 19), (37, 23), (37, 25), (37, 29), (37, 31), (37, 35), (37, 37), (37, 41), (37, 43), (37, 47), (37, 49), (37, 53), (37, 55), (37, 59), (37, 61), (37, 65), (37, 67), (37, 71), (37, 73), (37, 77), (37, 79), (37, 83), (37, 85), (37, 89), (37, 91), (37, 95), (37, 97), (37, 101), (37, 103), (37, 107), (37, 109), (37, 113), (37, 115), (37, 119), (37, 121), (37, 125), (37, 127), (37, 131), (37, 133), (37, 137), (37, 139), (37, 143), (49, 11), (49, 13), (49, 17), (49, 19), (49, 29), (49, 31), (49, 35), (49, 37), (49, 55), (49, 59), (49, 61), (49, 65), (49, 83), (49, 85), (49, 89), (49, 91), (49, 101), (49, 103), (49, 107), (49, 109), (49, 127), (49, 131), (49, 133), (49, 137), (61, 1), (61, 5), (61, 7), (61, 11), (61, 13), (61, 17), (61, 19), (61, 23), (61, 25), (61, 29), (61, 31), (61, 35), (61, 37), (61, 41), (61, 43), (61, 47), (61, 49), (61, 53), (61, 55), (61, 59), (61, 61), (61, 65), (61, 67), (61, 71), (61, 73), (61, 77), (61, 79), (61, 83), (61, 85), (61, 89), (61, 91), (61, 95), (61, 97), (61, 101), (61, 103), (61, 107), (61, 109), (61, 113), (61, 115), (61, 119), (61, 121), (61, 125), (61, 127), (61, 131), (61, 133), (61, 137), (61, 139), (61, 143), (73, 1), (73, 5), (73, 7), (73, 11), (73, 13), (73, 17), (73, 19), (73, 23), (73, 25), (73, 29), (73, 31), (73, 35), (73, 37), (73, 41), (73, 43), (73, 47), (73, 49), (73, 53), (73, 55), (73, 59), (73, 61), (73, 65), (73, 67), (73, 71), (73, 73), (73, 77), (73, 79), (73, 83), (73, 85), (73, 89), (73, 91), (73, 95), (73, 97), (73, 101), (73, 103), (73, 107), (73, 109), (73, 113), (73, 115), (73, 119), (73, 121), (73, 125), (73, 127), (73, 131), (73, 133), (73, 137), (73, 139), (73, 143), (85, 1), (85, 5), (85, 7), (85, 11), (85, 13), (85, 17), (85, 19), (85, 23), (85, 25), (85, 29), (85, 31), (85, 35), (85, 37), (85, 41), (85, 43), (85, 47), (85, 49), (85, 53), (85, 55), (85, 59), (85, 61), (85, 65), (85, 67), (85, 71), (85, 73), (85, 77), (85, 79), (85, 83), (85, 85), (85, 89), (85, 91), (85, 95), (85, 97), (85, 101), (85, 103), (85, 107), (85, 109), (85, 113), (85, 115), (85, 119), (85, 121), (85, 125), (85, 127), (85, 131), (85, 133), (85, 137), (85, 139), (85, 143), (97, 7), (97, 11), (97, 13), (97, 17), (97, 29), (97, 31), (97, 35), (97, 37), (97, 59), (97, 61), (97, 65), (97, 67), (97, 79), (97, 83), (97, 85), (97, 89), (97, 101), (97, 103), (97, 107), (97, 109), (97, 131), (97, 133), (97, 137), (97, 139), (109, 1), (109, 5), (109, 7), (109, 11), (109, 13), (109, 17), (109, 19), (109, 23), (109, 25), (109, 29), (109, 31), (109, 35), (109, 37), (109, 41), (109, 43), (109, 47), (109, 49), (109, 53), (109, 55), (109, 59), (109, 61), (109, 65), (109, 67), (109, 71), (109, 73), (109, 77), (109, 79), (109, 83), (109, 85), (109, 89), (109, 91), (109, 95), (109, 97), (109, 101), (109, 103), (109, 107), (109, 109), (109, 113), (109, 115), (109, 119), (109, 121), (109, 125), (109, 127), (109, 131), (109, 133), (109, 137), (109, 139), (109, 143), (121, 1), (121, 5), (121, 7), (121, 11), (121, 13), (121, 17), (121, 19), (121, 23), (121, 25), (121, 29), (121, 31), (121, 35), (121, 37), (121, 41), (121, 43), (121, 47), (121, 49), (121, 53), (121, 55), (121, 59), (121, 61), (121, 65), (121, 67), (121, 71), (121, 73), (121, 77), (121, 79), (121, 83), (121, 85), (121, 89), (121, 91), (121, 95), (121, 97), (121, 101), (121, 103), (121, 107), (121, 109), (121, 113), (121, 115), (121, 119), (121, 121), (121, 125), (121, 127), (121, 131), (121, 133), (121, 137), (121, 139), (121, 143), (133, 1), (133, 11), (133, 13), (133, 19), (133, 23), (133, 25), (133, 29), (133, 31), (133, 35), (133, 37), (133, 41), (133, 47), (133, 49), (133, 59), (133, 61), (133, 71), (133, 73), (133, 83), (133, 85), (133, 91), (133, 95), (133, 97), (133, 101), (133, 103), (133, 107), (133, 109), (133, 113), (133, 119), (133, 121), (133, 131), (133, 133), (133, 143)
Flip 3 coins 1/5184 = 0.00019290123456790122 (13, 13), (13, 35), (13, 41), (13, 43), (13, 49), (13, 71), (13, 85), (13, 107), (13, 113), (13, 115), (13, 121), (13, 143), (25, 5), (25, 19), (25, 29), (25, 31), (25, 41), (25, 55), (25, 65), (25, 67), (25, 77), (25, 91), (25, 101), (25, 103), (25, 113), (25, 127), (25, 137), (25, 139), (37, 1), (37, 5), (37, 7), (37, 11), (37, 13), (37, 17), (37, 19), (37, 23), (37, 25), (37, 29), (37, 31), (37, 35), (37, 37), (37, 41), (37, 43), (37, 47), (37, 49), (37, 53), (37, 55), (37, 59), (37, 61), (37, 65), (37, 67), (37, 71), (37, 73), (37, 77), (37, 79), (37, 83), (37, 85), (37, 89), (37, 91), (37, 95), (37, 97), (37, 101), (37, 103), (37, 107), (37, 109), (37, 113), (37, 115), (37, 119), (37, 121), (37, 125), (37, 127), (37, 131), (37, 133), (37, 137), (37, 139), (37, 143), (49, 19), (49, 29), (49, 55), (49, 65), (49, 91), (49, 101), (49, 127), (49, 137), (61, 5), (61, 19), (61, 25), (61, 35), (61, 41), (61, 55), (61, 61), (61, 71), (61, 77), (61, 91), (61, 97), (61, 107), (61, 113), (61, 127), (61, 133), (61, 143), (85, 5), (85, 7), (85, 13), (85, 35), (85, 41), (85, 43), (85, 49), (85, 71), (85, 77), (85, 79), (85, 85), (85, 107), (85, 113), (85, 115), (85, 121), (85, 143), (97, 29), (97, 31), (97, 65), (97, 67), (97, 101), (97, 103), (97, 137), (97, 139), (109, 1), (109, 5), (109, 7), (109, 11), (109, 13), (109, 17), (109, 19), (109, 23), (109, 25), (109, 29), (109, 31), (109, 35), (109, 37), (109, 41), (109, 43), (109, 47), (109, 49), (109, 53), (109, 55), (109, 59), (109, 61), (109, 65), (109, 67), (109, 71), (109, 73), (109, 77), (109, 79), (109, 83), (109, 85), (109, 89), (109, 91), (109, 95), (109, 97), (109, 101), (109, 103), (109, 107), (109, 109), (109, 113), (109, 115), (109, 119), (109, 121), (109, 125), (109, 127), (109, 131), (109, 133), (109, 137), (109, 139), (109, 143), (121, 5), (121, 7), (121, 19), (121, 29), (121, 41), (121, 43), (121, 55), (121, 65), (121, 77), (121, 79), (121, 91), (121, 101), (121, 113), (121, 115), (121, 127), (121, 137), (133, 19), (133, 25), (133, 35), (133, 41), (133, 61), (133, 71), (133, 91), (133, 97), (133, 107), (133, 113), (133, 133), (133, 143)
Flip 4 coins 1/5184 = 0.00019290123456790122 (37, 37), (37, 71), (37, 109), (37, 143), (109, 37), (109, 71), (109, 109), (109, 143)
Flip 5 coins 1/864 = 0.0011574074074074073 (13, 13), (13, 71), (13, 85), (13, 143), (133, 61), (133, 71), (133, 133), (133, 143)
Flip 6 coins 23/5184 = 0.004436728395061728 (61, 29), (61, 31), (61, 101), (61, 103), (85, 29), (85, 55), (85, 101), (85, 127)
Flip 7 coins 5/1152 = 0.004340277777777778 (61, 29), (61, 31), (61, 101), (61, 103), (85, 29), (85, 55), (85, 101), (85, 127)


To be continued...
Top
Piotr
Posted: Dec 27 2016, 12:07 PM


Unregistered









The problem is that 144 is very composite, so using LCG will make small cycles after doing the result mod divisor. For example all LCGs that use mod even either have all results odd or alternate odd and even. A good LCG can be done by using a mod of 2147483647 (prime) and multiplying by 16807 (no adding). It will make an uniformly random generator from 1 to 2147483646. Primes near powers of 2 are most convenient (3, 7, 17, 31, 127, 257, 8191, 65537, ...)
Top
« Next Oldest | General | Next Newest »
DealsFor.me - The best sales, coupons, and discounts for you

Topic Options



Hosted for free by zIFBoards* (Terms of Use: Updated 2/10/2010) | Powered by Invision Power Board v1.3 Final © 2003 IPS, Inc.
Page creation time: 0.0272 seconds · Archive