Exact Discrete pseudo-Random Sampling

[ttw]

There was an article (Daniel Lemire perhaps) recently on the inexactness of discrete sampling with pseudorandom number generators. (I'll use sample a single die toss.) He pointed out that a PRNG is usually based on a power of 2 but 6 is not such a power. THis leads to a sampling error in the the numbers from 1-6 are not sampled equally. For example, using an 8-bit generator (too short for real use but good for pedagogy) with an Linear Congruential Generator, X(n+1)=181*X(n)+181 mod 256, one has 42 samples of 1-6 and 4 more numbers leading to an error of 4/256. One suggested correction is to throw away numbers that are too long. Another (and probably good enough in practice) is to use a really large M (X=AX+b mod M is the recursion) like 64 or 128 or 256. There's still a bias but it's too small to matter.

I developed another method for my HP-50 dice game. Analogously to the Sophie-Germain-Safe Prime stuff, I use a prime modulus M=6*P+1, where P is another prime. The 6 is the size of the sample space. It could be replaced by whatever is useful in another problem. The point is that it's not hard (but it is expensive) to generate large M. One can use M=79=6*78+1 to see how the method works.

Next one has to find a multiplier A which is a primitive root of M. It's fairly easy on the HP50 or with a big number math package. This ensures the the LCG X=A*X mod M has cycle M-1 (=6*P). The output (as an integer) goes through P sets of each of 1 - 6.

To ensure small correlation, A/M as a continued fraction (a1,a2,a3...) should have each partial quotient Ai and the sum of partial quotients small. This takes a search. I did a short search: M=824633721547 and A=319872179753.

At this point I found that the HP50 has bugs with integers bigger than 2^40 (or something near there.) SORT does not work on large integers. FLOOR can fail; trying 'P/Q'
M
*
FLOOR
fails for large M.
However,
P/Q FXND SWAP
M SWAP
* / FLOOR
works.

If there's interest, I'll post my primitive root and quadratic residue and continued fraction codes later.