Weyl sequence overview

February 24th, 2016

Awhile back I ran across the stackoverflow question: What’s the origin of this GLSL rand() one-liner:

float rand(vec2 co){ return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453); }

to which I provided this guess of an answer. Writing that response reminded me of a closely related structure that doesn’t seem to be well known, the Weyl generator.

This post is intended to provide minimal background on this family and I’ll present some techniques based on them later.

The theory or this doesn’t really work

Weyl sequences (or generators) are named after Hermann Weyl for the equidistribution theorem and Weyl’s criterion SEE: Normal numbers and equidistributed sequences.

One possible way to generate a uniform sequence, where the members are named as follows:

\[\left\{ s_{0},s_{1},s_{2},s_{3},\cdots\right\}\]

is a polynomial mod 1 that contains at least one irrational coefficient other than the constant term. For degrees 1 and 2, where $ \alpha $ is an irrational constant, $ n \in \mathbb{N}_0 $ (positive integers including zero):

\[\begin{equation} \label{eq:rweyl1} s_{n}=\left(n\thinspace \alpha + \thinspace s_{0}\right)\bmod 1 \end{equation}\] \[\begin{equation} \label{eq:rweyl2} s_{n}=\left(n^2\thinspace \alpha + \thinspace s_{0}\right)\bmod 1 \end{equation}\]

both produce uniform sequences with $ s_{n} \in \left[0,1\right) $

The so far unspoken trick to the above is that it requires the operations to be performed in reals, as the properties are based on the infinite string of digits (in some base) of irrationals.

The real world without Reals

Working in finite math requires some tweaks, let’s start with floating point. We really want the low bits of the operations so let’s call our logical irrational $ k $, then our constant $ \alpha = k \bmod 1 $. We’ll call $ \eqref{eq:rweyl1} $ a “Weyl generator” and its formulation remains the same.

Equation $ \eqref{eq:rweyl1} $ can also be expressed as a recurrence relation (a.k.a. additive recurrence):

\[\begin{equation} \label{eq:fadditive} s_{n+1}=\left(s_{n}+\alpha\right)\bmod 1 \end{equation}\]

to state the obvious $ \eqref{eq:rweyl1} $ and $ \eqref{eq:fadditive} $ do not produce the same sequences in floating point. Also the modulo of $ \eqref{eq:fadditive} $ can be removed as it simply subtracts one if it has moved out of the interval (say if you have a fast select). We need to rework \eqref{eq:rweyl2} and we’ll call it a “nested Weyl”:

\[\begin{equation} \label{eq:fweyl2} s_{n}=n\left(\left(n\thinspace \alpha\right)\bmod 1\right)\bmod 1 \end{equation}\]

The generally better option is to convert to integers by scaling to some bit width $ b $, then $ \eqref{eq:rweyl1} $ and $ \eqref{eq:fadditive} $ become:

\[\begin{equation} \label{eq:iweyl1} s_{n}=\left(n\thinspace K + \thinspace s_{0}\right)\bmod 2^b \end{equation}\] \[\begin{equation} \label{eq:iadditive} s_{n+1}=\left(K + s_{n}\right)\bmod 2^b \end{equation}\]

restricting $ K $ to odd makes these one-to-one mappings (and bijections and permutations), so the sequence is periodic with period $ 2^b $. As such the additive constant $ s_0 $ isn’t needed. So the code for 32-bit might look like this:

static inline uint32_t weyl(uint32_t n)      { return K*n;   }  /* (5) */
static inline uint32_t weyl_next(uint32_t s) { return K+s;   }  /* (6) */

The constants are the secret sauce

The important part is the choice of constant(s) and that choice depends on how the sequence is going to be used. The only thing I’ll say here is that if you want full-period for integer, then the constant must be odd. For full period additive recurrence in floating-point we can follow the same reasoning:

\[\begin{equation} \label{eq:falpha} \alpha=\frac{i}{2^{s-1}} \end{equation}\]

where $ s $ is the significand precision for normal values and $ i $ is a odd integer on $\left[1,\thinspace2^{s-1}\right) $ (so $ s $ = 24 for singles).

In my opinion the main interest here is cheaply creating low-discrepancy sequences using the simple Weyl in both explicit and recurrence form. A rough overview for 1D and 2D can be found in this blog post.

Note that the comment:

The algorithm is symmetric about $ \alpha = 0.5 $

is due to the fact that $ \left(1-\alpha\right) $ produce the reversed sequence of $ \alpha $.

Since this has no figures or interactive demos, here’s a shadertoy. This breaks space into square cells and fills each with a fixed number of points using additive recurrence. This method will have defects along grid boundaries since each is independent.

Comments

math (32)