QUASI-MONTE-CARLO METHODS This paper appears in Mathematical and Computer Modelling 23 (1996), 9–23. AND

THE DISPERSION OF POINT SEQUENCES

Günter Rote and Robert F. Tichy

Abstract. Quasi-Monte-Carlo methods are well-known for solving different problems of numerical analysis such as integration, optimization, etc. The error estimates for global optimization depend on the dispersion of the point sequence with respect to balls. In general, the dispersion of a point set with respect to various classes of range spaces, like balls, squares, triangles, axis-parallel and arbitrary rectangles, spherical caps and slices, is the area of the largest empty range, and it is a measure for the distribution of the points. The main purpose of our paper is to give a survey about this topic, including some folklore results. Furthermore, we prove several properties of the dispersion, generalizing investigations of Niederreiter and others concerning balls. For several well-known uniformly distributed point sets we estimate the dispersion with respect to triangles, and we also compare them computationally. For the dispersion with respect to spherical slices we mention an application to the polygonal approximation of curves in space.

1. Introduction

It is a classical problem in numerical analysis to find the maximum of a function. For example let f be a continuous real function defined on the s-dimensional unit cube U^{s} = [0, 1]^{s}. The following is a simple algorithm for computing its maximum value M:

Take a point sequence (x_{n}) in U^{s}, define m_{1} = f(x_{1}) and recursively set

m_{n+1} := max (m_{n}, f(x_{n+1})).

This algorithm was analyzed by Niederreiter [1], cf. also [2] and chapter 6 in the monograph [3]: The numbers m_{n} tend to the maximum M if f is “sufficiently continuous” and the points x_{n} are “well distributed” in U^{s}. More precisely,

M − ω(d_{N}) ≤ m_{N} ≤ M,

where

ω(t) := sup_{d(x,y)≤t} |f(x) − f(y)|

is the modulus of continuity of f with respect to the Euclidean distance d(x, y) and d_{N} denotes the dispersion of the sequence x_{1}, …, x_{N}:

(1.1)  d_{N} = max_{x∈U^{s}} min_{1≤n≤N} d(x, x_{n}),

1991 Mathematics Subject Classification. 05B40, 11K45, 51M16. Key words and phrases. Quasi-Monte-Carlo methods, dispersion, point sequences, uniform

distribution, polygonal approximation, range spaces.

The second author was supported by the Austrian Science Foundation (Project P10223PHY)

Typeset by AMS-TEX

which was introduced by Hlawka [4] and later investigated in more general form in Niederreiter [5]. If d_{N} is small the points x_{1}, …, x_{N} are “well distributed” in the sense that they leave no big circular hole. This measure for the distribution behavior of a point sequence is particularly suited for the maximization problem.

Another problem that can be tackled by Monte-Carlo methods is integration in higher dimensions. Here, instead of the maximum of the function values one computes their average

I_{n} := \frac{1}{n} · ∑_{i=1}^{n} f(x_{i}),

which converges to ∫_{U^{s}} f(x) dx as n goes to infinity if (x_{n}) is uniformly distributed in U^{s}. The error estimates that one gets are then proportional to the well-known discrepancy of the point set A_{N} = {x_{1}, …, x_{N}}:

(1.2)  D_{N} = sup_{R} |\frac{#(A_{N} ∩ R)}{N} − μ(R)| ,

where μ denotes the Lebesgue measure on U^{s} and the supremum is taken over all axis-parallel boxes R, see [3, 6–13].

For certain other application, the discrepancy with respect to other classes of subsets R of U^{s} are important. For example, in case of a two-dimensional problem on U^{2}, one can consider all circular disks, all squares, all triangles, or all convex sets.

The notions of dispersion and discrepancy can be defined in a very general setting. Let (X, μ) be a probability space and let R be a family of measurable subsets of X, called ranges. We call (X, R, μ) a range space. For a finite set A ⊂ X of N points, we define the dispersion

(1.3)  d_{(X,R,μ)}(A) := sup { μ(R) : R ∈ R, A ∩ R = ∅ }

and the discrepancy

(1.4)  D_{(X,R,μ)}(A) := sup \{ |\frac{#(A ∩ R)}{N} − μ(R)| : R ∈ R \} .

Usually, X and μ will be understood from the context and we will only write D_{R}(A) or even D(A), and similarly for the dispersion. Intuitively, the dispersion is the largest empty range, i. e., the largest range containing no point of A.

We will consider range spaces where X is a submanifold of some Euclidean space and μ is the corresponding normalized surface measure. The ranges R ∈ R are usually “natural” geometric regions, as in the above examples. More general ranges in arbitrary compact metric spaces were considered in Drmota and Tichy [14] in connection with the discrepancy of continuous curves.

Note that traditionally (also in (1.1) above) the dispersion has been defined as the radius of the largest empty sphere, whereas we consider the volume. This allows us to consider more general ranges than spheres.

A different approach was taken by Haussler and Welzl [6]. They were motivated by range query problems in computational geometry. A typical range query problem would be as follows: A set of N given points in the plane is to be preprocessed in such a way that queries for the number of points contained in a triangle can be answered quickly. The N points are given in advance and remain fixed, and the preprocessing time to build up a data structure is not so important. However, there are many subsequent triangle queries, and they should be efficient. Haussler and Welzl introduced ε-nets to define such a data structure. An ε-net is a subset A of X with dispersion d_{(X,R,μ)}(A) at most ε, where μ is the uniform distribution on X and R is the family of all possible query regions (triangles in our example). Haussler and Welzl showed the existence of ε-nets of small size under a very general combinatorial condition on R. Their condition involves the concept of Vapnik-Chervonenkis dimension [15]: This is the cardinality of the largest subset Y of X such that the set system { R ∩ Y : R ∈ R } contains all subsets of Y.

If the range space R has finite Vapnik-Chervonenkis dimension d then a random subset of X of size \frac{d}{ε} log \frac{1}{ε} is an ε-net with high probability. Komlós, Pach and Woeginger [16] established that this bound on the size cannot be improved in general. Their construction for the lower bound is purely combinatorial. The interesting question remains open, whether such a lower bound holds for some natural geometric range space. Very recently, Matoušek, Welzl, and Wernisch [17] extended this approach to subsets of small discrepancy.

In section 2 of this article we collect some elementary inequalities between the dispersion of point sets for various range spaces in the unit square with respect to uniform distribution. Sections 3 and 4 are devoted to the investigation of special sequences which are used for quasi-Monte-Carlo applications, such as the Halton sequence [18], the Sobol sequence [10, 19], the Faure sequence [20], and the Hammersley sequence [3]. In section 3 we prove bounds for the dispersion of a few of these point sequences in the unit square with respect to uniform distribution μ and natural geometric range spaces. Procedures for generating quasi-random point sets are given in [21–25]. Some of these papers test the quality of the point sets by computational experiments. In section 4 we test the quality of these sets by computing their dispersion in two and three dimensions.

In the final section 5 we mention an application of the dispersion with respect to slices on the two-dimensional sphere to the piece-wise linear approximation of curves in space. This problem is of importance for example in motion control in robotics. For a detailed presentation we refer to [26].

2. General estimates for the dispersion in the unit cube U^{s}

In this section we consider the uniform distribution μ in the unit cube X = U^{s} in connection with various range spaces. For any particular range space the dispersion is clearly bounded by the discrepancy:

(2.1)  d_{(X,R,μ)}(A) ≤ D_{(X,R,μ)}(A).

This can be seen by restricting the supremum in definition (1.4) to those ranges R with A ∩ R = ∅.

In the following we consider three natural classes of range spaces. We will see that the dispersion is of the same order of magnitude within each class.

1. Range spaces of fixed shape. Examples are the range space of all balls, of all axisparallel or of arbitrary cubes. In general a range space of fixed shape consists of all sets which are homothetic or which are similar to a given convex body P. This includes the range space of all balls with respect to some norm in the s-dimensional Euclidean space (restricted to U^{s}).

2. Axis-parallel range spaces. This includes the range space of all axis-parallel rectangular boxes or of all axis-parallel ellipsoids. The characteristic feature of these ranges is that they can be arbitrarily “long and thin” in each coordinate direction. In general an axis-parallel range space consists of all sets which can be obtained from a given convex body P by independent scalings of the coordinate axes.

3. Isotropic range spaces. This includes the range space of all simplices, of all ellipsoids, of all rectangular boxes, or of all general convex regions. The characteristic feature of these ranges is that they can be arbitrarily “long and thin” in any direction. This class of range spaces might be defined as those ranges spaces which consist of convex sets and are closed under affine transformations. However, the range space of rectangular boxes would not fall under this definition, and therefore we have to use a broader definition: A family R of convex sets forms an isotropic range space if every ellipsoid S contains a range R ∈ R with μ(R) ≥ cμ(S), for some constant c.

Most of the simple geometric range spaces that we could think of (as long as they include only convex ranges) fall into one of the three classes.

Proposition 2.1. Let R and S be two range spaces of the same class. Then there are constants c_{1} and c_{2} such that the following relation holds for all point sets A:

(2.2)  c_{1}d_{R}(A) ≤ d_{S}(A) ≤ c_{2}d_{R}(A).

Proof. It is sufficient to show that each range R ∈ R contains a range S ∈ S with μ(S) ≥ cμ(R), for some constant c, and vice versa.

For ranges R in a range space R of fixed shape, it is clear that this relation holds between R and a Euclidean ball S: R is similar to the convex body P, which has a largest inscribed ball S_{1} and a smallest circumscribed ball S_{2}. Thus R contains a ball S with μ(S) ≥ \frac{μ(S_{1})}{μ(P)}μ(R) and similarly any ball S contains a range R ∈ R

with μ(R) ≥ \frac{μ(P)}{μ(S_{2})}μ(S). Between two different range spaces of the first class the relation follows via the balls.

For axis-parallel range spaces, the relation (2.2) follows in the same way by replacing balls by axis-parallel ellipsoids, and for isotropic range spaces we take arbitrary ellipsoids. Each range R of an isotropic range space R is convex and therefore it contains an ellipsoid S with μ(S) ≥ μ(R)/d^{d}; the converse statement (with a different constant) follows from the definition. For particular pairs of range spaces, the determination of the optimal constants involved in these relations is a problem of its own geometric interest. For example, every convex body R contains a triangle S with area μ(S) ≥ \frac{3\sqrt{3}}{4π}μ(R); every triangle R contains a rectangle S with area μ(S) ≥ μ(R)/2; and every rectangle R contains a triangle S with area μ(S) ≥ μ(R)/2, etc., see Fejes Tóth [27].

Theorem 2.2. Let R_{1}, R_{2}, and R_{3} be range spaces of classes 1, 2, and 3 defined above, respectively. Then the following inequalities hold for all point sets A:

(2.3)  c_{1}d_{R_{1}}(A) ≤ d_{R_{3}}(A) ≤ c_{2}d_{R_{1}}(A)^{1/s}

(2.4)  c_{3}d_{R_{1}}(A) ≤ d_{R_{2}}(A) ≤ c_{4}d_{R_{1}}(A)^{1/s}

(2.5)  c_{5}d_{R_{2}}(A) ≤ d_{R_{3}}(A) ≤ c_{6}d_{R_{2}}(A)^{1/s},

where c_{i} are appropriate positive constants.

Proof. By Proposition 2.1 it is sufficient to consider a fixed representative of each class of range spaces. Thus we take as R_{1} the class of all Euclidean balls, as R_{2} the class of all axis-parallel ellipsoids, and as R_{3} the class of all convex bodies. Since R_{1} ⊂ R_{2} ⊂ R_{3}, the left-hand sides of the above relations follow trivially, and for the right-hand sides it is enough to prove (2.3). For this purpose we make use of the following lemma.

Lemma 2.3. For an s-dimensional convex body P with volume V, inradius r, and circumradius R, the following inequality holds:

V ≤ 2ω_{s−1} · rR^{s−1},

where ω_{s−1} is the volume of the (s − 1)-dimensional unit ball.

Proof. We show that the inequality is true even if we replace the inradius r by the radius r_{C} of the largest inscribed ball B which is centered at the center of gravity C of P. Let D be a point where B touches the boundary of P. Since P is convex, there is a supporting hyperplane h of P through D; since P contains B, the supporting hyperplane h must be perpendicular to the radius CD (see the left part of Figure 1, where the line through CD is taken as the x-axis).

Figure 1.

Now consider the cylinder Z which is bounded by h and by the parallel hyperplane h′ which lies symmetric to h with respect to C, and by the tangents of the circumsphere which are parallel to the x-axis. The height of this cylinder is 2r_{C}, and the volume of its basis is ω_{s−1}R^{s−1}. Hence, the volume of the cylinder is just our desired upper bound on the volume V of P.

We prove the lemma indirectly and assume the volume V is bigger than this bound. Then some part of P must lie to the right of h′. Now we construct a modified body P̃ by “moving” some part of the volume of P from the right side of h′ into the cylinder Z until this cylinder is full. Formally, we take any subset of the right part of P which has the correct measure and remove it, and we replace the left part of P by Z, see Figure 1. (We need not care about convexity at this stage of the proof.)

Since mass has been moved from the right side of h′ to the left of h′, the center of gravity of P̃ has its x-coordinate to the left of C. On the other hand, P̃ consists of two parts: Z, whose center has the same x-coordinate as C; and an additional part to the right of h′. Thus P̃’s center of gravity must lie to the right of C, a contradiction.

Now let us conclude the proof of the theorem. For the point set A we have an empty convex body of volume V = d_{R_{3}}(A) − ε, for any ε > 0. Since this body is contained in the unit cube, its circumradius R is at most \sqrt{s}/2. By the lemma, there is thus an empty ball of radius

r ≥ \frac{V}{2ω_{s−1} · (\sqrt{s}/2)^{s−1}} = c · V.

Hence the volume of this ball is at least c^{s}ω_{s}·(d_{R_{3}}(A)−ε)^{s}, for any ε > 0. From this the right-hand side of (2.3) follows, and the proof of the theorem can be completed as described above.

3. The dispersion of special sequences in the unit cube U^{s}

In this section we will establish several bounds for the dispersion of some special point sequences with lattice structure. All these sequences are multi-dimensional extensions of the well-known van der Corput sequence, which is defined by γ_{2}(n) = ∑_{j≥0} a_{j}2^{−j−1} for n = ∑_{j≥0} a_{j}2^{j} in dyadic representation. For the discrepancy D of the first N elements of this sequence with respect to intervals the following bound holds.

(3.1)  c_{1} \frac{log N}{N} ≤ D ≤ c_{2} \frac{log N}{N} ,

where c_{1}, c_{2} are suitable positive constants, see Kuipers and Niederreiter [9]. Considering general q-ary digit representations and setting

γ_{q}(n) = ∑_{j≥0} b_{j}q^{−j−1} for n = ∑_{j≥0} b_{j}q^{j},

the following s-dimensional extensions of the van der Corput sequence are known. Definition. Let q_{i} be co-prime integers ≥ 2. The Halton sequence in s dimensions [18] is the sequence

E_{N} = { (γ_{q_{1}}(n), …, γ_{q_{s}}(n)) : 0 ≤ n < N },

and the Hammersley sequence in s dimensions is the sequence

H_{N} = \{ (\frac{n}{N} , γ_{q_{1}}(n), …, γ_{q_{s−1}}(n)) : 0 ≤ n < N \}.

For the discrepancy (with respect to axis-parallel rectangles) the following estimates are known:

(3.2)  D(E_{N}) ≤ c_{3} \frac{log^{s} N}{N} ,

and

(3.3)  D(H_{N}) ≤ c_{4} \frac{(log N)^{s−1}}{N} ,

where c_{3} and c_{4} are suitable constants. For the Halton sequence E_{N} in s = 1 dimension the converse inequality of (3.2) holds also (with a different value of the constant c_{3}), and similarly, a lower bound matching the upper bound of (3.3) is known for the Hammersley sequence H_{N} in dimensions s = 1 and s = 2. Furthermore Larcher [28] proved that the isotropic discrepancy of these point sequences is of the order O(N^{−1/s}). The constants in the above estimates are superexponentially increasing in the dimension s. For getting smaller constants Sobol [10], Faure [20] and Niederreiter [29, 30] considered far-reaching extensions of these sequences, the so-called net-sequences, see also [31]. Definition. Let 0 ≤ t ≤ m be integers. A (t, m, s)-net in base q is a point set A of q^{m} points in U^{s} such that #(I ∩ A) = q^{t} for every s-dimensional interval

(3.4)  I = ∏_{i=1}^{s} [b_{i}q^{−d_{i}}, (b_{i} + 1)q^{−d_{i}})  (b_{i}, d_{i} integral)

with volume q^{− ∑ d_{i}} = q^{t−m}. Let t be a positive integer. Then a point sequence x_{1}, x_{2}, … ∈ U^{s} is called a (t, s)-sequence in base q if for all non-negative integers k and m ≥ t, the set {x_{n} : kq^{m} < n ≤ (k + 1)q^{m}} is a (t, m, s)-net in base q. Remark. For (t, m, s)-nets and for (t, s)-sequences the above bounds (3.2) and (3.3) are true with constants tending to 0 for s → ∞ if q is chosen appropriately depending on s, see [29].

Proposition 3.1. Let A be the Hammersley sequence or the Halton sequence with bases q_{1}, q_{2}, …, or any (t, m, s)-net with basis q, consisting of N points. Then the dispersion with respect to axis-parallel boxes satisfies the following estimate:

c_{5} \frac{1}{N} ≤ d(A) ≤ c_{6} \frac{1}{N}

with suitable positive constants c_{5}, c_{6}.

Proof. The lower bound d(A) ≥ 1/(N + 1) holds for any sequence. For the case of (t, m, s)-nets any interval of type (3.4) with volume q^{t}/N contains exactly q^{t} points. For the Halton sequence the elementary cells are defined in analogy to (3.4) with different bases q_{i}. In the case of the Hammersley sequence an elementary cell is also defined analogously, except that the first dimension of I is bounded by arbitrary multiples of 1/N. In any case, an elementary cell of volume ≥ \frac{1}{N} contains at least one point. Since any axis-parallel box of volume v contains an interval of type (3.4) with volume bigger than v/(2q)^{s} (or v/ ∏(2q_{i}) or (v − 2/N)/ ∏(2q_{i}), respectively), a box with volume greater than q^{t}(2q)^{s}/N (or 2^{s} ∏ q_{i}/N or (2^{s} ∏ q_{i}+ 2)/N, respectively) cannot be empty (cf. Larcher [32]).

Theorem 3.2. Let H_{N} = \{ (\frac{n}{N} , γ_{2}(n)) : n = 0, …, N −1 \} be the two-dimensional Hammersley sequence. Then the dispersion d(H_{N}) with respect to arbitrary rectangles satisfies

(3.5)  \frac{c_{7}}{\sqrt{N}} ≤ d(H_{N}) ≤ \frac{c_{8}}{\sqrt{N}}

with some positive constants c_{7}, c_{8}.

Proof. The right-hand side of (3.5) follows immediately from the above proposition and Theorem 2.2. For the left-hand side we proceed as follows:

W. l. o. g. we take N = 2^{2s} and write n < N in binary representation

n = ∑_{j=1}^{2s} e_{j}2^{j−1}.

Hence we get

x = \frac{n}{N} = ∑_{j=1}^{2s} e_{j}2^{j−1−2s} and y = γ_{2}(n) = ∑_{j=1}^{2s} e_{j}2^{−j}.

If e_{k} = e_{2s−k+1} for all k the point lies on the main diagonal y = x. Otherwise let k be the smallest index such that e_{k} ≠ e_{2s−k+1}. Then we have

|y − x| = |\frac{n}{N} − γ_{2}(n)| = |\frac{e_{k} − e_{2s−k+1}}{2^{k}} + \frac{e_{k+1} − e_{2s−k+}}{2^{k+1}} + ⋯ + \frac{e_{2s−k+1} − e_{k}}{2^{2s−k+1}} |

= |(e_{k} − e_{2s−k+1}) (\frac{1}{2^{k}} − \frac{1}{2^{2s−k+1}} ) + (e_{k+1} − e_{2s−k}) (\frac{1}{2^{k+1}} − \frac{1}{2^{2s−k}} ) +

⋯ + (e_{s} − e_{s+1}) (\frac{1}{2^{s}} − \frac{1}{2^{s+1}} )|

≥ (\frac{1}{2^{k}} − \frac{1}{2^{2s−k+1}} ) − (\frac{1}{2^{k+1}} − \frac{1}{2^{2s−k}} ) − ⋯ − (\frac{1}{2^{s}} − \frac{1}{2^{s+1}} ) ≥ \frac{1}{2^{s+1}} .

(The last sum achieves its minimum for k = s.) Thus the strip 0 < y − x < 2^{−s−1} contains no point of the point set H_{N}. It contains a rectangle of area 2^{−s−1} − 2^{−2s−2} ≈ 1/(4\sqrt{N}), and the proof of the theorem is complete.

Remark. The above theorem can easily be generalized to the two-dimensional q-ary Hammersley sequence for arbitrary q. Problem 1. It is an interesting question to ask for dispersion bounds for the Hammersley sequence in dimension s ≥ 3 as well as for the Halton sequence in dimension s ≥ 2. Definition. For N = 2^{t} the Sobol sequence S_{N} in two dimensions can be defined via the matrix M_{t} = (\binom{j}{i} mod 2)_{0≤i,j<t} of binomial coefficients modulo 2. We set

S_{N} = {( ∑_{0≤i<t} x_{i}2^{−i−1}, ∑_{0≤i<t} y_{i}2^{−i−1})},

[figure]

Figure 2. The Sobol sequence with 2^{7} points.

where the vector X = (x_{i})_{0≤i<t} runs through all 0,1-sequences of length t and

Y = (y_{i})_{0≤i<t} = M_{t}X mod 2.

Figure 2 shows the Sobol sequence with 128 points.

Theorem 3.3. Let t = 2^{k+1} −1 and let S_{N} be the two-dimensional Sobol sequence with N = 2^{t} elements. Then the dispersion d(S_{N}) with respect to rectangles satisfies

(3.6)  \frac{c_{9}}{\sqrt{N}} ≤ d(S_{N}) ≤ \frac{c_{10}}{\sqrt{N}}

with some positive constants c_{9}, c_{10}.

Proof. Since S_{N} is a (0, t, 2)-net the right-hand side of (3.6) follows immediately from Proposition 3.1 and Theorem 2.2. For the proof of the left-hand side, note that the upper-triangular matrix M_{t} for t = 2^{k+1} can be written recursively as follows:

M_{2^{k+1}} = ( \begin{matrix}M_{2^{k}} & M_{2^{k}} \\ 0 & M_{2^{k}} \end{matrix} ) .

From this (or directly from the definition) it follows easily that (M_{t})^{2} = I, for all t. Let us write the vector X in two groups: X_{0} = (x_{0}, x_{1}, …, x_{2^{k}−1}) and X_{1} = (x_{2^{k}}, …, x_{2^{k+1}−3}, x_{2^{k+1}−2}, 0) and similarly for (Y_{0}, Y_{1}) = (Y, 0). (For reasons of symmetry, we have inserted a last component x_{t} = y_{t} = 0.) Then the equations defining the Sobol sequence can be written as

Y_{0} = M_{2^{k}}X_{0} + M_{2^{k}}X_{1}
Y_{1} = + M_{2^{k}}X_{1}

Since the last component of X_{1} is 0 and the matrix M_{2^{k}} is upper triangular it follows from the first equation that the last components of X_{0} and Y_{0} are equal: x_{2^{k}−1} = y_{2^{k}−1}. If we partition the unit square into squares of size 2^{−2^{k}−2} ×2^{−2^{k}−2}, x_{2^{k}−1} = y_{2^{k}−1} means that the empty squares and the non-empty squares form a checkerboard pattern, see figure 2, where k = 2. We will show that the points in the squares on the main diagonal (X_{0} = Y_{0}) lie exactly on the main diagonal, i. e., they also fulfill X_{1} = Y_{1}: Multiplying the first equation by M_{2^{k}} and using (M_{2^{k}})^{2} = I and X_{0} = Y_{0}, we get X_{1} = (M_{2^{k}} − I)X_{0}, and the second equation gives

X_{1} − Y_{1} = (I − M_{2^{k}})X_{1} = −(M_{2^{k}} − I)^{2}X_{0} = 0.

Thus a strip of height 1/2^{2^{k}−2} = 1/2^{(t−3)/2} above and below the main diagonal is empty. It contains a rectangle of area

\frac{1}{2^{(t−3)/2}} − \frac{1}{2^{t−3}} ≈ \frac{1}{\sqrt{2N}} ,

and the lower bound is proved.

Remark. Since S_{2^{t−1}} is a subset of S_{2^{t}}, the empty strip in the above proof does not contain a point in any set S_{2^{t}}, for t ≤ 2^{k+1} − 1, and yields a lower bound for the dispersion. However, the area of this strip becomes smaller and smaller when compared to the number of points as t decreases, and thus the quality of the bound becomes worse. For powers of 2, t = 2^{k}, we only obtain the trivial bound d(S_{N}) ≥ c_{1}/N. Problem 2. The dispersion of the two-dimensional Sobol sequence S_{2^{t}} for general values of t is not known. The Sobol sequence has been generalized to other bases than 2 and to higher dimensions by Faure [20], see also [29] and the following definition in the next section. It is open what happens in these cases.

4. Computational results

In this section we give some numerical results on the dispersion of the following point sets in the plane: the Hammersley sequence for q = 2 and q = 3, the Sobol sequence, and the Halton sequence for bases q_{1} = 2 and q_{2} = 3. In three dimensions, we consider the Halton sequence for bases q_{1} = 2, q_{2} = 3 and q_{3} = 5, the Hammersley sequence for bases q_{1} = 2 and q_{2} = 3, Sobol’s dyadic net sequences [11], and the Faure sequence for q = 3. The Faure sequence is a generalization of the Sobol sequence to base q and arbitrary dimension s. Definition. Let q be a prime and q ≥ s. For N = q^{t} the Faure sequence F_{N} in s dimensions to base q can be defined via the matrix M_{t} = (\binom{j}{i} mod q)_{0≤i,j<t} of binomial coefficients modulo q. We set

F_{N} = {( ∑_{0≤j<t} x_{j}^{1}q^{−j−1}, ∑_{0≤j<t} x_{j}^{2}q^{−j−1}, …, ∑_{0≤j<t} x_{j}^{s}q^{−j−1}, )},

where the vector X^{1} = (x_{j}^{1})_{0≤j<t} runs through all sequences of length t with digits from the set {0, 1, …, q − 1}, and

X^{i} = (x_{j}^{i})_{0≤j<t} = M_{t}X^{i−1} mod q, for i = 2, …, s.

In the plane, we constructed the respective sequences for all powers of the appropriate bases (2 or 3) up to N = 3^{12} = 531 441. For each point set we computed the largest empty triangle that is formed by points of the sequence, because this is easier to compute than the largest empty triangle T with arbitrary corners. Since a largest empty triangle T with arbitrary corners always contains a triangle whose vertices belong to the given set and whose area is at least a quarter of the area of T, the numbers in tables 1 and 2 give the dispersion with respect to triangles up to a constant factor. By Proposition 2.1 this also holds for arbitrary isotropic range spaces. The numbers in the tables are the areas multiplied by N to make them comparable for different values of N.

We describe now the algorithm that we use to compute the largest empty triangle. We go through all empty triangles and keep the largest one. In order to enumerate all empty triangles, we fix the leftmost vertex v_{1} of the triangle. We declare that the segment v_{2}v_{3} between two points to the right of v_{1} is visible (from v_{1}) if the triangle v_{1}v_{2}v_{3} is empty. In fact, if we look at the simple polygon formed by v_{1} and the points v to the right of v_{1} in angular order, the visible segments together with the segments from v_{1} to all other vertices form the visibility graph of this polygon, a well-known concept in computational geometry. The visibility graph of a general simple polygon can be found in time proportional to its size [33, 34]. For computing the visibility graph one needs a triangulation of the polygon. In our case the polygon is star-shaped with respect to v_{1}, and a triangulation is readily available. Thus the algorithm becomes quite simple. We give here only an intuitive description of the ideas behind the algorithm. Imagine spanning a rubber string from v_{2} via v_{1} to v_{3} and letting it loose at v_{1}. The string will form a convex piecewise linear curve between v_{2} and v_{3}, and it will be a straight segment if and only if v_{2}v_{3} is visible. The union of all these strings from a fixed vertex v_{2} is a tree rooted at v_{2}.

We sweep through the points v_{2} to the right of v_{1} in angular order around v_{1}. As we go, we maintain the tree of rubber strings from v_{2} to all points v_{3} preceding v_{2} in the angular order. If we go from a vertex v′_{2} to its successor v_{2} in the angular order, the new tree edges in the tree rooted at v_{2} that were not contained in the old tree rooted at v′_{2} are precisely the tree edges out of v_{2}. It is possible to update the tree in time that is proportional to the number of these new edges, by exploring a subtree of the old tree starting at v′_{2}. On the other hand, the tree edges out of v_{2} correspond to the empty triangles v_{1}v_{2}v_{3}, and thus, the total work involved in updating the trees is proportional to the number of empty triangles.

Before we start the sweep, we have to perform an angular sort around each vertex v_{1}, at a total cost of N · O(N log N). (It is possible to compute the sorted angular lists around all points in O(N^{2}) total time, but only at the cost of O(N^{2}) storage.) Thus the overall time for our algorithm is proportional to N^{2} log N plus the number of empty triangles. Since we expect our sequences to be “well” distributed, we may reasonably assume that the number of empty triangles is small compared to the maximum possible number of O(N^{3}). In our experiments, the number of empty triangles appeared to be growing at a rate between N^{2.02} and N^{2.1}, which is reasonably close to N^{2}. (The exponents were obtained by a logarithmic regression.) This is in accordance with a result of Bárány and Füredi [35] who showed that, for N random points, the expected number of empty triangles is O(N^{2}).

The programs were written in Pascal. The codes are available from the authors.

N × (area of the largest empty triangle)

N  Hammersley sequence  Sobol Halton seq.
q = 3  q = 2  sequence q_{1}=2, q_{2}=3
3 = 3^{1}  0.000000  0.375000
4 = 2^{2}  0.375000 1.000000 0.500000
8 = 2^{3}  1.312500 1.750000 1.166667
9 = 3^{2} 1.777778 1.812500
16 = 2^{4}  2.343750 3.500000 2.888889
27 = 3^{3} 3.851852 2.296875
32 = 2^{5}  4.359375 4.000000 2.901235
64 = 2^{6}  4.359375 6.562500 4.475309
81 = 3^{4} 5.135802 4.875000
128 = 2^{7}  5.953125 7.000000 5.419753
243 = 3^{5} 11.950617  5.966797
256 = 2^{8}  6.445312 8.562500 6.286008
512 = 2^{9}  11.976563 15.750000 9.607682
729 = 3^{6} 11.950617  9.216797
1024 = 2^{10}  11.976563 22.359375 10.764060
2048 = 2^{11}  23.988281 30.312500 12.506401
2187 = 3^{7} 35.983539  11.141357
4096 = 2^{12}  23.988281 52.939453 14.043439
6561 = 3^{8} 35.983539  14.584961
8192 = 2^{13}  47.994141 51.767578 15.343520
16384 = 2^{14}  47.994141 90.667969 19.415231
19683 = 3^{9} 107.994513  20.362000
32768 = 2^{15}  95.997070 89.997253 18.413529
59049 = 3^{10} 107.994513  25.303238
65536 = 2^{16}  95.997070 39.998779 28.082999
131072 = 2^{17}  191.998535 80.078125 23.567032
171147 = 3^{11} 323.998171  29.122341
262144 = 2^{18}  191.998535 159.370117 30.958357
524288 = 2^{19}  383.999268 317.802734 28.360017
531441 = 3^{12}  28.746941

Table 1. Largest-area empty triangles in quasi-random point sets

We now discuss the results of table 1. We see that for the Hammersley sequence the area of the largest triangle decreases by a factor of q when N increases from an odd power of q to the subsequent even power (except for small N). In table 1, the corresponding entries are equal. On the other hand, when N increases from an even power of q to the subsequent odd power, the area of the largest empty triangle remains almost unchanged. (The corresponding entries in the table are approximately multiplied by q.) We conjecture that this is true in general. The data seem to support the hypothesis that the area of the largest empty triangle for the Hammersley sequence with N = 3^{n} and N = 2^{n} points is asymptotically equal to (4/9) · 3^{−⌊n/2⌋} and (3/8) · 2^{−⌊n/2⌋}, respectively. This would be in accordance with Theorem 3.2.

The figures for the Sobol sequence exhibit a similar O(\sqrt{N}) growth rate, but the behavior is more erratic. The numbers make large jumps from N = 2^{8} to 2^{9}, from 2^{13} to 2^{14}, and from 2^{16} to 2^{17}. On the other hand there is a small decrease from 2^{12} to 2^{13} and a huge decrease from 2^{15} to 2^{16}. Note that the case N = 2^{15} is covered by the lower bound of Theorem 3.3, and thus it is not surprising that this sequence should have quite large empty triangles.

The Halton sequence is clearly the best-distributed one. The growth of the numbers is very slow, but it seems that it is slightly faster than logarithmic. Also, the largest empty triangle was always unique for the Halton sequence, whereas there were many largest empty triangles for the Sobol sequence and even more for the Hammersley sequence.

Table 2 shows the area of the largest empty triangle for random (0, m, 2)-nets in base 2. Such a net is obtained by starting with an arbitrary (0, m, 2)-net whose point coordinates are multiples of 2^{−m}. Then we look at each pair of adjacent “elementary” vertical strips [b2^{−d}, (b + 1/2)2^{−d}] × [0, 1) and [(b + 1/2)2^{−d}, (b + 1)2^{−d}] × [0, 1), for 0 ≤ d < t and 0 ≤ b < 2^{d}, and we exchange the points in these two strips with probability 1/2 (independently of each other). This amounts to a change of x-coordinates of the points. It is clear that the properties of a net are maintained, and one can see that every possible (0, m, 2)-net is obtained with the same probability.

N × (largest empty triangle): average, minimum, maximum

N  random nets  random permutations
4 0.700000 (0.375000–0.875000)  0.600000 (0.875000–1.000000)
8 1.450000 (1.062500–2.000000)  1.331250 (0.875000–1.875000)
16 2.075000 (1.812500–2.343750)  2.812500 (1.875000–3.687500)
32 3.350000 (2.859375–3.796875)  3.985937 (2.500000–6.000000)
64 4.354687 (3.640625–5.210938)  5.685156 (4.921875–8.265625)
128 5.639062 (4.960938–6.492188)  7.344141 (5.675781–9.113281)
256 6.941992 (5.949219–8.439453)  9.574609 (7.968750–12.210938)
512 8.413672 (7.693359–9.249023) 11.033594 (10.239258–12.608398)
1024 10.262158 (9.305664–11.684082) 13.008350 (10.472656–15.835449)
2048 12.475562 (11.070313–15.442871) 15.529590 (13.556152–16.931152)
4096 13.878772 (12.431396–15.653809) 16.201721 (15.084106–17.378418)
8192 16.067871 (14.014038–19.222839) 18.151520 (16.203369–20.549500)
16384 18.094440 (16.801270–19.333923) 19.526987 (17.706146–22.538086)

Table 2. Largest-area triangles in random nets and random permutations. For each value of N we carried out 10 independent runs.

For comparison, the second column shows the result for point sets which form a random permutation between the x-coordinates {0, 1/N, 2/N, …, (N − 1)/N} and the y-coordinates from the same set. It is well-known that the largest empty triangle in a set of n points drawn independently from a uniform distribution in U^{s} has an expected area of Θ(log N/N), as follows from general results, see for example [15] or [6]. Since our program was designed to work only for permutations, we chose random permutations as an approximation of independent uniform points. Each entry represents the average of 10 runs, and the maximum and minimum values are also indicated. One can see that random nets give a slightly better distribution than the Halton sequence. Random permutations are not quite as good, but are still definitely better than the Hammersley sequences. The largest empty triangles were usually unique, except for a few cases with small numbers of points.

Table 3 presents the analogous results for three dimensions. We computed the largest empty tetrahedron with vertices from the given point set, for the same reasons as in the planar case. For generating the point sequences we used codes from Sobol and Shukhman. In three dimensions we apply a more brute-force method than in the plane. We generate all \binom{n}{4} tetrahedra and check for each one whether it is empty. For this check, we use the multidimensional binary search tree data structure (3-D tree) of Bentley [36], see also [37]: The element whose x-coordinate is the median of all x-coordinates is stored in the root node, and the remaining points are split into those with smaller x-coordinates and those with larger x-coordinates. The two halves become the two sub-trees of the root. At the next level, each half is again split, this time according to y-coordinates, and the median element is stored in the root of the subtree. At the third level, the points are split according to z-coordinates, and in the lower levels, the splitting direction continues to cycle through the three coordinate directions. Finally, if only a single point remains from the splitting it forms a leaf of the data structure. For each node we also store the coordinates of the axis-parallel box enclosing all points in its subtree. A search for a point contained in a given tetrahedron proceeds from the root to the leaves. In general, when we search a tree node, we have to recursively search both subtrees. When we find that the enclosing box of a node is disjoint from the tetrahedron in question, we can abandon the search of the whole subtree. We first search that subtree whose enclosing box contains the center of gravity of the tetrahedron. This is intended to direct the search quickly onto those regions where a point can be found. On average, less than 9 tree nodes had to be visited in our experiments before a point in the tetrahedron was found or the tetrahedron was established to be empty. (We checked only those tetrahedra which were at least as large as the largest tetrahedron found so far.)

N × (area of the largest empty tetrahedron)

N  Halton Hammersley  Sobol  N  Faure
q⃗ = (2, 3, 5) q⃗ = (2, 3)  q = 2  q = 3
4 0.08888889 0.11111111 0.00000000
8 0.44444444 0.33333333 0.66666667
16 1.77185185 0.94444444 1.26041667
32 2.80059259 2.71141975 2.47916667 9 0.50000000
64 4.83634568 4.17129630 5.29166667 27 3.16666667
128 9.57353086 5.56854424 7.07812500 81 7.77777778
256 15.97181893 8.51485340 8.36421712 243 14.91906722
512 32.89752757 12.73225309 11.51106771 729 50.84339278

Table 3. Largest-volume empty tetrahedra in three-dimensional quasi-random point sets. q⃗ stands for the sequence of bases (q_{1}, q_{2}, …).

The Hammersley and Sobol sequences seem to be best, and they beat the Halton and Faure sequences clearly. In contrast to the planar case, the Sobol sequence also seems to be competitive, and it tends to become better as the number of points grows. The largest empty tetrahedron was usually unique, except when the number of points was small.

5. Spherical dispersions with an application to polygonal approximation of curves

In this section we consider the s-dimensional unit sphere X = S^{s} with the uniform distribution μ. The range space that is usually considered on the sphere is the range space C of spherical caps, i. e., intersections of the sphere with half-spaces. For quasi-Monte-Carlo methods on the sphere we refer to Tichy [38].

In the following we will consider the dispersion with respect to the range space S of spherical slices and we will explain an application of spherical dispersions for this range space to the piecewise linear approximation of curves in space.

Obviously, as in the case of the unit cube, the dispersion does not exceed the corresponding discrepancy. The following elementary relation holds between cap and slice dispersion.

Proposition 5.1. For any point set A on the sphere S^{s}, we have

d_{S}(A) ≤ c · d_{C}(A)^{1/s},

for some positive constant c.

Theorem 5.2. For every N there is a point set A on the sphere S^{s} with slice dispersion

d_{S}(A) = O(1/N).

A detailed proof is given in [26]. One simply has to distribute the N points as evenly as possible over the \binom{s+1}{2} “coordinate circles” which are obtained as the intersections of the coordinate planes with the unit sphere, and place the corresponding points equidistantly on each circle.

In the remainder of this section we will briefly show how the slice dispersion on the sphere in three-dimensional space arises in a problem of piecewise linear approximation of curves in space.

For instance, in robotics it is an important problem to approximate a “general” curve by simple curves like straight lines, circles etc., because the arm of the robot can only run along such simple curves. The most important case is the approximation by a polygonal line through suitably chosen interpolation points on the curve.

Let us first consider a plane curve. Suppose we have already selected a starting sequence of points P_{1}, …, P_{k−1}. We can then try to advance the new point P_{k} on the curve as far as possible while still maintaining the property that the curve between P_{k−1} and P_{k} lies in an ε-strip around the the segment P_{k−1}P_{k}, for a given error bound ε. The algorithmic details of this procedure are quite straightforward and can be found in [26]. If one wants to carry over this method to three dimensions one runs into some geometric difficulties. However, we can reduce this problem to the planar case by considering suitably chosen projections. We project the curve orthogonally onto N different planes and ensure that in each projection the projected curve between P_{k−1} and P_{k} lies in an ε′-strip around the the projected segment P_{k−1}P_{k}, for a suitably chosen error bound ε′ < ε. This guarantees that the original three-dimensional curve between P_{k−1} and P_{k} lies in an aε′-strip around the the segment P_{k−1}P_{k}, for some constant a which we would like to be as small as possible.

Let us see how a depends on the choice of the N projection directions π_{1}, …, π_{N}. Let g be the line through P_{k−1} and P_{k}. We know that the curve segment between P_{k−1} and P_{k} lies in the set

Z = Z(π_{1}, …, π_{N}; g) = ⋂_{i=1}^{N} { x : dist(π_{i}(x), π_{i}(g)) ≤ ε′ },

where π_{i}(x) denotes the projection of x along the direction π_{i} onto an orthogonal plane, and dist(P, ℓ) denotes the distance of the point P from the line ℓ. The set Z is thus an intersection of N parallel slabs. To calculate the maximum distance of a point in Z from g, we intersect Z with an orthogonal plane g^{⊥} of g and obtain a convex symmetric polygon Z whose edges are parallel to the projections of the directions π_{i} onto g^{⊥}. The distance of the edges to the center of Z is ε′. The maximal distance of a vertex of Z to the center is

\frac{ε′}{cos α/2},

where α is the maximal angle between two adjacent edges. If we consider the projection directions π_{1}, …, π_{N} and the line g as point pairs on the sphere S^{2}, α is the opening angle of the largest empty slice with corners at the two points corresponding to g. Since we want α to be small for all directions g, we have to choose these N points exactly in such a way that the slice-dispersion d_{S} is minimal.

Thus, in order to ensure that the spatial curve lies in an ε-neighborhood of the piece-wise linear curve through the points P_{1}, P_{2}, …, we have to set ε′ = ε · cos(d_{S}({±π_{1}, …, ±π_{N}})/2).

6. Concluding Problem

As mentioned in the introduction, there is a lower bound of Ω(\frac{1}{ε} log \frac{1}{ε}) for the size of an ε-net in range spaces of finite Vapnik-Chervonenkis dimension. It is conjectured that this lower bound is true even for simple geometric range spaces, like triangles with the uniform distribution μ. We remark that this would imply a solution to a well-known problem of Danzer (see Beck and Chen [39], p. 285).

Acknowledgement. We are indebted to G. Larcher and Raimund Seidel for valuable discussions and comments and to J. Thuswaldner for help with implementing the programs.

References

[1] H. Niederreiter, Quasi-Monte Carlo methods for global optimization, Mathematical Statistics and Applications, Vol. B, Proc. 4^{th} Pannonian Symp., Bad Tatzmannsdorf, 1983 (W. Grossmann, G. C. Pflug, I. Vincze, and W. Wertz, ed.), Reidel, Dordrecht, 1985, pp. 251–267.

[2] H. Niederreiter and P. Peart, Localization of search in quasi-Monte Carlo methods for global optimization, SIAM J. Sci. Stat. Comput. 7 (1986), 660–664.

[3] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM, Philadelphia, 1992.

[4] E. Hlawka, Abschätzung von trigonometrischen Summen mittels diophantischer Approximation, Österreich. Akad. Wiss. Math.-Natur. Kl. Sitzungsber. II 185 (1976), 43–50.

[5] H. Niederreiter, A quasi Monte-Carlo method for the approximate evaluation of the extreme values of a function, Studies in Pure Mathematics (P. Erdös, ed.), Birkhäuser, Basel, 1983, pp. 523–529.

[6] D. Haussler and E. Welzl, ε-nets and simplex range queries, Discrete Comput. Geom. 2 (1987), 127–151.

[7] N. M. Korobov, The approximate computation of multiple integrals, Dokl. Akad. Nauk SSSR 124 (1959), 1207–1210. (Russian)

[8] N. M. Korobov, Number-Theoretic Methods in Approximate Analysis, Fizmatgiz, Moscow, 1963. (Russian)

[9] L. Kuipers and H. Niederreiter, Uniform Distribution of Sequences, Wiley, New York, 1974.

[10] I. M. Sobol′, On the distribution of points in a cube and the approximate evaluation of integrals, Zh. Vychisl. Mat. i Mat. Fiz. 7 (1967), 784–802 (Russian); English translation in U.S.S.R. Comput. Math. and Math. Phys. 7 (1967), 86–112.

[11] I. M. Sobol′, Multidimensional quadrature formulas and Haar functions, Nauka, Moscow, 1969. (Russian)

[12] E. Hlawka, Funktionen von beschränkter Variation in der Theorie der Gleichverteilung, Ann. Mat. Pura Appl. 54 (1961), 325–333.

[13] E. Hlawka, Zur angenäherten Berechnung mehrfacher Integrale, Monatsh. Math. 66 (1962), 140–151.

[14] M. Drmota and R. F. Tichy, C-uniform distribution on compact metric spaces, J. Math. Anal. Appl. 129 (1988), 284–292.

[15] V. N. Vapnik and A. Ya. Chervonenkis, On the uniform convergence of relative frequencies of events to their probabilities, Teor. Veroyatnost. i Primenen. 16 (1971), 264–279 (Russian); English translation in Theory Probab. Appl. 16 (1971), 264–280.

[16] J. Komlós, J. Pach, and G. Woeginger, Almost tight bounds for epsilon-nets, Discrete Comput. Geom. 7 (1992), 163–174.

[17] J. Matoušek, E. Welzl, and L. Wernisch, Discrepancy and approximations for bounded VCdimension, Combinatorica 13 (1993), 455–466.

[18] J. H. Halton, On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals, Numer. Math. 2 (1960), 84–90.

[19] I. M. Sobol′, Uniformly distributed sequences with an additional property of uniformity, Zh. Vychisl. Mat. i Mat. Fiz. 16 (1976), 1332–1337, 1375 (Russian); English translation in U.S.S.R. Comput. Math. and Math. Phys. 16 (1976), 236–242.

[20] H. Faure, Discrépance de suites associées à un système de numération (en dimension s), Acta Arith. 41 (1982), 337–351.

[21] J. H. Halton and G. B. Smith, Algorithm 247; radical-inverse quasi-random point sequence, Commun. ACM 7 (1964), 701–702.

[22] B. L. Fox, Algorithm 647: Implementation and relative efficiency of quasirandom sequence generators, ACM Trans. Math. Software 12 (1986), 362–376.

[23] P. Bratley and B. L. Fox, Algorithm 659: Implementing Sobol’s quasi-random sequence generator, ACM Trans. Math. Software 14 (1988), 88–100.

[24] P. Bratley, B. L. Fox, and H. Niederreiter, Implementation and tests of low-discrepancy sequences, ACM Trans. Modeling and Computer Simulation 2 (1992), 195–213.

[25] I. M. Sobol and B. V. Shukhman, On computational experiments in uniform distribution, Österreich. Akad. Wiss. Math.-Natur. Kl. Sitzungsber. II 210 (1992), 161–167.

[26] G. Rote and R. F. Tichy, Spherical dispersion with an application to polygonal approximation of curves, Anz. Österreich. Akad. Wiss. Math.-Natur. Kl. Abt. II 132 (1995), 3–10.

[27] L. Fejes Tóth, Lagerungen in der Ebene, auf der Kugel und im Raum, 2nd ed., Springer-Verlag, Berlin, Heidelberg, 1972.

[28] G. Larcher, Über die isotrope Diskrepanz von Folgen, Arch. Math. (Basel) 46 (1986), 240–249.

[29] H. Niederreiter, Point sets and sequences with small discrepancy, Monatsh. Math. 104 (1987), 273–337.

[30] H. Niederreiter, Low-discrepancy and low-dispersion sequences, J. Number Th. 30 (1988), 51–70.

[31] I. M. Sobol′, On functions satisfying a Lipschitz condition in multidimensional problems of numerical mathematics, Dokl. Akad. Nauk SSSR 293 (1987), 1314–1319 (Russian); English translation in Soviet Math. Doklady 35 (1987), 466–470.

[32] G. Larcher, The dispersion of a special sequence, Arch. Math. (Basel) 47 (1986), 347–352.

[33] J. Hershberger, Finding the visibility graph of a simple polygon in time proportional to its size, Proc. 3rd Annu. Sympos. Comput. Geom., Assoc. Comput. Mach., New York, 1987, pp. 11–20.

[34] J. Hershberger, An Optimal Visibility Graph Algorithm for Triangulated Simple Polygons, Algorithmica 4 (1989), 141–155.

[35] I. Bárány and Z. Füredi, Empty simplices in Euclidean spaces, Canad. Math. Bull. 30 (1987), 436–445.

[36] J. L. Bentley, Multidimensional binary search trees used for associative searching, Commun. ACM 18 (1975), 509–517.

[37] F. P. Preparata and M. I. Shamos, Computational Geometry: an Introduction, Springer-Verlag, New York, 1985.

[38] R. F. Tichy, Random points on the sphere with applications to numerical analysis, Z. angew. Math. Mech. 70 (1990), 642–646.

[39] Jószef Beck and William Chen, Irregularities of Distribution, Cambridge University Press, Cambridge, 1987.

Institut für Mathematik, Technische Universität Graz, Steyrergasse 30, A-8010 Graz, Austria

e-mail: rote@opt.math.tu-graz.ac.at, tichy@weyl.math.tu-graz.ac.at