apollonian.fm

Appeared in Fractals, Vol. 2, No. 4 (1994) 521-526

THE FRACTAL DIMENSION OF THE

APOLLONIAN SPHERE PACKING

M. BORKOVEC and W. DE PARIS

Institute of Terrestrical Ecology, Federal Institute of Technology ETH,

Grabenstrasse 3, 8952 Schlieren, Switzerland

R. PEIKERT

Interdisciplinary Project Center for Supercomputing,

Federal Institute of Technology ETH, Clausiusstrasse 59, 8092 Zürich

Received January 15, 1993; Accepted April 12, 1993

Abstract

The fractal dimension of the Apollonian sphere packing has been computed numerically up
to six trusty decimal digits. Based on the 31 944 875 541 924 spheres of radius greater than
2

-19

contained in the Apollonian packing of the unit sphere, we obtained an estimate of

2.4739465, where the last digit is questionable. Two fundamentally different algorithms
have been employed. Outlines of both algorithms are given.

1. INTRODUCTION

For the study of the classical Apollonian packing, it
is convenient to treat the more general osculatory
packings. In two dimensions, an osculatory packing
of the unit disk S

is given by three pairwise exter-

nally tangent disks S

, S

which are all inter-

nally tangent to S

. More (open) disks S

, ... are

added in a way that each S

is contained in S

, dis-

joint from S

, ... S

−

and of maximal size. If S

, S

and S

are chosen to be of equal size, the Apollo-

nian packing is obtained (see Fig. 1). The disks S

, S

are called the initial circles of the pack-

ing.

In the n-dimensional setting, the number of initial

circles is n + 2. It is known that osculatory packings
are complete, i.e. the sum of the area of S

, S

, ...

converges to the area of S

. Mandelbrot

defines the

fractal dimension D of the packing as the Haus-
dorff-dimension of the residual set

. A

remarkable theorem due to Boyd

says that D is the

same for all osculatory packings. It is generally
believed that D cannot be found analytically.

∞

∪

Rigorous bounds are given by Melzak

and Boyd

for two dimensions and by Boyd and Larman

for

three dimensions.

Numerical estimates can be obtained by restrict-

ing the packing to spheres of curvature (inverse
radius) not exceeding a given value

. Let

)

denote this set, N(

) its cardinality, s(

) the sum of

perimeters and p(

) the area of the set

. Then, two theorems due to Boyd

5,6

imply the following asymptotic behavior of the
three functions:

(1)

Manna and Herrmann

used these expressions to

obtain numerically the value for the fractal dimen-
sion of 1.30568... in two dimensions. They remark
that the values based on s(

) converge slower. This

is, however, just the effect of not having included
the circle S

Boyd was the first to attack numerically the 3-

dimensional problem (see Fig. 2). In Ref 8, he con-
jectures a value of approx. 2.42. It is our signifi-
cantly different result which motivated us to write
this note.

Although the definition of osculatory packings is

constructive, an efficient algorithm cannot be
derived from it. Instead, two fast algorithms are
sketched in the next two sections.

2. THE “INVERSION ALGORITHM”

As Mandelbrot points out in Ref. 1, 2-dimen-

sional osculatory packings can be constructed by
repeatedly applying circular inversion to the initial
circles. A set of four inversion circles is sufficient,
namely those keeping fixed three out of four initial
circles (see Fig. 1).

We notice that every packing circle which is not

an initial circle lies completely inside one and out-
side all other inversion circles. Therefore, it has
exactly one larger image. This observation allows
to construct a simple algorithm for generating each
circle exactly once: Starting with the initial circles,
apply to each circle exactly those of the four inver-
sions that reduce its size. If only circles of a certain
minimum size are to be generated, this process can
be made finite by cutting each computation branch
whenever a smaller circle has been produced.
Finally, by assigning distinct colors to the initial
circles, the packing can be divided into four clans,
i.e. subsets of constant color (as shown in Fig. 2).

In three dimensions, inversion spheres overlap so

that some of the packing spheres have two larger
images. The naïve adaptation of the above algo-
rithm would therefore generate all spheres, how-
ever some of them more than once (but with
consistent colors). To restore the uniqueness we
allow not all size-reducing inversions, but only
those which produce spheres having their center in
the “target region” of the respective inversion
sphere. The target regions are the insides of the
inversion spheres minus parts of the regions of
overlap. They are chosen such that their disjoint
union equals the union of the inversion spheres. It is
this simple algorithm that we have used to obtain
numerically the fractal dimension of the three-
dimensional Apollonian packing.

As a matter of efficiency, it is advisable to repre-

sent spheres not with Cartesian coordinates x

,...,

and a radius r but with inversive coordinates

(2)

In these coordinates, the inversion becomes a lin-

ear transformation of the coordinate vector. Note
that also the curvature can be expressed linearly
(namely a

n+2

−

n+1

), and so can the criterium for

the target regions.

To further simplify calculations, coordinates can

be altered again by the linear transformation

(3)

leading to integer coordinates for the initial spheres

(4)

( )

–

∪

( ) κ

∼

( ) κ

–

∼

( ) κ

–

∼

-----

…

, ,

-----

…

–

--------------------------------------------------

…

–

--------------------------------------------------

-------

–

-------

–

------- 0

-------

–

-------

–

------- 0

-------

–

-------

–

-------

------- 0

–

-------

------- 0

0 2 2 2 1

, , , ,

〈

〉

2 0 2 2 1

, , , ,

〈

〉

2 2 0 2 1

, , , ,

〈

〉

2 2 2 0 1

, , , ,

〈

〉

0 0 0 0

–

, , , ,

〈

〉

and simpler linear mappings for the spherical inver-
sions:

(5)

3. THE “SODDY ALGORITHM”

In n-space, the maximum number of mutually tan-
gent spheres is n + 2, if one demands that no three
spheres have a common point. We will call such a
set of spheres a Soddy set. Soddy's formula

is an

identity for the curvatures in a Soddy set:

(6)

If it is solved for one of the

, we get exactly two

distinct solutions. Hence, it is always possible to
replace any sphere in a Soddy set by its “second
solution”, the result is again a Soddy set. We might
call this a Soddy operation. Since the sum of the
solutions of Eq. (6) is:

(7)

Soddy operations can be computed as linear func-
tions on n + 2-vectors (assuming that coordinates
are of no interest).

Applying these definitions to osculatory pack-

ings, we notice that the set of initial spheres is a
Soddy set. A method to produce all spheres of the
packing is to apply to the set of initial spheres all
n + 2 Soddy operations and to all its successors
those n + 1 which do not invert the last operation.

In two dimensions, this procedure has two prop-

erties:

(i) Each Soddy operation generates a circle which

is smaller than the four circles in the Soddy
set.

(ii) The same circle will not be generated twice. If

the goal is to enumerate the circles of a certain
minimum size only, we can cut each branch of
the computation as soon as a smaller circle has
been generated. This algorithm has been used

by Melzak

In three dimensions, neither (i) nor (ii) hold.

Based on geometrical evidence, we found a modifi-
cation to enforce (i) and (ii). Unlike Boyd’s
algorithm

, it remains coordinate-free. We will

make use of the coloring introduced in the last sec-
tion. One can prove by induction that the spheres in
a Soddy set are always of distinct colors.

The modification basically consists in cutting

additional computation branches. The cases to be
treated specially are:

• If the newly generated sphere S

has larger

radius than the smallest sphere S

in the Soddy

set, cut the branch.

• If rad(S

) = rad(S

), do not count the new

sphere, but cut the branch only if color(S

) <

color(S

). This is to cut exactly one of two iden-

tical branches. The current Soddy set will be
produced in a second branch with the rôles of S

and S

swapped.

• If rad(S

) < rad(S

), do not cut the branch.

Count the new sphere in general. In order to find
the exceptions, perform a “look-ahead” Soddy
operation replacing S

and producing a new

sphere S

- If rad(S

) > rad(S

), do not count the sphere.

- If rad(S

) = rad(S

), count the new sphere

only as a half.

Although we do not have a rigorous correctness

proof, this 3-dimensional Soddy-algorithm has
been helpful for double-checking the results of the
inversion algorithm.

4. NUMERICAL RESULTS

We generated all spheres of curvature less than

calculating on a cluster of 20 DEC 3000-

400AXP workstations during 5.5 days. The spheres
have been used to sample the three functions N(

) ,

) and p(

) at the values

= 2

0.1

, 2

0.2

,..., 2

19.0

–

1 1 1 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

…

1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

1 1 1

–

0 0 0 0 1

2 1 1 1

–

1 2 1 1

–

1 1 2 1

–

1 1 1 2

–

1 1 1 1

–

∑













–

------------

1 j

≠

∑

Figure 3: Estimates for the fractal dimension as a function of the log

of the maximal curvature.

Table 1 displays them for some values of

. We

also confirm Boyd’s number (see Ref. 8) of
1 693 595 spheres having a curvature less than 601.

For each ten successive samples we did a least-

square-fit using the asymptotic Eq. (1) in logarith-
mic form. The obtained estimates for D which we
denote by D

(

), D

(

), D

(

) are plotted in Fig. 3.

It is interesting to observe that the fluctuations of
the latter two almost cancel each other. Their arith-

metic mean is also shown in Fig. 3 and yields the
estimate of 2.4739465 ± 0.0000001 for the fractal
dimension of the three-dimensional osculatory
packing.

Acknowledgement

We are grateful to Prof. R. Wuillemin of the EPF
Lausanne for generously granting the use of 20
new-installed workstations for a full week.

Table 1: Excerpt from the Computer Output

log

(

)

(

)

(

)

(

κ)

26.94641279

2.258953709774

1236

76.96194118

0.717030382959

6335400

396.34565660

0.116244349769

33532722660

2048.53328254

0.018774688805

31944875541924

7623.11018362

0.004366569649

D( )

log

2.473940

2.473941

2.473942

2.473943

2.473944

2.473945

2.473946

2.473947

2.473948

2.473949

2.473950

2.473951

2.473952

2.473953

2.473954

2.473955

2.473956

2.473957

Legend:

D ( )

(D ( )+D ( ))/2