Quite a while ago, while discussing Monte Carlo Integration with my students, the topic of choosing sample locations came up, and we discussed low-discrepancy sequences (a.k.a. quasi-random sequences). In a low-discrepancy sequence, values generated look kind of uniform-random, but avoids clumping. A closer examination reveal that they are suspiciously well-spaced. That’s what we want in Monte Carlo integration.

But how do we generate such sequences? Well, there are many ways to do so. Some more amusing than other, some more structured than others. One of the early example, Halton sequences (c. 1964) is particularly well behaved: it generates 0, 0.5, then 0.25 and 0.75, then 0.125, 0.375, 0.625, and 0.875, etc. It does so with a rather simple binary trick.

Halton observed that if you count in binary from 0 to infinity, reverse the bits, and made it fractional, you’ll end up filling the interval [0,1):

n | binary | reversed | fraction |

0 | 0 | 0 | 0 |

1 | 1 | .1 | 0.5 |

2 | 10 | .01 | 0.25 |

3 | 11 | .11 | 0.75 |

4 | 100 | .001 | 0.125 |

5 | 101 | .101 | 0.625 |

6 | 110 | .011 | 0.375 |

7 | 111 | .111 | 0.875 |

Or, if we look at iterations graphically (up to 7 bits):

*

* *

The first element is then the bit-reversing function. It is a bit different from the usual version, because bit 0 doesn’t end-up in bit 31, but rather in a position that depends on the magnitude of the number. Well, fortunately, there’s an easy way to compute this (albeit not that efficient):

T p_halton_reverse(T x) { T t=0; while (x) { t<<=1; t|=(x&1); x>>=1; } return t; }

The next part is to compute the next power-of-two. There are quite many ways to do so, but in fact, we can maintain it incrementally as we generate the counting sequence. The limit is initially 1. When the counter reaches it, we double it. This entails a check every iteration, but it’s likely less expensive (and, anyway, a lot simpler) than recomputing the next power of two from scratch every time.

So the method is then:

- Next power of two is initially 1, counter is 0.
- Increment counter, adjust the next power of two.
- Output counter scaled by the next power of two.
- Goto 2, unless tired.

*

* *

Now, to make it fun and practical, it should be implemented as a range, or iterator, so that we can write something like:

for (auto t : halton_sequence(128)) std::cout << t << std::endl;

Fortunately, C++11-style range-based loops doesn’t require much from an object: `begin()` and `end()` functions that return an iterator object that knows how to to `!=`, `operator *()` (indirection) and `operator++`.

The “host object” must know how to do `begin()` and `end()`. It must contain the two functions:

const_iterator begin() const { return const_iterator(0,max_val); } const_iterator end() const { return const_iterator(max_val,max_val); }

Where `const_iterator(start_value,end_value)` is the iterator object itself. It doesn’t do much: it takes a `start_value` from which it will start to count, and an `end_value` that limits the count. In my implementation, the iterator is equivalent to `for (T i=start_value; i<end_value;i++)`. It may not be as general as possible, but having `<` rather than `<=` deals with a number of problems such as overflow and repeated values.

The iterator itself is very simple, as it contains only three important functions (or, more precisely, operators):

bool operator!=(const const_iterator & other) { return current_val!=other.current_val; } F operator*() const { return p_halton_reverse(current_val)/(F)next; } const const_iterator & operator++() { if (++current_val==next) next<<=1; return *this; }

The `operator!=` is needed so that the loop can compare to `end()`. The indirection operator is what will give you the value through the iterator. Finally, `operator++` advances the iterator to the next value.

*

* *

To make things fun, or at least practical, we could make this template code. The template for the host object should have two arguments, one for the type of integer used for the counter, and one for the float type used for the generated sequence. We could, as the STL often does, hide a couple of special cases behind typedef/using aliases. In any case, the class shouldn’t be hard to use:

for (auto t : halton_sequence<int,double>(16)) std::cout << t << std::endl;

*

* *

The complete PoC:

#include <iostream> #include <bitset> #include <cstdint> #include <limits> template <typename T=int, typename F=float> class halton_sequence { protected: const T max_val; static T p_halton_reverse(T x) { T t=0; while (x) { t<<=1; t|=(x&1); x>>=1; } return t; } static T p_next(T x) { T i=1; while (i<x) i<<=1; return i; } public: class const_iterator { protected: const T max_val; T current_val, next; public: bool operator!=(const const_iterator & other) { return current_val!=other.current_val; } F operator*() const { return p_halton_reverse(current_val)/(F)next; } const const_iterator & operator++() { if (++current_val==next) next<<=1; return *this; } const_iterator(const T & s, const T & m) : max_val(m), current_val(s), next(p_next(s)) {} const_iterator(const T & other) : max_val(other.max_val), current_val(other.current_val), next(other.next) {} ~const_iterator()=default; }; const_iterator begin() const { return const_iterator(0,max_val); } const_iterator end() const { return const_iterator(max_val,max_val); } halton_sequence(const T & max_=std::numeric_limits<T>::max()) : max_val(max_) {} ~halton_sequence()=default; }; int main() { for (auto t : halton_sequence<int,double>(16)) std::cout << t << std::endl; return 0; };