How can one generate an open-ended sequence of low-discrepancy points in 3D?

I’d like a low-discrepancy sequence of points over a 3D-hypercube [1,1]3, but don’t want to have to commit to a fixed number n of points beforehand, that is just see how the numerical integration estimates develop with increasing numbers of low-discrepancy points.

I’d like to avoid have to start all over again, if the results with a fixed n are unsatisfactory. Of course, one could just employ random numbers, but then the convergence behavior would be poorer.

“A sequence of n-tuples that fills n-space more uniformly than uncorrelated random points, sometimes also called a low-discrepancy sequence. Although the ordinary uniform random numbers and quasirandom sequences both produce uniformly distributed sequences, there is a big difference between the two.” (

This question has also just been put on the

Since in his answer below, Martin Roberts advances a very interesting, appealing approach to the open-ended low-discrepancy problem, I’d like to indicate an (ongoing) implementation of his approach I’ve just reported in . In sec. XI (p. 19) and Figs. 5 and 6 there, I analyze two problems—one with sampling dimension d=36 and one with d=64—both using the parameter α0 set to 0 and also to 12. To convert the quasi-uniformly distributed points yielded by the Roberts’ algorithm to quasi-uniformly distributed normal variates, I use the code developed by Henrik Schumacher in his answer to


As the OP cross-posted this question from Math stackexchange, I have also cross-posted the answer that I wrote there.


The simplest traditional solution to the d-dimensional which provides quite good results in 3-dimensions is to use the Halton sequence based on the first three primes numbers (2,3,5). The Halton sequence is a generalization of the 1-dimensional Van der Corput sequence and merely requires that the three parameters are pairwise-coprime. Further details can be found at the Wikipedia article: “Halton Sequence”.

An alternative sequence you could use is the generalization of the Weyl / Kronecker sequence. This sequence also typically uses the first three prime numbers, however, in this case they are chosen merely because the square root of these numbers is irrational.

However, I have recently written a detailed blog post, “The Unreasonable effectiveness of Quasirandom Sequences, on how to easily create an open-ended low discrepancy sequences in arbitrary dimensions, that is:

  • algebraically simpler
  • faster to compute
  • produces more consistent outputs
  • suffers less technical issues

than existing existing low discrepancy sequences, such as the Halton and Kronecker sequences.

The solution is an additive recurrence method (modulo 1) which generalizes the 1-Dimensional problem whose solution depends on the Golden Ratio. The solution to the d-dimensional problem, depends on a special constant ϕd, where ϕd is the unique positive root of xd+1=x+1.

For d=1ϕ1=1.618033989…, which is the canonical golden ratio.

For d=2, ϕ2=1.3247179572…, which  is often called the plastic constant, and has some beautiful properties. This value was conjectured to most likely be the optimal value for a related two-dimensional problem [Hensley, 2002]. Jacob Rus has posted a beautiful visualization of this 2-dimensional low discrepancy sequence, which can be found here.

And finally specifically relating to your question, for d=3, ϕ3=1.2207440846…

With this special constant in hand, the calculation of the n-th term is now extremely simple and fast to calculate:


Of course, the reason this is called a recurrence sequence is because the above definition is equivalent to

In nearly all instances, the choice of αα0 does not change the key characteristics, and so for reasons of obvious simplicity, αα0=00 is the usual choice. However, there are some arguments, relating to symmetry, that suggest that αα0=1/21/2 is a better choice.

Specifically for d=3, ϕ3=1.2207440846… and so for αα0=(1/2,1/2,1/2),
αα=(0.819173,0.671044,0.549700) and so the first 5 terms of the canonical 3-dimensional sequence are:

  1. (0.319173, 0.171044, 0.0497005)
  2. (0.138345, 0.842087, 0.599401)
  3. (0.957518, 0.513131, 0.149101)
  4. (0.77669, 0.184174, 0.698802)
  5. (0.595863, 0.855218, 0.248502)

Of course, this sequence ranges between [0,1], and so to convert to a range of [-1,1], simply apply the linear transformation x:=2x1. The result is

  1. (-0.361655, -0.657913, -0.900599)
  2. (-0.72331, 0.684174, 0.198802)
  3. (0.915035, 0.0262616, -0.701797)
  4. (0.55338, -0.631651, 0.397604)
  5. (0.191725, 0.710436, -0.502995),…

The Mathematica Code for creating this sequence is as follows:

f[n_] := x /. FindRoot[x^(1 + n) == x + 1, {x, 1}];

d = 3;
n = 5

gamma = 1/f[d];
alpha = Table[gamma^k , {k, Range[d]}]
ptsPhi =  Map[FractionalPart, Table[0.5 + i alpha, {i, Range[n]}], {2}]

Similar Python code is

# Use Newton-Rhapson-Method
def gamma(d):
    for i in range(20):
        x = x-(pow(x,d+1)-x-1)/((d+1)*pow(x,d)-1)
    return x


g = gamma(d)
alpha = np.zeros(d)                 
for j in range(d):
    alpha[j] = pow(1/g,j+1) %1
z = np.zeros((n, d))    
for i in range(n):
    z = (0.5 + alpha*(i+1)) %1


Hope that helps!

Source : Link , Question Author : Paul B. Slater , Answer Author : Martin Roberts

Leave a Comment