I have 32 bytes of random data.
I want to generate random numbers within variable ranges between 0-9 and 0-100.
If I used an arbitrary precision arithmetic (bignum) library, and treated the 32 bytes as a big number, I could simply do:
random = random_source % range;
random_source = random_source / range;
as often as I liked (with different ranges) until the product of the ranges nears 2^256.
Is there a way of doing this using only (fixed-size) integer arithmetic?
Certainly you can do this by doing base 256 long division (or push up multiplication). It is just like the long division you learnt in primary school, but with bytes instead of digits. It involves doing a cascade of divides and remainders for each byte in turn. Note that you also need to be aware how you are consuming the big number, and that as you consume it and it becomes smaller, there is an increasing bias against the larger values in the range. Eg if you only have 110 left, and you asked for a rnd(100), the values 0-9 would be 10% more likely than 10-99 each.
But, you don't really need the bignum techniques for this, you can use ideas from arithmetic encoding compression, where you build up the single number without actually ever dealing with the whole thing.
If you start by reading 4 bytes to an unsigned uint_32 buffer, it has a range 0..4294967295 , a non-inclusive max of 4294967296. I will refer to this synthesised value as the "carry forward", and this exclusive max value is also important to record.
[For simplicity, you might start with reading 3 bytes to your buffer, generating a max of 16M. This avoids ever having to deal with the 4G value that can't be held in a 32 bit integer.]
There are 2 ways to use this, both with accuracy implications:
Stream down:
Do your modulo range. The modulo is your random answer. The division result is your new carry forward and has a smaller range.
Say you want 0..99, so you modulo by 100, your upper part has a range max 42949672 (4294967296/100) which you carry forward for the next random request
We can't feed another byte in yet...
Say you now want 0..9, so you modulo by 10, and now your upper part has a range 0..4294967 (42949672/100)
As max is less than 16M, we can now bring in the next byte. Multiply it by the current max 4294967 and add it to the carry forward. The max is also multiplied by 256 -> 1099511552
This method has a slight bias towards small values, as 1 in the "next max" times, the available range of values will not be the full range, because the last value is truncated, but by choosing to maintain 3-4 good bytes in max, that bias is minimised. It will only occur at max 1 in 16million times.
The computational cost of this algorithm is the div by the random range of both carry forward and max, and then the multiply each time you feed in a new byte. I assume the compiler will optimise the modulo
Stream up:
Say you want 0..99
Divide your max by range, to get the nextmax, and divide carryforward by nextmax. Now, your random number is in the division result, and the remainder forms the value you carry forward to get the next random.
When nextmax becomes less than 16M, simply multiply both nextmax and your carry forward by 256 and add in the next byte.
The downside if this method is that depending on the division used to generate nextmax, the top value result (i.e. 99 or 9) is heavily biased against, OR sometimes you will generate the over-value (100) - this depends whether you round up or down doing the first division.
The computational cost here is again 2 divides, presuming the compiler optimiser blends div and mod operations. The multiply by 256 is fast.
In both cases you could choose to say that if the input carry forward value is in this "high bias range" then you will perform a different technique. You could even oscillate between the techniques - use the second in preference, but if it generates the over-value, then use the first technique, though on its own the likelihood is that both techniques will bias for similar input random streams when the carry forward value is near max. This bias can be reduced by making the second method generate -1 as the out-of-range, but each of these fixes adds an extra multiply step.
Note that in arithmetic encoding this overflow zone is effectively discarded as each symbol is extracted. It is guaranteed during decoding that those edge values won't happen, and this contributes to the slight suboptimal compression.