Daniel Lemire's blog

, 27 min read

Visiting all values in an array exactly once in “random order”

36 thoughts on “Visiting all values in an array exactly once in “random order””

  1. Federico says:

    For a “more random” solution, why not using exponentiation? x \mapsto x^h is injective modulo every prime p for which gcd(p-1, h) = 1. For a small enough h, for instance h = 3, it is reasonably fast to compute.

    1. Daniel Lemire says:

      I don’t want to assume that the size of the interval is a prime. You can, of course, patch things up, but I also want the code to be simple and easy to predict.

  2. Federico says:

    Another neat trick in the same direction is the following: for each r>2, the map k -> (5^k mod 2^r) runs exactly once through all numbers whose binary expansion ends in 01 before looping. So you can visit all values in an array of size 2^(r-2) with the function ((5^k mod 2^r) – 1) / 4. It should be reasonably fast, because it’s just a multiplication by 5 and a few shifts.

    1. Daniel Lemire says:

      I don’t want to assume that the interval is a power of two.

      1. Marc Reynolds says:

        If you have range of ‘n’ with a configurable power-of-two permutation then you can generate on ceil(log_2(n)) and reject elements >= n. So worst case rejection is ~1/2 or ~2 elements per draw on average. If I’m doing the math right then average rejection rate (assuming ‘n’ is uniform) is ~.18. Of course only worth considering if you need to speed up the sampling of elements.

        1. Daniel Lemire says:

          Yes, you are right… with the possible upside that you can get well known statistical properties… but with the downsides that…

          1) you cannot easily randomly access the indexes… e.g, in my implementation, I can ask what was the kth index,

          2) in your proposal, the running time of accessing the next integer is O(n), not constant time…

          3) my implementation is branchless… so no branch misprediction…

          If performance is a priority, I think that my approach is going to be much faster in realistic scenarios, albeit, we need to benchmark it.

  3. Federico says:

    Also, in your code, why don’t you update your index with answer = (old_answer + prime) % maxrange, instead of doing the multiplication? Or do you need to have non-sequential access to your array as well? If so, then my tricks with the multiplicative order don’t work. 🙂

    1. Daniel Lemire says:

      Yes, thanks. You don’t need the multiplication, you are right. But it is probably the least of your concerns here. The modulo reduction is probably want you want to optimize away first.

      I did point out in my post that this code can be further optimized.

  4. Julian Hyde says:

    Even though it has limitations, I’m very fond of the shuffle algorithm that appears in Knuth. It is essentially selection sort using a random comparison: for each i, choose a random element at position i or higher, emit it, swap it with the element at position i.

    For this problem, the shuffled array is not the goal, but is merely a by-product; we need to move elements out of the way so that we do not consider them more than once.

    So, we can implement the algorithm with no extra storage, and better randomness than your algorithm, if: (a) we allow the array to be permuted while the algorithm is running, (b) we have a random number generator that can run backwards, (c) 2x running time is acceptable. The idea is to reverse the permutation at the end of the algorithm, running the random number generator backwards to put each element back into its original position.

    Of course, this is a significant extra cost over your algorithm. But acceptable, I think, for applications that require good randomness.

    I don’t know whether there exist “good” random number generators that are reversible, but I assume that there are, since a random number generator that loses information is not a very good random number generator.

    1. Daniel Lemire says:

      Good point.

      This suggests a follow-up blog post.

  5. Vadim Kantorov says:

    Also a few weeks ago a post on the same problem appeared: http://fabiensanglard.net/fizzlefade/index.php mentioning Linear Feedback Shift Registers

    1. Daniel Lemire says:

      Thanks. I have not read in detail the post you offer, but I suspect it is a different problem.

  6. Simon Hardy-Francis says:

    Some years ago I was experimenting with one of the murmur hash algorithms and discovered that commenting out part of the algorithm resulted in randomly visiting each array element once if the array is a power of two large. I’ve used it a couple of times over the years. Very handy.

    1. Daniel Lemire says:

      The bit mixing function executed at the end of murmur hash is invertible so it is definitively the case that it will visit all values exactly once.

      1. llogiq says:

        This is actually not the case. The mixing function is:

        int mixing(int h, int len) {
        h ^= len;
        h ^= h >> 16;
        h *= 0x85ebca6b;
        h ^= h >> 13;
        h *= 0xc2b2ae35;
        h ^= h >> 16;
        return h % len;
        }

        This will *not* generate all values exactly once, as a simple experiment with 16 values shows: 12 8 7 4 1 5 6 12 15 10 13 8 0 10 13 11

        1. Marc Reynolds says:

          If you remove the xor len and mod len then the remaining part is a bijection (or invertible/permutation) function. To make it a sequence of all integers you need to feed it the output of some full-period sequence. Just 0,1,2…. will result in pretty respectable results.

          1. Daniel Lemire says:

            Yes. Marc is correct. What I refer to as “the bit mixing function” is equivalent to function you use (at least for 32-bit integers), without the “len” parameter. You will find it in widespread use.

            There are other possibilities, see section 7 of this paper : https://arxiv.org/pdf/1609.09840.pdf

  7. Sagar says:

    A Feistel Network? http://antirez.com/news/113

    1. Daniel Lemire says:

      Again, as other proposals above, it assumes that the size of the interval is a power of two.

      1. Sagar says:

        I haven’t looked into it, but the author of that blog says it is possible to generalize the Feistel network to any radix (see the comment by antirez in the link).

        1. Daniel Lemire says:

          the author of that blog says it is possible to generalize the Feistel network to any radix

          Well, we do know one way to do it, it has been described in the comments here, and it involves branching and O(n) complexity for the “next” calls.

          It is possible that he has something more clever in mind… but if it is described in the blog post in question, then I missed it.

  8. KWillets says:

    If c = (x – ax)%m then x maps to itself. However I believe sufficient conditions are given in theorem A here: http://www.math.cornell.edu/~mec/Winter2009/Luo/Linear%20Congruential%20Generator/linear%20congruential%20gen1.html , or you can stay with c=0 or a=0 (the other coefficient being relatively prime).

    1. Daniel Lemire says:

      If c = (x – ax)%m then x maps to itself.

      My proposed technique is not a recurrence formula. So while what you write is correct (albeit I think I use ‘b’ and not ‘c’), it does not matter.

      1. Kendall Willets says:

        You are correct, I misremembered it as iterating x = bx+c. Iteration may be another way to do it then.

      2. Marc Reynolds says:

        Well it’s additive recurrence…just in closed-form for the ‘i’th element. >;)

    2. Daniel Lemire says:

      I have added a benchmark against LCG: https://lemire.me/blog/2017/09/26/benchmarking-algorithms-to-visit-all-values-in-an-array-in-random-order/

  9. David says:

    Any solution that works for all powers of two trivially works for any size. Choose the smallest power of two greater than or equal to the size you want. If the algorithm gives an index too big just skip that one and go on to the next “random” index. This could possibly be faster than a solution that natively handles arbitrary sizes, because modulo instruction is quite slow on most hardware.

    1. Daniel Lemire says:

      As pointed out earlier, you do not need the modulo instruction.

      1. Daniel Lemire says:

        I have added benchmarks in a separate post: https://lemire.me/blog/2017/09/26/benchmarking-algorithms-to-visit-all-values-in-an-array-in-random-order/

  10. Devin says:

    You can also use an RNG to choose a permutation-index [0, n!). Each index represents a specific ordering of all of the elements. For example, 0 might represent the (original) order 0, 1, …, n – 2, n – 1; 1 might represent the order 0, 1, …, n – 1, n – 2; n!-1 might represent the order n – 1, n – 2, …, 1, 0; etc.

    The interesting part is defining a specific algorithm that satisfies the above properties in an efficient way. It’s similar to algorithms that generate all possible permutations… but in this case, we want to jump to a specific one. I’ve coded something like this up before.

    One such algorithm could be related to an efficient solution to Project Euler: https://projecteuler.net/index.php?section=problems&id=024. https://r.prevos.net/euler-problem-24/

    1. Daniel Lemire says:

      I think you refer to Lehmer codes.

      Did you post your code?

      1. Devin says:

        Ah, I think Lehmer code is on point. Although I’m not sure I knew that it was called that. I don’t have the specific code I wrote… think that was about 10 years ago. I’ve whipped something up though that is roughly equivalent (although not specifically Lehmer) https://gist.github.com/devinrsmith/0320756f5fc9b38fbd23c464dd516de9

        This is very similar in spirit the solution that Julian Hyde proposed.

  11. Petr says:

    Hi Daniel,
    I tried you C implementation of shuffle() / shuffleReallyFast()
    from here :
    https://github.com/lemire/Code-used-on-Daniel-Lemire-s-blog/tree/master/2017/09/25

    and result of shuffling of array initialized with items
    0, 1, 2, 3, .. 99
    is this:
    20, 71, 22, 73, 24, 75, 26, 77, 28, 79, ..

    Is this expected behaviour ? Output is just two interleaved growing
    sequences with stride of 2 (result of shuffleLCG() looks reasonably random).

    Thanks !

    1. Daniel Lemire says:

      The C code is not randomized, see the remark in the first line of visit.c:

      https://github.com/lemire/Code-used-on-Daniel-Lemire-s-blog/blob/master/2017/09/25/visit.c#L1

      I only wrote the C code to benchmark speed.

      In any case, as I make clear in my blog post, this fast approach is crude and will not fool anyone who is looking for true randomness.

  12. Petr says:

    Thanks a lot for explanation Daniel !

  13. Duke says:

    Great stuff! This really helped me a lot. I was looking for a randomizing algorithm that functions as a “Stateless Deck Shuffler for Very Casual Use” that:

    Must return each unique value within a set exactly once
    Must return the
    values in a pseudorandom order based on a seed
    Must allow access the values in an arbitrary
    order without generating or storing prior or available values first
    Does not need a hair of cryptographic security

    My particular use case is to randomize the order of minor award items given to a player in a game. The player is not likely to notice or care that the items are metered out by a weak randomization algorithm as long as the randomization is good enough to give a thin illusion of randomization.

    Using your article as a base, I was able to add some more “randomness” to my implementation by using a transposition and shifting function multiple times on the index before using it. The number of transpositions/shifts used can be tailored based on desired “randomness” and calculation time available.

    The transposition function I used was a simple de-interleave with the second half of the values reversed. This function could also be tailored to needs.

    It does takes a bit of finesse and a strong eye for patterns to tease out artifacts that can result from bad choices for the shift parameters.

    const int Prime = 53;
    const int Range = 100;
    static int Rando(int index, int seed)
    {
    var transformedIndex = Transpose(index, Range);
    transformedIndex = Transpose(Shift(transformedIndex, 39, Range), Range);
    transformedIndex = Transpose(Shift(transformedIndex, 7, Range), Range);
    transformedIndex = Transpose(Shift(transformedIndex, 68, Range), Range);
    transformedIndex = Transpose(Shift(transformedIndex, 27, Range), Range);

    return ((transformedIndex * Prime + seed) % Range);
    }

    static int Transpose(int i, int range)
    {
    var odd = (i & 1) == 1;
    var i2 = i >> 1;
    if(odd)
    {
    i2 = range - i2 - 1;
    }
    return i2;
    }

    static int Shift(int i, int shift, int range)
    {
    return (i + shift) % range;
    }