How does this algorithm for the Van der Corput sequence work?

Question

For any natural number $n$ write its binary expansion as $n = \sum_{i=0}^{k(n)} n_i 2^i$. Then the $n$th entry of the binary Van der Corput sequence is defined to be the dyadic rational

$$V(n) = \sum_{i=0} ^{k(n)} n_i 2^{-i-1}.$$

In words, $V(n)$ are the binary integers, reflected across the decimal point. The first few binary numbers are $$1,10,11,100,101,110,111,1000,1001,1010,1011,\dots$$ and correspondingly the first few entries of $V(n)$ are (in binary) $$ 0.1, 0.01, 0.11, 0.001, 0.101, 0.011, 0.111,0.0001,0.1001,0.0101,0.1101,\dots$$ or equivalently $$ \frac12, \frac14, \frac34,\frac18,\frac58,\frac38,\frac78,\frac1{16},\frac7{16},\frac5{16},\frac{13}{16},\dots$$

Apparently, this algorithm from the Wikipedia page for Halton sequences (a generalisation of $V$) computes the sequence $V(n)$:

def halton_sequence(b):
    """Generator function for Halton sequence."""
    n, d = 0, 1
    while True:
        x = d - n
        if x == 1:
            n = 1
            d *= b
        else:
            y = d // b
            while x <= y:
                y //= b
            n = (b + 1) * y - x
        yield n / d

I'm trying to understand this algorithm. Clearly, $n$ and $d$ denote the numerator and denominator, and $x=d-n>0$.

If $x=1$ i.e. $n/d = (b^{k}-1)/b^{k}$ then we need to increment the power of the denominator to get $1/b^{k+1}$; this is clear.

In the remaining case, we leave the denominator as is and focus on the next numerator. The number $y$ is the largest power of $b$ that is strictly smaller than $x=d-n$. Then, the next numerator is claimed to be $ n_{new} = (b+1) y - x = (b+1)y - d + n.$ In summary my pseudocode for the algorithm is:

1. given $(n, d)$ (initialised with $n=0,d=1$)
2. if $n = d-1$ then return $(1, db)$ else
3. define $y = \max_{k \ge 0} \{ b^k : n + b^k < d \}$
4. return $( (b+1)y - d + n, d)$

Q: How do I link this back to the original definition of $V(n)$?

I have checked the paper (Berblinger, Michael; Schlier, Christoph (1991). "Monte Carlo integration with quasi-random numbers: some experience". Computer Physics Communications. 66 (2–3): 157–166. doi:10.1016/0010-4655(91)90064-R. ISSN 0010-4655.) where this algorithm originates but did not find any clarifying remarks.

From steps 3 and 4, it seems the quantity $y-d+n \in [-d,0)$ could be useful. I've computed the first few iterations of $(n,d;y-d+n)$ for $b=2$. In the case when $n=d-1$ I fill the $y-d+n$ value with the symbol $X$ instead: $$ \text{init with }(0,1) \to (1,2; X) \to (1,4;X) \to (3,4;-1) \to (1,8;X) \to (5,8;-3) \to (3,8;-1) \to (7,8;-1) \to (1,16;X) $$ Its not so easy to search for this sequence on Math.SE; one other post on this sequence is Why is the formula for generating Van der Corput sequences called an Inverse Radical Function?

Answer 1

The algorithm works as follows:

1. We always have $d$ a power of $b$, $d=b^D$, say, and $0\le n< b^D$. So, write $n$ as a sequence of $D$ base-$b$ digits.
2. If $n$'s digits all equal $b-1$, $n=(b-1)\cdots (b-1)$, then we have exhausted all possible digit sequences of length $D$ and need to increase $D$. Then, $n=b^D-1$. Since this is equivalent to $x=d-n=b^D-(b^D-1)=1$, the test if x == 1 checks for this.
3. In this case, we multiply $d$ by $b$, meaning that we increase $D$ by 1, and set $n$ to $0\cdots01$ to start making digit sequences of length $D+1$.
4. Otherwise, we need to increment $n$. To do this, we need to find the leftmost digit of $n$ which is not $b-1$.
5. This is done by the loop on the variable $y$, which is a power of $b$, $b^i$, say, corresponding to a digit position in $n$. Since $x=d-n$, the test x <= y is equivalent to testing whether $d\le n + y$, i.e., whether $n\ge b^D - b^i$. Since $b^D-b^i$ has the base-$b$ expansion $(b-1)\cdots(b-1)0\cdots 0$, this will be true just when the digit positions in $n$ which correspond to $y$ or are to the left of $y$ all contain the digit $b-1$.
6. When we exit the while loop, we have found a digit which is not $b-1$. We must increment the digit, which we can do by adding $y=b^i$ to $n$, and clear the digits $b-1$ to the left of the digit, which we can do by adding $by=b^{i+1}$ to $n$. Clearing the digits $b-1$ will cause an overflow carry out of the leftmost digit of $n$, which we can remove by subtracting $d=b^D$ from $n$.
7. Therefore, we need to add $(b+1)y-d$ to $n$. Since $x=d-n$, this can be done by the assignment statement n = (b + 1) * y - x.
8. Finally, in both cases, we yield the generated value.