# Uniform PRNG for long integer structures

Good morning!

I don't actually know where to attribute this question (maybe it's better to publish it on StackOverflow), but it's more related to math theory than to actual realization.

Since the integer type is only 32-bit, I've built a LongInt type which is actually a vector of integer digits; the base of the digits is variable and may be set by a programmer.

So I want to generate long integers from (inclusive)

$A = a_0 * BASE^0 | a_1 * BASE^1 | a_2 * BASE^2 | ... | a_n * BASE^n | 0*BASE^{n+1} | ... | 0*BASE^m$

to (exclusive)

$B = b_0 * BASE^0 | b_1 * BASE^1 | b_2 * BASE^2 | ... | b_m * BASE^m$

$m>=n,\quad \exists j \in 1..m : a_j < b_j \quad and \quad \forall i:\quad j<i<=m: \quad a_i <= b_i$

( Digits are not actually multiplied by the corresponding power of BASE, I just written that for clarity. The actual structure is a_0 | a_1 | ... | a_n )

All I currently have is a plenty of uniform PRNGs for integers. We assume that for any 32-bit integer ranges, the distribution is uniform.

Can I somehow build a PRNG for my LongInt by making use of an integer PRNG - and retain uniformity? I've never seen any literature or solutions on the topic.

Idea. I may somehow generate digits one by one, starting with most significant ones. Let $C = B - 1$.

The algorithm is as follows (starting with most significant digits): 1. For every i such that $a_i = c_i$, just write the digit into the result. 2. Then start generating digits in $[a_i; b_i]$. 3. Once the generated digit is strictly less than to $c_i$, all less significant digits will be generated in the range $[0; BASE-1]$.

But this method seems not to guarantee any uniformity distribution (I'll conduct some tests and maybe publish the histograms).

Do you have any ideas how to solve this problem?

-

If you can sample uniformly from the real interval [0,1) then you can adjust your idea.

Sampling digit by digit uniformly does not work: If BASE=10, a = 1, and b = 13, then 75% of the numbers have first digit 0, and 25% have first digit 1.

75%: 01, 02, 03, 04, 05, 06, 07, 08, 09

25%: 10, 11, 12

However, if you use your method, 50% of generated numbers will have first digit 0 and 50% will have first digit 1.

However, this just tells you what distribution you need to sample from: the first digit is chosen from {0,1} with weights {75%, 25%}, and then the second digit is either chosen from {1,2,3,4,5,6,7,8,9} uniformly if the first digit was 0, otherwise from {0,1,2} uniformly.

So all you need to do is count how many numbers (as a percentage) have each first digit. Generally there are three cases: (1) first digit is smallest, so "a" matters, (2) first digit is largest, so "b" matters, (3) first digit is safely in the middle. (An actual implementation has to deal with the problem that (3) may not occur, and that both (1) and (2) can occur simultaneously, but these cases are easier to analyze.)

If you are ok with approximately uniform, then just use the difference (subtraction) of the leading digits to get an approximation of the ratio for use with your short PRNG.

If you need exactly uniform, then I think this idea just hides the problem in requiring a really good [0,1) PRNG. I suspect though that most uses are fine with the approximation (if you are using 32-bit PRNGs then the difference from uniform is less than 1/232).

If you need exactly uniform, then you can use the inefficient but unbiased and uniform method: generate all the digits uniformly at random from 0 to BASE. If the number lies in your range, return it. Otherwise, generate another. If you have a bignum-division, then you can use a modulo operation to speed this up in the same fashion as converting a 32-bit PRNG to a uniform distribution on {0,1,2,3,4,5,6,7,8,9}.

In more detail: Let C = B−A. Let n be the number of digits of C. Generate n random digits (each one uniformly from 0 to BASE−1 inclusive), and call this number D. If D is bigger than floor(BASEn/C)*C (so D is very large and near the top of the range), discard D and try again (or suffer some minor non-uniformity of size C/232). Otherwise, reduce D modulo C and add the result to A, that is return (D % C) + A.

For example, to sample exactly uniformly from {0,1,2,3,4,5,6,7,8,9} you cannot use a 32-bit PRNG and just take the result mod 10 since the resulting distribution is gives each of 0,1,2,3,4,5 weight 10.00000001% and 6,7,8,9 weight 9.999999986%. So what you do is ask the PRNG for a number D. If D ≤ 4294967290, then return D % 10, otherwise discard D and try again.

Again, for a lot of people these 0.00000001% deviations from uniform don't matter.

Why do we discard large D? Imagine we want to sample uniformly from {0,1,2,3,4}, but we only have a six-sided die. Our algorithm can be simple: roll the die and take the result mod 5. Unfortunately, we quickly see this is not uniform: We get 0 more than 33% of the time, while each of 1,2,3,4 occur less than 17% of the time. What went wrong?

We tried to use too much of the random number generator's range. The first 5 numbers are fine, but the next "5" numbers only had one number in it!

Imagine this algorithm instead: We roll the die. If the number is ≤ 5, then we return it mod 5, otherwise we roll again.

This returns each number exactly 20% of the time, but sometimes when we run the algorithm it takes a long time to return because it keeps re-rolling. The average number of rolls is 1.2, so its not too bad on average.

Ok, what if we have a 12-sided die? The correct algorithm is then: roll the die. If the number is ≤ 10, then return it mod 5, otherwise roll again. Again, this returns each of 0,1,2,3,4 exactly 20% of the time. The average number of rolls needed is still 1.2, so again not too bad on average.

We have the same problem generating our big random-numbers from 0 to C (exclusive). We need to leave out the "top-most" part of the random number generator, otherwise the generator will not be even. We can use any full multiple of C, but the "last" multiple of C, the one that is partial, we have to discard.

-
Jack, thanks for your answer! Am I understanding you right: I generate all digits from 0 to BASE-1 and then take the remainder modulo B - will this provide me exact uniformity? – wh1t3cat1k Apr 30 '11 at 9:17
More or less, yes. I added a more detailed explanation. – Jack Schmidt Apr 30 '11 at 9:42
I'm sorry I still can't quite understand what does $floor(BASE^n/C)*C$ actually mean. – wh1t3cat1k Apr 30 '11 at 9:50
Ok, I tried to explain. Basically, we need to divide up the PRNG's output evenly, but there will be a "leftover" part at the top that we have to discard. – Jack Schmidt Apr 30 '11 at 10:08
It appears that GMP uses the exact method: generate each digit uniformly, if the result is in the leftover piece, retry (up to a certain number of times), otherwise return the result mod C. The code is in mpz/urandomm.c – Jack Schmidt Apr 30 '11 at 10:22