April 18th, 2016 Overview This public domain single header file library (Sobol.h) implements basic functions for sampling Sobol sequences. It implements two types of generators: Progressive, which is the standard sampling method. Specifically a generator can be queried for any number of samples. The library supports 1-3 dimensional generators. Stratified generators. These generate a sequence one dimension lower, require specifying the number of samples which will be generated and the remaining dimension is equispaced over the interval. Documentation The library consists of a data table and mostly inline functions. Any functions which are not defined as inline are noted below. As such one file must define a macro (#define SOBOL_IMPLEMENTATION) prior to including the header. Data structures Progressive generators contains an implied position element $ i $ followed by a state element per dimension. typedef struct { uint32_t i, d0; } sobol_1d_t; typedef struct { uint32_t i, d0,d1; } sobol_2d_t; typedef struct { uint32_t i, d0,d1,d2; } sobol_3d_t; Stratified generators start with the same data layout as progressive generators of one dimension lower followed by the element $r$ which is $\frac{1}{number\ of\ samples}$ . typedef struct { uint32_t i, d0; float r; } sobol_fixed_2d_t; typedef struct { uint32_t i, d0,d1; float r; } sobol_fixed_3d_t; typedef struct { uint32_t i, d0,d1,d2; float r; } sobol_fixed_4d_t; Initalizing Initializing a generator sets its position to zero and caller provides the initial value for the state of each dimension. This is logically equivalent to XORing this value to the given dimension result at each query. (SEE: implementation details) Standard generators void sobol_1d_init(sobol_1d_t* s, uint32_t hash) void sobol_2d_init(sobol_1d_t* s, uint32_t hash0, uint32_t hash1) void sobol_3d_init(sobol_1d_t* s, uint32_t hash0, uint32_t hash1, uint32_t hash2) Progressive generators additionally need the number of samples which will be generated. void sobol_fixed_2d_init(sobol_fixed_2d_t* s, uint32_t len, uint32_t hash) void sobol_fixed_3d_init(sobol_fixed_3d_t* s, uint32_t len, uint32_t hash0, uint32_t hash1) void sobol_fixed_4d_init(sobol_fixed_4d_t* s, uint32_t len, uint32_t hash0, uint32_t hash1, uint32_t hash2) Position in sequence manipulation Returns the current position in the sequence of the generator. This accepts any type of generator. uint32_t sobol_tell(void* s) Moves the current position by the relative value offset. Only standard generators are directly supported. These functions are not defined as inline. void sobol_1d_seek(sobol_1d_t* s, uint32_t off) void sobol_2d_seek(sobol_2d_t* s, uint32_t off) void sobol_3d_seek(sobol_3d_t* s, uint32_t off) Reading sequence members Each dimension returns results on the unit interval: $ \left[0,\thinspace1\right) $ by default. (SEE: implementation details) The 1D version returns the next member of the sequence as either a single or double and increments the position. float sobol_1d_next_f32(sobol_1d_t* s) double sobol_1d_next_f64(sobol_1d_t* s) nD versions place the $n$ dimensional elements sequentially at $d$, in either as singles or doubles, and increments the position. void sobol_2d_next_f32(sobol_2d_t* s, float* d) void sobol_3d_next_f32(sobol_3d_t* s, float* d) void sobol_2d_next_f64(sobol_2d_t* s, double* d void sobol_3d_next_f64(sobol_3d_t* s, double* d) void sobol_fixed_2d_next_f32(sobol_fixed_2d_t* s, float* d) void sobol_fixed_3d_next_f32(sobol_fixed_3d_t* s, float* d) void sobol_fixed_4d_next_f32(sobol_fixed_4d_t* s, float* d) void sobol_fixed_2d_next_f64(sobol_fixed_2d_t* s, double* d) void sobol_fixed_3d_next_f64(sobol_fixed_3d_t* s, double* d) void sobol_fixed_4d_next_f64(sobol_fixed_4d_t* s, double* d) Exposed internals The state updating functions. (SEE: implementation details) void sobol_1d_update(sobol_1d_t* s) void sobol_2d_update(sobol_2d_t* s) void sobol_3d_update(sobol_3d_t* s) Optional features The following funtions are only available with #define SOBOL_EXTRA. These additionally require a single precision square root function be defined as sqrtf. None of these functions are defined as inline. Uniform sampling of (open) unit disk ($\mathbb{D}^1$),half-disk ($x\geq 0$) and quater-disk ($x,y\geq 0$)are performed using a rejection method. This has no impact on the well-space properties. The average number of iterations is $\frac{4}{\pi} \approx 1.27324$. They both take a 2D generator $s$ and place the point in $p$. The returned value is $p \cdot p$. float sobol_uniform_d1(sobol_2d_t* s, float* p); float sobol_uniform_hd1(sobol_2d_t* s, float* p); float sobol_uniform_qd1(sobol_2d_t* s, float* p); Uniform sampling of (open) unit 3D sphere ($\mathbb{S}^2$) and half-sphere ($z\geq 0$) are generated using Marsaglia’s transform1 of the unit disk. void sobol_uniform_s2(sobol_2d_t* s, float* p); void sobol_uniform_hs2(sobol_2d_t* s, float* p); Implementation details and configuration options State update The basic construction is the classic bit-reversed gray codes2. The state in each explicit dimension is a 32-bit integer which at each update has some number of bits flipped based on the trailing one count of the position. Instead of storing and updating the explicit position, the library instead stores its bit-not and adding one becomes subtracting one ($ \sim x = -x-1 $). This allows the state update to directly use a trailing zero count (without an extra bit-not op). The input to the trailing count is only zero when an update is wrapping the position from -1 to 0, so intel’s BSF (bit scan forward) can be used and is in fact the libraries default. Configuration if a macro #define SOBOL_BIAS is defined (default is undefined) then single precision queries do not perform a bit-shift to convert the integer to floating point. This causes single results to be on the interval $ \left[0,\thinspace1\right] $, double results are uneffected. The library needs to be able to perform a trailing zero count. Currently this is only set-up for x86-a-likes. #define SOBOL_NTZ(X) number_of_trailing_zeroes_function(X) References and Footnotes Choosing a Point from the Surface of a Sphere, George Marsaglia, 1972 (PDF) ↩ Algorithm 659: Implementing Sobol’s quasirandom sequence generator, P. Bratley and B. L. Fox, 1988. ↩ Comments SHF (3) low-discrepancy (3) , sampling (7) , sequence (6)