Unique Random Iterator

On quite some occasions I have been needing to iterate a range in a random order, without visiting the same position twice.

1. Naive Attempt

I think that the usual approach involves creating an array with all positions, and then shuffling it. It ensures that each position is only iterated once, and you can easily rely on your standard library which should provide you with a good shuffle. 

// initialize values in range - O(n) memory and complexity
int range = 10000;
std::vector v;
v.reserve(range);
for (auto i = range; i >= 0; --i) {
    v.push_back(i);
}

// shuffle all positions - at least O(n) complexity
std::shuffle(v.begin(), v.end(), g);

// example usage
for (auto &position : v) {
    if (is_suitable(position)) {
        return position;
    }
}

However, this involves unnecessarily allocating a lot of memory for big ranges, especially if you expect to only need a small portion of the range.

2. Prime Addition

Iterating over a range in a random fashion can be done adding an arbitrary prime bigger than the array size to the current position and using modulo on the result. This is a known trick, however finding a good prime is not easy. 

A is the array
S is the size of the array
Generate a random prime number (P), such that P > S
Q = P % S
for i = 1 to S
    process A[Q]
    Q = (Q + P) % S

3. LCG Based Random Iterator

I have been implementing a LCG wrapped as an interface to iterate over a specified integer range. We make use of its unique properties respecting the Hull-Dobell Theorem (Hull, T. E. and Dobell, A. R., 1962) to find a suitable value for our modulo value. Note that due to its nature, the generator will alternate between odd and even numbers. If we can accept the inferior results of a LCG, we can reduce the complexity and completely avoid the allocation of an array. 

#include <cstdint>

inline constexpr uint64_t next_pow2_64(uint64_t i) {
    --i;
    i |= i >> 1;
    i |= i >> 2;
    i |= i >> 4;
    i |= i >> 8;
    i |= i >> 16;
    i |= i >> 32;
    return ++i;
}

/*
 * LCG based random iterator implementation
 * Requirements (Hull-Dobell Theorem):
 * 1) mod and incr are relatively prime
 * 2) mult-1 is divisible by all prime factors of range
 * 3) mult-1 is divisible by 4 if mod is divisible by 4
 * We round up to the next power of two for mod to guarantee 2) and 3).
 * Outputs bigger than range (<50%) are discarded.
 */
class UniqueRandomIterator {
private:
    uint64_t state, count;
    const uint64_t mult, incr, mod_mask, range, offset;

public:
    inline UniqueRandomIterator(
        uint64_t min,
        uint64_t max,
        uint64_t seed = 0xC789F4F722444853,
        uint64_t stream = 0x5084B57B66D07A23,
        uint64_t mult = 0xC360248E09A02DF5
    )
    : mult(mult), incr((stream << 1) | 1)
    , range(max - min + 1)
    , offset(min), count(0)
    , mod_mask(next_pow2_64(max - min + 1) - 1) {
        state = seed % range;
    }

    inline uint64_t next() {
        do {
            state = (state * mult + incr) & mod_mask;
        } while (state >= range);
        if (++count > range) [[unlikely]] {
            count = 1;
        }
        return state + offset;
    }

    inline uint64_t current() const {
        return state + offset;
    }

    inline uint64_t left() const {
        return range - count;
    }

    inline bool next_available() const {
        return left() > 0;
    }
};

We have been avoiding modulo by using a power of two for our LCG modulo value, which allows us to use a bit-mask instead of an expensive division. Luckily this fits well with our requirements. By using a base offset, we support arbitrary integer ranges. Values that fall out of our power of two modulo range are discarded, and are less than 50% even in the worst case by definition.

3.1. Example Usage

You only need to specify your range, and an optional seed as third parameter. In this example we just print out our positions line by line. 

#include <cstdio>

int main() {
    UniqueRandomIterator it(1, 10);
    while (it.next_available()) {
        printf("%i\n", it.next());
    }
    return 0;
}

3.2. Example Output

While the approach shines for huge ranges that are not processed in whole, for showcase I have set the range from one to ten. 

5
9
3
2
4
7
6
1
8
10

3.3. Generated Assembly

This is using x86-64 GCC 13.2 using -O3. You can see that the resulting assembly is very compact and basic. 

.LC0:
        .string "%i\n"
main:
        push    r13
        mov     r13d, 10
        push    r12
        movabs  r12, -4368451446083932683
        push    rbp
        movabs  rbp, -6842820550568381369
        push    rbx
        mov     ebx, 9
        sub     rsp, 8
.L2:
        imul    rbx, r12
        add     rbx, rbp
        and     ebx, 15
        cmp     rbx, 9
        ja      .L2
        xor     eax, eax
        lea     rsi, [rbx+1]
        mov     edi, OFFSET FLAT:.LC0
        call    printf
        sub     r13, 1
        jne     .L2
        add     rsp, 8
        xor     eax, eax
        pop     rbx
        pop     rbp
        pop     r12
        pop     r13
        ret

4. Extension for 2D Ranges

By wrapping our iterator, we can split the range into a two dimensional area. This is useful, if for example, you need to iterate over a set of cells or tiles in a random fashion. 

struct UniqueRandomIterator2DPos {
    uint64_t x, y;
};

class UniqueRandomIterator2D : public UniqueRandomIterator {
private:
    const uint64_t offset_x;
    const uint64_t size_x;

public:
    UniqueRandomIterator2D(
        uint64_t x1,
        uint64_t x2,
        uint64_t y1,
        uint64_t y2,
        uint64_t seed,
        uint64_t stream = 0x5084B57B66D07A23
    )
    : UniqueRandomIterator(
        (x2 - x1 + 1) * y1,
        (x2 - x1 + 1) * (y2 + 1) - 1,
        seed, stream)
    , offset_x(x1), size_x(x2 - x1 + 1) {}

    auto to_pos(const uint64_t state) const {
        return UniqueRandomIterator2DPos {
            .x = (state % size_x) + offset_x,
            .y = (state / size_x),
        };
    }

    auto current_pos() const {
        return to_pos(current());
    }

    auto next_pos() {
        return to_pos(next());
    }
};

random_iterator_2d.gif

Figure 1: Visualization

5. Extension for 3D Ranges

For completion, the same can be done for three dimensions. Feel free to write your own n-dimensional template implementation. 

struct UniqueRandomIterator3DPos {
    uint64_t x, y, z;
};

class UniqueRandomIterator3D : public UniqueRandomIterator {
private:
    const uint64_t offset_x, offset_y;
    const uint64_t size_x, size_y;

public:
    UniqueRandomIterator3D(
        uint64_t x1,
        uint64_t x2,
        uint64_t y1,
        uint64_t y2,
        uint64_t z1,
        uint64_t z2,
        uint64_t seed,
        uint64_t stream = 0x5084B57B66D07A23
    )
    : UniqueRandomIterator(
        (x2 - x1 + 1) * (y2 - y1 + 1) * z1,
        (x2 - x1 + 1) * (y2 - y1 + 1) * (z2 + 1) - 1,
        seed, stream)
    , offset_x(x1), offset_y(y1)
    , size_x(x2 - x1 + 1), size_y(y2 - y1 + 1) {}

    auto to_pos(uint64_t state) const {
        return UniqueRandomIterator3DPos {
            .x = ((state % size_x) + offset_x),
            .y = ((state / size_x) % size_y) + offset_y,
            .z = ((state / size_x) / size_y),
        };
    }

    auto current_pos() const {
        return to_pos(current());
    }

    auto next_pos() {
        return to_pos(next());
    }
};