Why do changing the limits of RandomInteger not affect the result?

Question

Recently I noticed (on Windows 10, Mathematica 12.1) that changing the limits of RandomInteger, while keeping SeedRandom fixed, does not affect the result and produces a weird ~~non-random~~ non-pseudo random pattern. For example,

Table[SeedRandom[12345678]; RandomInteger[{1, i}], {i, 2, 50}]

returns

{2, 1, 4, 4, 4, 7, 7, 2, 2, 2, 2, 13, 13, 13, 13, 9, 9, 9, 9, 21, 21,21, 21, 21, 26, 26, 26, 26, 26, 26, 26, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18}

Naively I would expect that by changing the range of the random numbers would affect the result, as the chance to get a random integer in the range [1,4] is not the same as in [1,5] even with a fixed SeedRandom.

Moreover, this seems to be affected by the method used in the SeedRandom function:

CountDistinct[Table[SeedRandom[12345678, Method -> #];RandomInteger[{1, i}], {i, 2, 9999}]] & /@ {"Congruential","ExtendedCA", "Legacy", "MersenneTwister", "MKL", "Rule30CA"}

returns the number of distinct numbers

{35, 12, 25, 25, 7514, 25}

which implies that using the MKL method (only available for Intel processors) offers the largest variety.

However, RandomReal gives very different results:

CountDistinct[Table[SeedRandom[12345678, Method -> #];RandomReal[{1, i}], {i, 2, 9999}]] & /@ {"Congruential","ExtendedCA", "Legacy",MersenneTwister", "MKL", "Rule30CA"}
{9998, 9998, 9998, 9998, 9998, 9998}

So, what is going on here? Is there some sort of weird bug in RandomInteger or am I missing something?

Lukas Lang · Accepted Answer · 2021-07-16T13:45:14.633

I'm pretty sure this is expected. Consider the following implementation of RandomInteger:

randomBit[] := RandomInteger[{0, 1}]
rand[{min_, max_}] :=
 Module[
  {diff, n, res = Infinity},
  diff = max - min;
  n = Ceiling@Log2[max - min + 1];
  While[res > diff,
   res = FromDigits[Table[randomBit[], n], 2]
   ];
  min + res
  ]

Note that we use RandomInteger[{0,1}] internally to generate single random bits using Mathematica's RNG. The numbers are then generated by generating n random bits, where n is the smallest number of bits that can represent the required number of different outcomes (e.g. n==3 when we need 8 outcomes). We then convert these bits into a number. If the number is too large, we try again, otherwise we return it (this is called rejection sampling). E.g. if we want a number between 0 and 5, the following could happen:

We need 3 bits, since 2 bits can only represent 4 different outcomes
So we generate 3 bits, e.g. {1, 1, 0}, which is 6
This is too high, so we try again and get e.g. {1, 0, 0}, which is 4
This is ok, so we return it.

We can verify that this seems to be the internal mechanism used by RandomInteger:

SeedRandom[12345678];
res1 = RandomInteger[{3, 10}, 10]
(* {3, 8, 6, 5, 4, 4, 10, 7, 9, 3} *)
SeedRandom[12345678];
res2 = Table[rand[{3, 10}], 10]
(* {3, 8, 6, 5, 4, 4, 10, 7, 9, 3} *)
res1 == res2
(* True *)

We can now also explain what you are seeing: Let say we generate a number in the ranges {0,5}, {0,6} and {0,7}, all with the same seed. In all cases, we need 3 bits, so let's say we got {1, 0, 0} again, which is 4. Since in all 3 cases, 4 is a valid answer, we return it, resulting in the same result in all cases.

As for the other methods: It looks like both "Legacy" and "MKL" use different strategies for RandomInteger, where obviously the one used by "MKL" does not have this particular behavior.

Update

It seems we can reproduce the results from "MKL" using the following strategy:

randomBits[n_] := IntegerDigits[RandomInteger[{0, 2^31 - 1}], 2, 31][[;; n]]
rand[{min_, max_}] :=
 Module[
  {diff, n, res = Infinity},
  diff = max - min;
  n = Ceiling@Log2[max - min + 1];
  While[res > diff,
   res = FromDigits[randomBits[n], 2]
   ];
  min + res
  ]
SeedRandom[12345678, Method -> "MKL"];
res1 = RandomInteger[{3, 10}, 10]
(* {9, 6, 7, 9, 3, 7, 8, 6, 5, 7} *)
SeedRandom[12345678, Method -> "MKL"];
res2 = Table[rand[{3, 10}], 10]
(* {9, 6, 7, 9, 3, 7, 8, 6, 5, 7} *)
res1 == res2
(* True *)

The overall method is the same, we simply extract all n bits from a single 31 bit integer (no idea why not 32).

Crystal clear answer, many thanks. The behavior I noticed now makes perfect sense and I learned a lot about how RandomInteger works internally. I will wait a bit before accepting the answer, as is customary. — Hans Olo, Jul 15 '21 at 13:10

Why do changing the limits of RandomInteger not affect the result?

1 Answers1

Update