HUTCHINSON OPERATORS IN ℝ^{3}

CARLA M. RIGGI

Abstract. In this project we study the mathematics behind Hutchinson Operators in ℝ^{3} with primary attention given to Iterated Function Systems. We have developed software necessary for defining and rendering the attractor of an IFS and studying its properties. To this end we present a contraction factor theorem for affine maps in ℝ^{n} along with corollaries for special cases in ℝ^{3}. We examine the relationships between IFS attractors in ℝ^{3} and their orthogonal projections onto the xy, yz, and xz planes; specifically the case where a 2d projection is itself the attractor of an IFS in ℝ^{2}. We also define the group of symmetry operations for the Sierpinski Triangle and the Sierpinski Tetrahedron. Our work is ongoing in the study of a family of 48^{4} IFSs related to the Sierpinski Right Tetrahedron and the 6 symmetry groups of their attractors. Our goal is to count and classify all distinct attractors in this family; a task that requires a solution to what remains an open question in this field as far as we know. Finally, the chaos3d module is written for Maple in order to compute the contraction factors of affine maps, define various Hutchinson Operators in ℝ^{3}, and render their attractors.

1. Introduction

In my undergraduate Chaos and Fractals course at the University of Scranton we studied two-dimensional fractal structures generated using a type of Hutchinson Operator known as an Iterated Function System (IFS). The Chaos module for Maple written by Ken Monks allowed us to render two-dimensional IFS attractors, but knowing that Hutchinson’s IFS theory holds for ℝ^{n}, we were naturally curious to see a fractal in 3d. By generalizing the definitions and theorems applicable to Hutchinson Operators in ℝ^{2} that were presented in the course, we can determine the formulas necessary for computing the affine maps used to describe fractals in three dimensions as well as a contraction factor for an n-dimensional affine map. These results allow us to write the Chaos3d Maple module for computing and displaying fractals in ℝ^{3}. This program will be a resource for future students in the Chaos and Fractals course and provides us with a tool for further investigation into the properties of 3d attractors having 2d fractal projections, and for studying the symmetry groups of IFS attractor families.

The paper is organized as follows. We first present the mathematical background and notation needed for our research. In the second section we present our main results. Here is found a contraction factor theorem for affine maps in ℝ^{n} followed by corollaries for special cases in ℝ^{3}. We then go on to investigate some conditions on affine maps in an IFS that guarantee its ℝ^{3} attractor to have a fractal projection into ℝ^{2}. Next we study the symmetry groups of two famous fractals and then investigate the classification of an IFS family based on the symmetry groups of their attractors. Peitgen, Jürgens, and Saupe [4] counted and classified all the distinct attractors produced by the family of 512 IFS’s related to the Sierpinski right triangle. In this paper we study a similar family in ℝ^{3}: the family of 48^{4} IFS’s related to the Sierpinski right tetrahedron. All proofs of theorems and corollaries stated in our main results are presented in the following section. The appendices contain additional notation, Maple programming documentation, and data resulting from our research. In Appendix A we define our notation for the symmetry group of the cube and list sets of elements useful in defining IFS’s having symmetric attractors including the center, and a few centralizer, and normalizer subgroups. In order to write the Chaos3d Maple module, it was first necessary to create the data structures with which to define affine maps in ℝ^{3}. These five affine data forms are defined in detail in Appendix B. Appendix C holds the Chaos3d source code itself and Appendix D contains help on using these procedures for uniting maps into Hutchinson operators and then rendering their attractors. Appendix E lists the IFS definitions for 8 distinct attractors related to the Sierpinski right tetrahedron that are each symmetric with respect to all six possible symmetry operations on this family of attractors. When using Maple to investigate properties of elements in our symmetry group of the cube, it is very convenient to maintain the notation we have defined. Therefore, an extension to Maple’s group package is written by Ken Monks specifically for this symmetry notation and these procedures are listed in Appendix F.

Date: April 15, 2001.

2000 Mathematics Subject Classification. 28A80.

Key words and phrases. fractals, IFS, Hutchinson operators, symmetry groups, computer graphics.

I would like to thank Dr. Monks without whose guidance and mathematical expertise I could not have completed this project.

2. Mathematical Background and Notation

A set theoretic discrete dynamical system is an ordered pair, (X, f), where X is a set and f : X → X. For any f : X → X we define f^{0} = i_{X}and f^{n} = f ∘ f^{n−1}, for all n ∈ ℕ^{+} where i_{X} is the identity map i_{X} : X → X given by i_{X}(s) = s, for all s ∈ X.

A metric space is a pair (X, d) where X is a set and d : X ×X → ℝ such that ∀x, y, z ∈ X, the following properties hold: d (x, y) ≥ 0, d (x, y) = 0 ⇔ x = y, d (x, y) = d (y, x), and d (x, y) + d (y, z) ≥ d (x, z). Here d is called a metric (or distance function) on X, and the elements of X are points in the metric space. An infinite sequence x_{0}, x_{1}, x_{2}, ... of points in X is a Cauchy sequence if and only if ∀ϵ ∈ ℝ^{+} there exists an N ∈ ℤ^{+} such that ∀i > N and ∀j > N, d (x_{i}, x_{j}) < ϵ. A sequence x_{0}, x_{1}, x_{2}, ... converges to x ∈ X if and only if ∀ϵ ∈ ℝ^{+} there exists an M ∈ ℕ such that ∀m > M, d (x_{m}, x) < ϵ (in this case we write lim_{m→∞} x_{m} = x). A metric space is complete if and only if every Cauchy sequence in (X, d) converges to a limit x ∈ X. One family of complete metric spaces is (ℝ^{n}, d_{Euc}) where d_{Euc} is the Euclidean metric given by d_{Euc} ((x_{1}, x_{2}, ..., x_{n}) , (y_{1}, y_{2}, ..., y_{n})) = \sqrt{∑_{i=1}^{n}(y_{i} − x_{i})^{2}}. However, many of the definitions and theorems in this paper hold for any complete metric space.

Given a metric space (X, d), f : X → X we define the scaling factor of f to be any value c ∈ ℝ^{+} such that ∀x, y ∈ X, d (f (x) , f (y)) ≦ cd (x, y) if such a value exists. We say f is a contraction mapping if and only if f has a scaling factor c ∈ (0, 1) . In this situation c is called a contraction factor of f. Given a set theoretic discrete dynamical system (X, f), the sequence f^{0}(x_{0}), f^{1}(x_{0}), f^{2}(x_{0}), ... of iterations of f acting on a seed value, x_{0} ∈ X, is called the f-orbit of x_{0}. For ease of notation, we will also sometimes write the f-orbit of x_{0} as the sequence, x_{0,}x_{1,}x_{2,}... where x_{n} = f (x_{n−1}) for all n ∈ ℕ^{+}. The Contraction Mapping Theorem states that for any contraction mapping f : X → X on a complete metric space (X, d) with contraction factor c, (1) f has a unique fixed point, q, (2) ∀ x ∈ X, lim_{n→∞} f^{n}(x) = q (i.e. the f-orbit of every x ∈ X converges to the fixed point q), and (3) if x_{0}, x_{1}, x_{2}, ... is the f-orbit of x_{0} ∈ X then d(x_{n}, q) ≤ \frac{c^{n}}{1−c}d(x_{0}, x_{1}).

Define K_{n} = {A ⊆ ℝ^{n} : A is compact}. The Hausdorff metric d_{H} on K_{n} is defined as follows. Let d_{H} : K_{n} × K_{n} → ℝ and S, T ∈ K_{n}. Define the closed collar of radius δ about S to be the set \overline{B} (S; δ) consisting of all points contained in the union of closed balls of radius δ about all α ∈ S. We define d_{H} (S, T) = inf \{δ : S ⊆ \overline{B} (T; δ) and T ⊆ \overline{B} (S; δ)\} . It is known that (K_{n}, d_{H}) is a complete metric space.

Let w_{1}, w_{2}, ..., w_{k} be contraction mappings on ℝ^{n} with contraction factors c_{1}, c_{2}, ..., c_{k} respectively. The Hutchinson operator, W : K_{n} → K_{n}, associated with w_{1}, w_{2}, ..., w_{k} is defined such that ∀A ∈ K_{n},

W (A) = w_{1} (A) ∪ w_{2} (A) ∪ ... ∪ w_{k} (A)  (1)

Hutchinson showed that such a map W is always a contraction map on (K_{n}, d_{H}) with contraction factor c = max{c_{1}, c_{2}, ..., c_{k}}. Thus by the contraction mapping theorem, W has a unique fixed point and the W -orbit of any compact set converges to this fixed point. The unique fixed point of W is called the attractor of W and is denoted F_{W}. We will often abbreviate an operator of the form (1) by writing W = w_{1} ∪ w_{2} ∪ ... ∪ w_{k}. We also will find it convenient to specify such an operator by its action on an arbitrary one point set, {v} , and further abbreviate W ({v}) = w_{1} ({v}) ∪ w_{2} ({v}) ∪ ... ∪ w_{k} ({v}) as W (v) = w_{1} (v) ∪ w_{2} (v) ∪ ... ∪ w_{k} (v) as usual in this situation.

From the above definitions we see that

F_{W} = W (F_{W}) = w_{1} (F_{W}) ∪ w_{2} (F_{W}) ∪ ... ∪ w_{k} (F_{W})

Therefore, F_{W} is the unique compact set in K_{n} that is composed entirely of contracted w_{i}-copies of itself. A shape satisfying such a geometric condition is often called a fractal. For all S ⊆ ℝ^{n}, define Sym (S) to be the symmetry group of S.

Define the sequence space on n letters, Σ_{n}, to be the set of infinite sequences composed of natural numbers in 1 thru n, i.e. Σ_{n} = {s : s = s_{1}s_{2}... and s_{1}, s_{2}, ... ∈ {1..n}} . Define

the function Φ_{W} : Σ_{n} → ℝ^{n} by Φ_{W} (s) = lim_{m→∞} w_{s_{1}} ∘ w_{s_{2}} ∘ ... ∘ w_{s_{m}} (\overrightarrow{0_{n}}) where \overrightarrow{0_{n}} =

(0, 0, ..., 0) ∈ ℝ^{n}. It is known that F_{W} = Φ_{W} (Σ_{n}) . So ∀a ∈ F_{W,} ∃r ∈ Σ_{n}, a = Φ_{W} (r) = lim_{m→∞} w_{r_{1}} ∘w_{r_{2}} ∘...∘w_{r_{m}} (\overrightarrow{0_{n}}) . Such a value r ∈ Σ_{n} is called an address of the point a ∈ F_{W} and we say that a is the point corresponding to the address r.

In this paper we focus on Hutchinson operators of a special type known as an Iterated Function System (IFS). An IFS is a Hutchinson operator composed entirely of affine contraction maps [3]. A function T : ℝ^{n} → ℝ^{n} is an affine map (or affine transformation) if and only if T (x) = Mx + B for some matrix M ∈ M_{n}(ℝ) and some B ∈ ℝ^{n}. Here B defines the translation of the origin under T, and the columns of M define the image of each basis vector in ℝ^{n} under T before translation.

3. Statement of the Main Results

3.1. Contraction factor for Affine Maps in ℝ^{n}. A very basic question to ask regarding any operator of the form (1) is whether it is an IFS. Since an IFS can only be formed using affine maps which are also contraction maps, it is necessary to be able to determine if a given affine map is a contraction map. In this situation, knowing a contraction factor for such maps would allow us to determine the rate of convergence of the IFS as well. The most accurate convergence rate for an affine map is given by its smallest contraction factor. This value is the best possible contraction factor for that map. For affine maps in ℝ^{2}, the use of complex analysis makes this computation relatively simple, but naturally the question arises: How can we compute the best possible contraction factor for affine maps in other dimensions? These problems are solved by the following theorem.

Theorem 3.1. (Contraction Factor) Let T : ℝ^{n} → ℝ^{n} by T (x) = Mx + B where M ∈ M_{n}(ℝ). Let c = max \{λ : λ is an eigenvalue of M^{T} M\}. T is a contraction map iff \sqrt{c} < 1. If T is a contraction map then \sqrt{c} is the best possible contraction factor.

As a special case we can derive specific formulas for the contraction factor of affine maps in ℝ^{3}. Corollary 3.2. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B, where M is a diagonal matrix. Then a scaling factor for T is

max {|m| : m is a value on the main diagonal of M} .

If max {|m|} < 1 then T is a contraction mapping and this is the best possible contraction factor. Corollary 3.3. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B where M^{T} M is diagonal. Then a scaling factor for T is

max \{\sqrt{α} : α is a value on the main diagonal of M^{T} M\} .

If max \{\sqrt{α}\} ∈ (0, 1) then T is a contraction mapping and this is the best possible contraction factor.

It is easy to show that M^{T} M is diagonal if and only if M is a matrix that rotates and scales the unit cube such that T of the unit cube is still an orthogonal parallelepiped. A matrix, M, of this type can be defined according to Geometric form given by

M = [ \begin{matrix}cos γ & − sin γ & 0 \\ sin γ & cos γ & 0 \\ 0 & 0 & 1 \end{matrix} ] [ \begin{matrix}cos ϕ & 0 & sin ϕ \\ 0 & 1 & 0 \\ − sin ϕ & 0 & cos ϕ \end{matrix} ] [ \begin{matrix}1 & 0 & 0 \\ 0 & cos θ & − sin θ \\ 0 & sin θ & cos θ \end{matrix} ] [ \begin{matrix}r & 0 & 0 \\ 0 & s & 0 \\ 0 & 0 & t \end{matrix} ] (2)

= [ \begin{matrix}r cos ϕ cos γ & s sin θ sin ϕ cos γ − s cos θ sin γ & t cos θ sin ϕ cos γ + t sin θ sin γ \\ r cos ϕ sin γ & s sin θ sin ϕ sin γ + s cos θ cos γ & t cos θ sin ϕ sin γ − t sin θ cos γ \\ −r sin ϕ & s sin θ cos ϕ & t cos θ cos ϕ \end{matrix} ]

This matrix first scales the i⃗, j⃗, and k⃗ vectors by real-valued constants r, s,and t respectively. The resulting orthogonal parallelepiped rotates about the x, y, and z axes consecutively by angles θ, ϕ,and γ respectively. When M is defined in this manner it is easy to determine the contraction values for its associated affine map according to the following corollary. Corollary 3.4. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B, where M is a matrix in the form (2) with r, s, t ∈ ℝ^{+}. Then a scaling factor for T is

max {r, s, t}

If max {r, s, t} < 1 then T is a contraction mapping and this is the best possible contraction factor.

3.2. 2d Projections of 3d IFS Attractors. We can analyze much of the structure of a given 3d Euclidean geometric shape based upon the properties of its orthogonal projections. It is natural then to ask what a projection into ℝ^{2} of a fractal in ℝ^{3} can tell us about the structure of the whole and vice versa. Specifically, we ask, when might a 2-dimensional projection of some IFS attractor in ℝ^{3} itself be the attractor of an IFS in ℝ^{2}? By examining the standard projections onto the xy, yz, and xz planes we find that a fractal projection onto one of these planes is guaranteed if every affine map in the IFS is in one of the forms stated in the following theorems.

Definition 3.5. Define P_{x} : ℝ^{3} → ℝ^{2}, P_{y} : ℝ^{3} → ℝ^{2}, and P_{z} : ℝ^{3} → ℝ^{2} by

P_{x} (x, y, z) = (y, z)
P_{y} (x, y, z) = (x, z)
P_{z} (x, y, z) = (x, y)

for all x, y, z ∈ ℝ. Theorem 3.6. Let w : ℝ^{3} → ℝ^{3} and v : ℝ^{2} → ℝ^{2} be affine maps. Then P_{z} ∘ w = v ∘ P_{z} if and only if there exists a, b, c, d, e, f, g, h, i ∈ ℝ such that

w (x, y, z) = [ \begin{matrix}a & b & 0 \\ c & d & 0 \\ e & f & g \end{matrix} ] [ \begin{matrix}x \\ y \\ z \end{matrix} ] + [ \begin{matrix}h \\ i \\ j \end{matrix} ]

and

v (x, y) = [ \begin{matrix}a & b \\ c & d \end{matrix} ] [ \begin{matrix}x \\ y \end{matrix} ] + [ \begin{matrix}h \\ i \end{matrix} ]

for all (x, y, z) ∈ ℝ^{3}. Theorem 3.7. Let W = w_{1} ∪ ⋯ ∪ w_{k} and V = v_{1} ∪ ⋯ ∪ v_{k} be IFS’s on ℝ^{3} and ℝ^{2} respectively such that P_{z} ∘ w_{i} = v_{i} ∘ P_{z} for all i ∈ {1, …, k} .Then

P_{z} ∘ Φ_{W} = Φ_{V} .

Theorem 3.8. Let W = w_{1} ∪ ⋯ ∪ w_{k} and V = v_{1} ∪ ⋯ ∪ v_{k} be IFS’s on ℝ^{3} and ℝ^{2} respectively such that P_{z} ∘ w_{i} = v_{i} ∘ P_{z} for all i ∈ {1, …, k} .Then

P_{z} (F_{W}) = F_{V} .

Theorems similar to each of the three above hold when P_{z} is replaced by P_{y} and when w in Theorem 3.6 is replaced by

w (x, y, z) = [ \begin{matrix}a & 0 & b \\ c & d & e \\ f & 0 & g \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}h \\ i \\ j \end{matrix} ] and v (x, z) = [ \begin{matrix}a & b \\ f & g \end{matrix} ] ( \begin{matrix}x \\ z \end{matrix} ) + [ \begin{matrix}h \\ j \end{matrix} ]

and again in the case where P_{z} is replaced by P_{x} and when w in Theorem 3.6 is replaced by

w (x, y, z) = [ \begin{matrix}a & b & c \\ 0 & d & e \\ 0 & f & g \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}h \\ i \\ j \end{matrix} ] and v (x, y) = [ \begin{matrix}d & e \\ f & g \end{matrix} ] ( \begin{matrix}x \\ y \end{matrix} ) + [ \begin{matrix}i \\ j \end{matrix} ] .

As a natural consequence of these theorems, we know that an IFS whose affine maps all

have matrices of the form [ \begin{matrix}a & 0 & 0 \\ 0 & b & 0 \\ 0 & 0 & c \end{matrix} ] for some a, b, c ∈ ℝ, will have a fractal projection

onto each of the three planes defined according to the following corollary. Corollary 3.9. Let W = w_{1} ∪ ⋯ ∪ w_{k} be an IFS on ℝ^{3} and let V = v_{1} ∪ ⋯ ∪ v_{k}, U = u_{1} ∪ ⋯ ∪ u_{k}, and T = t_{1} ∪ ⋯ ∪ t_{k} be IFS’s on ℝ^{2} such that for all i ∈ {1, …, k}

P_{z} ∘ w_{i} = v_{i} ∘ P_{z}
P_{x} ∘ w_{i} = u_{i} ∘ P_{x}
P_{y} ∘ w_{i} = t_{i} ∘ P_{y}

Then

P_{z} (F_{W}) = F_{V}
P_{x} (F_{W}) = F_{U}
P_{y} (F_{W}) = F_{T}

One such IFS produces the Sierpinski Tetrahedron, whose three standard projections in ℝ^{2} are each the Sierpinski Triangle. This IFS is given by

W (v) = Mv ∪ Mv + b_{1} ∪ Mv + b_{2} ∪ Mv + b_{3}  (3)

for all v ∈ ℝ^{3} where M = [ \begin{matrix}.5 & 0 & 0 \\ 0 & .5 & 0 \\ 0 & 0 & .5 \end{matrix} ], b_{1} = [ \begin{matrix}0.5 \\ 0 \\ 0 \end{matrix} ] , b_{2} = [ \begin{matrix}0 \\ 0.5 \\ 0 \end{matrix} ] , and b_{3} = [ \begin{matrix}0 \\ 0 \\ 0.5 \end{matrix} ] .

Its projection onto the xy plane is the Sierpinski Triangle, which is the attractor of

V (u) = M′u ∪ M′u + b′_{1} ∪ M′u + b′_{2}  (4)

for all u ∈ ℝ^{2} where M′ = [ \begin{matrix}.5 & 0 \\ 0 & .5 \end{matrix} ] , b′_{1} = [ \begin{matrix}0.5 \\ 0 \end{matrix} ] and b′_{2} = [ \begin{matrix}0 \\ 0.5 \end{matrix} ] .

Since the above theorems are not if and only if cases, it is natural to ask if there exist IFS’s on ℝ^{3} of a form other than that suggested by the previous theorems whose attractors would also have projections into ℝ^{2} that were themselves the attractors of 2d IFS’s. Furthermore we would like to know if it is possible to create a fractal in ℝ^{3} having any three twodimensional IFS attractors as its projections onto each one of the three coordinate planes. In the general case, the answer is clearly no, as can be seen by the following counterexample. Given that a single point can be the attractor of an IFS, we can ask if a 3-dimensional fractal exists which has a single point as its projection onto the xy, and yz planes and having the Sierpinski Triangle as its projection onto the xz plane. This clearly cannot be done as the only shape which can have a single point for two of its projections is a single point.

However, a typical attractor of an IFS is an uncountable set of points. But even if we restrict ourselves to such attractors for the choice of the three projections, the answer is still clearly no. It is clear that any object in ℝ^{3} must be contained completely within the intersection of the three sets formed by extending each projection into its respective third dimension. Similarly, a 3d IFS attractor would have to be entirely within the intersection of the x, y, or z extensions formed by each attractor in ℝ^{2}, but this intersection can be empty. Thus it is not possible in general to find an IFS on ℝ^{3} whose attractor has a given set of projections, even if those projections are the attractors of two-dimensional IFS’s.

3.3. Relatives of the Sierpinski Tetrahedron. There are many other IFS’s besides V given by (4) which produce the Sierpinski Triangle as their attractor. Peitgen, Jürgens, and Saupe [4] (pg 246) studied and classified the attractors of a family of IFS’s related to V. Similarly, we now investigate a family of IFS’s related to W given by (3).

In the case of the Sierpinski triangle, notice that each affine map in V has the effect of scaling each of the i⃗ and j⃗ basis vectors by 1/2 and then translating them to one of three fixed positions while keeping the vectors orthogonal to one another without any rotations or reflections. A graphical representation of the first iteration of this IFS using MrFace (Appendix C) as a seed shows the relation of each affine map to one another. Each face corresponds to a single affine map.

First iteration of V using MrFace as seed  6^{th} iteration of V using a solid square as seed  (5)

We can reorient any of these faces according to the 8 symmetry operations of the square without changing their relative positions to one another. Each such orientation and position corresponds to a unique affine map. Therefore it is possible to define 8^{3} = 512 distinct IFS’s just by using this Sierpinski-related L-shaped configuration of maps. Peitgen, Jürgens, and Saupe studied and classified the attractors of this IFS family. In particular, they found within this family alone 8 distinct IFS’s that each generate the Sierpinski Triangle as their attractor. Moreover, there are 7 other distinct attractors each of which is generated by 8 distinct IFS’s from this family.

Naturally we would like to know if analogous results hold in ℝ^{3}. Therefore we consider a family of IFS’s related to the Sierpinski Tetrahedron as produced by the IFS W given by (3). Rendering the first iteration of W using MrBlockHead (Appendix C) as our seed we can see a geometric representation of each affine map involved.

First iteration of W using MrBlockHead as seed  (6)

Just as the MrFaces above are all contained within the unit square, here all four blocks are similarly contained within the unit cube. Each of the four cubes corresponds to one of the four affine maps in the definition of W. Since there are 48 symmetry operations for the cube, there are 48^{4} different possible ways of orienting these four blocks within this basic tetragonal formation, with each configuration corresponding to a unique IFS. We describe this family of IFS’s precisely in the following definition.

Definition 3.10. Let v_{1}, v_{2}, v_{3}, and v_{4} be the affine maps on ℝ^{3} defined as follows

v_{1} (x, y, z) = (1/2x, 1/2y, 1/2z)
v_{2} (x, y, z) = v_{1} (x, y, z) + (1/2, 0, 0)
v_{3} (x, y, z) = v_{1} (x, y, z) + (0, 0, 1/2)
v_{4} (x, y, z) = v_{1} (x, y, z) + (0, 1/2, 0)

for all (x, y, z) ∈ ℝ^{3}. Let C denote the group of symmetry operations of the cube. Then we define the set of IFS relatives of the Sierpinski Tetrahedron to be the set

Ψ = {W : W = v_{1}d_{1} ∪ v_{2}d_{2} ∪ v_{3}d_{3} ∪ v_{4}d_{4} and d_{1}, d_{2}, d_{3}, d_{4} ∈ C}

Note that in this definition, we are identifying an element d of C with the unique affine map on ℝ^{3} whose restriction to the cube [0, 1]^{3} is the symmetry operation d. Our notation for the elements in the group of symmetry operations on the cube is defined in Appendix A. Six of these symmetry operations are needed in order to investigate important reflective and rotational symmetries according to which we can categorize the attractors of Ψ. We define this set of reflections and rotations as follows.

Δ = \{ \begin{matrix}e = the identity map \\ α_{1} = −t = reflection across the x = y plane \\ α_{2} = −r^{2}s^{3} = reflection across the y = z plane \\ α_{3} = −r^{3} = reflection across the x = z plane \\ α_{4} = rt = rotation about the line x = y = z by 120^{∘} directed from the x-axis to z-axis \\ α_{5} = r^{3}s^{3} = rotation about the line x = y = z by 120^{∘} directed from the x-axis to y-axis \end{matrix} \}

(7)

The following definition clarifies a subscripting notation that will be used frequently in the following theorems.

Definition 3.11. Let α ∈ Δ and ⋃_{i=1}^{4} v_{i}d_{i} ∈ Ψ. Let σ : Δ × {1, 2, 3, 4} → {1, 2, 3, 4} be a

map such that for each α, define σ_{α} : {1, 2, 3, 4} → {1, 2, 3, 4} such that σ_{α} (i) = σ (α, i) is given by

v_{σ(α,i)} = α^{−1}v_{i}α.

For all α, σ_{α} ∈ S_{4}. By direct calculation, one can verify that sigma is a well defined function. We will usually abbreviate σ (α, i) = \overline{i} when the value of α is clear from context. The values of \overline{i} are listed in the following table.

\begin{matrix}i & 1 & 2 & 3 & 4 \\ σ (e, i) & 1 & 2 & 3 & 4 \\ σ (α_{1}, i) & 1 & 4 & 3 & 2 \\ σ (α_{2}, i) & 1 & 2 & 4 & 3 \\ σ (α_{3}, i) & 1 & 3 & 2 & 4 \\ σ (α_{4}, i) & 1 & 4 & 2 & 3 \\ σ (α_{5}, i) & 1 & 3 & 4 & 2 \end{matrix}

Given any W ∈ Ψ we can find a corresponding W ′ ∈ Ψ whose attractor is a reflection or rotation of F_{W} with respect to some α ∈ Δ according to the following theorem.

Theorem 3.12. Let α ∈ Δ and W = ⋃_{i=1}^{4} v_{i}d_{i} ∈ Ψ. Define W′ = ⋃_{i=1}^{4} v_{i}d′_{i} ∈ Ψ by

d′_{i} = αd\overline{i}α^{−1} for all i ∈ \{1, 2, 3, 4\}

then α (F_{W} ) = F_{W ′}.

The Sierpinski tetrahedron and triangle are the prototypes for the families of IFS’s being considered. As a starting point for the study of symmetries of attractors in these families, we can determine its group of symmetry operations for these two fractals.

Theorem 3.13. Let (Sym (t) , ∘) be the group of all symmetry operations on the Sierpinski right triangle and let (Sym (T) , ∘) be the group of all symmetry operations on the Sierpinski right tetrahedron. Then

a) Sym (t) is isomorphic to ℤ_{2} (and hence D_{2})

b) Sym (T) is isomorphic toS_{3} (and hence D_{3})

As mentioned above, in the case of the Sierpinski triangle there are 8 distinct V -related (4) IFS’s that will produce it. So naturally we would like to know how many distinct W ∈ Ψ produce the Sierpinski right tetrahedron as their attractor. In general, we would like to be able to classify all the attractors of the set Ψ according to the following criteria.

Given W ∈ Ψ such that W (A) = A we wish to know

1) What is the order of the symmetry group of A? ( | Sym (A) | =?)

2) How many distinct W′ ∈ Ψ also produce A as their attractor?

We can answer the second question if we know the answer to the first as shown in the next theorem.

Theorem 3.14. Let W ∈ Ψ such that W (A) = A. Let Sym (A) be the group of symmetry operations on A.

Case I : If | Sym (A) | = 1 then only 1 IFS in Ψ has attractor A.
Case II : If | Sym (A) | = 2 then 2^{4} = 16 distinct IFS’s in Ψ have attractor A.
Case III : If | Sym (A) | = 3 then 3^{4} = 81 distinct IFS’s in Ψ have attractor A.
Case IV : If | Sym (A) | = 6 then 6^{4} = 1296 distinct IFS’s in Ψ have attractor A.

For example, by this theorem, we know that there must be 1296 distinct IFS’s in Ψ that all produce the Sierpinski right tetrahedron. Using our Maple program for rendering attractors, we have found 8 distinct attractors, A, of IFS’s in our family, Ψ, having | Sym (A) | = 6, each of which is the attractor of 1296 distinct IFS’s. The IFS definitions for each of these attractors are listed in Appendix E. The proof of this theorem, given in the following section, reveals how once we know the symmetry group belonging to the attractor of a given IFS W ∈ Ψ we can define all other IFS’s in our family that also converge to F_{W} . Moreover, if we know that W produces an asymmetric attractor, i.e. | Sym (F_{W}) | = 1, then Theorem 3.12 allows us to define the five other distinct W_{i} ∈ Ψ that each converge to α_{i} (F_{W}) for i ∈ {1..5}.

Clearly, in order to count the total number of distinct attractors in our family and to define the IFS’s that produce them, we require a condition on all d_{i} in a given W ∈ Ψ such that W (A) = A is guaranteed to have a symmetry group of a certain order. The following theorems make some progress toward solving this problem, but further work is required.

Theorem 3.15. Let α ∈ Δ and let W ∈ Ψ such that A = W (A) = ⋃_{i=1}^{4} v_{i}d_{i}(A).

α ∈ Sym (A) ⇔ ∀i, d_{i}^{−1}αd\overline{i} ∈ Sym (A) .

Theorem 3.16. Let α ∈ Δ and let W ∈ Ψ such that A = W (A) = ⋃_{i=1}^{4} v_{i}d_{i}(A).

If ∀i, αd\overline{i} = d_{i}α then α ∈ Sym (A) .

Theorem 3.17. Let α ∈ Δ and let W, W′ ∈ Ψ such that W = ⋃_{i=1}^{4} v_{i}d_{i}, W′ (A) =

⋃_{i=1}^{4} v_{i}di^{′} (A) and A = W′ (A) .

A = W (A) ⇔ ∀i, d_{i}^{−1}d′_{i} ∈ Sym (A)

We leave these as reference for those continuing to research the full classification of this IFS family and also in search of the answer to the broader question which asks how to determine the symmetry group of the attractor of a given IFS.

4. Proofs

4.1. Contraction Factor Theorem. The proof of the following contraction factor theorem relies upon the following definitions and theorems taken from Lay’s linear algebra book [2] where the proofs of Theorems 4.3 and 4.4 can be found. Definition 4.1. A quadratic form on ℝ^{n} is a function Q : ℝ^{n} → ℝ^{n} defined by Q(x⃗) = x^{T} Ax⃗ for all x⃗ ∈ ℝ^{n} where A ∈ M_{n}(ℝ) is symmetric (A^{T} = A). Definition 4.2. A quadratic form Q is

a) positive definite if Q(x⃗) > 0 for all x⃗ ≠ 0

b) negative definite if Q(x⃗) < 0 for all x⃗ ≠ 0

c) indefinite if Q(x⃗) assumes both positive and negative values Theorem 4.3 (A). A quadratic form x^{T} Ax⃗ is

a) positive definite iff the eigenvalues of A are all positive.

b) negative definite iff the eigenvalues of A are all negative.

c) indefinite iff A has both positive and negative eigenvalues. Theorem 4.4 (B). Let A be a symmetric matrix and define m = min {x^{T} Ax⃗ : ||x⃗|| = 1} and M = max {x^{T} Ax⃗ : ||x⃗|| = 1}. Then M is the largest eigenvalue, λ_{1}, of A and m is the smallest eigenvalue of A. The value of x^{T} x⃗ is M when x⃗ is a unit eigenvector, v⃗_{1}, corresponding to M.The value of x^{T} Ax⃗ is m when x⃗ is a unit eigenvector corresponding to m.

With this background we can now prove Theorem 3.1 which we restate here for convenience. Theorem 3.1. (Contraction Factor) Let T : ℝ^{n} → ℝ^{n}by T (x) = Mx + B where M ∈ M_{n}(ℝ). Let

c = max {λ : λ is an eigenvalue of M^{T} M}.

T is a contraction map iff \sqrt{c} < 1. If T is a contraction map then \sqrt{c} is the best possible contraction factor.

Proof. Let v⃗, w⃗ ∈ ℝ^{n} such that v⃗ ≠ w⃗. Let q⃗ = v⃗ − w⃗ and let u⃗ = \frac{q⃗}{||q⃗||}. The vector u⃗ is well defined since q ≠ 0. Furthermore, by definition of unit vector, ||u⃗|| = 1. Let r = ||q⃗||. Then by substitution, q⃗ = ru⃗. Let d be the Euclidean metric in.ℝ^{n}. Then by definition, d(v⃗, w⃗) = ||v⃗ − w⃗|| = ||q⃗|| = r. Now, let T : ℝ^{n} → ℝ^{n} such that T (x) = Mx + B⃗ where M ∈ M_{n}(ℝ) and B⃗ ∈ ℝ^{n}. Then

d (T \binom{}{w} , T \binom{}{v}) = ||T (v⃗) − T (w⃗)||
= || (Mv⃗ + B⃗) − (Mw⃗ + B⃗) || by definition of d
= ||Mv⃗ − Mw⃗||
= ||M(v⃗ − w⃗)||
= ||Mq⃗||
= ||M(ru⃗)||
= ||rMu⃗||
= |r| ||Mu⃗||
= r ||Mu⃗||
= r\sqrt{Mu⃗ · Mu⃗} by definition of vector length
= r\sqrt{(Mu⃗)^{T} Mu⃗} by the property of the dot product
= r\sqrt{u^{T} M^{T} Mu⃗} by property of the transpose

Let A = M^{T} M. It is well known that A is a symmetric matrix. Since Mu⃗ · Mu⃗ is always positive, u^{T} Au⃗ > 0,for all u⃗. Define c = max {u^{T} Au⃗ : ||u|| = 1} . Then by theorem 4.4, c = max {λ : λ is an eigenvalue of A} . Moreover, since u^{T} Au⃗ is positive definite, we know that c > 0 by Theorem 4.3 part a. According to the equations above, we have

d (T \binom{}{v} , T \binom{}{w})^{2} = r^{2}(u^{T} Au⃗)
≤ r^{2}c

and so since √ is an increasing function,

d (T \binom{}{v} , T \binom{}{w}) ≤ r\sqrt{c}
= \sqrt{c}d(v⃗, w⃗)

(<=) Assume \sqrt{c} < 1 and since we have shown that d (T \binom{}{v} , T \binom{}{w}) ≤ \sqrt{c}d(v⃗, w⃗) for all v⃗, w⃗ ∈ ℝ^{n}, we know by definition that T is a contraction map with contraction factor \sqrt{c}. (=>) To show that \sqrt{c} is the best possible contraction factor, assume that T is a contraction map with contraction factor s < 1. Let x⃗ be the unit eigenvector corresponding to eigenvalue c.

By Theorem 4.4 we know that x^{T} Ax⃗ = c. We therefore can show

d(T \binom{}{0} , T \binom{}{x})^{2} = ||T \binom{}{x} − T \binom{}{0}||^{2}
= ||Mx⃗ + B⃗ − (M0 + B⃗)||^{2}
= ||Mx⃗ − M0||^{2}
= ||M(x⃗ − 0)||^{2}
= ||Mx⃗||^{2}
= x^{T} M^{T} Mx⃗
= x^{T} Ax⃗
= c
= c(1^{2})
= c |\left|{ \atop x}\right||^{2}
= c ||x⃗ − 0||^{2}
= cd(0, x⃗)^{2}

and so

d(T (0), T (x⃗)) = \sqrt{c}d(0, x⃗)
≤ sd(0, x⃗)

which implies that \sqrt{c} ≤ s. By definition we know that s < 1 so \sqrt{c} < 1 and since\sqrt{c} ≤ s, we have shown that \sqrt{c} is the best possible contraction factor for T. □

Corollary 3.2. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B, where M is a diagonal matrix. Then a scaling factor for T is

max {|m| : m is a value on the main diagonal of M} .

If max {|m|} < 1 then T is a contraction mapping and this is the best possible contraction factor.

Proof. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B, where M is a diagonal matrix. From Theorem 3.1 we see that the scaling factor for any affine map T (x) = Mx+B is \sqrt{c} where c = max \{λ : λ is an eigenvalue of M^{T} M\} . If M is a diagonal matrix, then T has only the effect of scaling each of the basis vectors i⃗, j⃗, k⃗ by some constants x, y, z ∈ ℝ respectively before translation, and M is given by

M = [ \begin{matrix}x & 0 & 0 \\ 0 & y & 0 \\ 0 & 0 & z \end{matrix} ] .

Therefore,

M^{T} M = [ \begin{matrix}x & 0 & 0 \\ 0 & y & 0 \\ 0 & 0 & z \end{matrix} ] [ \begin{matrix}x & 0 & 0 \\ 0 & y & 0 \\ 0 & 0 & z \end{matrix} ] = [ \begin{matrix}x^{2} & 0 & 0 \\ 0 & y^{2} & 0 \\ 0 & 0 & z^{2} \end{matrix} ] .

Solving for the eigenvalues we get

0 = det(M^{T} M − λI) = det [ \begin{matrix}x^{2} − λ & 0 & 0 \\ 0 & y^{2} − λ & 0 \\ 0 & 0 & z^{2} − λ \end{matrix} ] = (x^{2} − λ)(y^{2} − λ)(z^{2} − λ)

Therefore, the eigenvalues of M^{T} M are x^{2}, y^{2}, and z^{2} and c = max {x^{2}, y^{2}, z^{2}}. Since the scaling factor is \sqrt{c} we know that for an affine map T having a diagonal matrix M the scaling factor is simply max {|m| : m is a value on the main diagonal of M}. Then by Theorem 3.1 and the definition of contraction map, if max {|m|} < 1 then T is a contraction mapping and max {|m|} is the best possible contraction factor. □

Corollary 3.3. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B where M^{T} M is diagonal. Then the scaling factor for T is

max \{\sqrt{α} : α is a value on the main diagonal of M^{T} M\} .

If max \{\sqrt{α}\} ∈ (0, 1) then T is a contraction mapping and this is the best possible contraction factor.

Proof. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B and M^{T} M is diagonal. From linear algebra we know that the eigenvalues of a diagonal matrix are equal to values on the main diagonal of that matrix. Therefore by Theorem 3.1 we have

c = max {λ : λ is an eigenvalue of M^{T} M}
= max {α : α is a value on the main diagonal of M^{T} M} .

Thus the scaling factor for T is \sqrt{c} = max \{\sqrt{α} : α is a value on the main diagonal of M^{T} M\} . Thus by definition of contraction map, if max \{\sqrt{α}\} ∈ (0, 1) then T is a contraction mapping and max \{\sqrt{α}\} is the best possible contraction factor. □

Corollary 3.4. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B, where M is a matrix in the form (2) with r, s, t ∈ ℝ^{+}. Then a scaling factor for T is

max {r, s, t}

If max {r, s, t} < 1 then T is a contraction mapping and this is the best possible contraction factor.

Proof. Let T : ℝ^{3} → ℝ^{3} be the affine map given by T (x) = Mx + B. Consider an affine map having a matrix M = (v⃗_{i}, v⃗_{j}, v⃗_{k}) composed of column vectors

v⃗_{i} = [ \begin{matrix}a_{i} \\ b_{i} \\ c_{i} \end{matrix} ] , v⃗_{j} = [ \begin{matrix}a_{j} \\ b_{j} \\ c_{j} \end{matrix} ] , v⃗_{k} = [ \begin{matrix}a_{k} \\ b_{k} \\ c_{k} \end{matrix} ]

that are images of i⃗, j⃗, k⃗ respectively under the affine map T before translation. We have

M^{T} M = [ \begin{matrix}a_{i} & b_{i} & c_{i} \\ a_{j} & b_{j} & c_{j} \\ a_{k} & b_{k} & c_{k} \end{matrix} ] [ \begin{matrix}a_{i} & a_{j} & a_{k} \\ b_{i} & b_{j} & b_{k} \\ c_{i} & c_{j} & c_{k} \end{matrix} ]

is diagonal if and only if

0 = a_{i}a_{j} + b_{i}b_{j} + c_{i}c_{j} = v⃗_{i} · v⃗_{j}
0 = a_{i}a_{k} + b_{i}b_{k} + c_{i}c_{k} = v⃗_{i} · v⃗_{k}
0 = a_{j}a_{k} + b_{j}b_{k} + c_{j}c_{k} = v⃗_{j} · v⃗_{k}.

Since the dot product of any column of M with the each of the other two columns must be zero, then the columns of M must be orthogonal and linearly independent. Therefore v⃗_{i}, v⃗_{j},and v⃗_{k} form an orthogonal basis for ℝ^{3}. Geometrically, since v⃗_{i}, v⃗_{j},and v⃗_{k} must all be orthogonal to one another we see that M^{T} M is diagonal if and only if M is a matrix that rotates and scales the unit cube such that T of the unit cube is still an orthogonal parallelepiped. Certainly then, an affine map determined using the Geometric form ( 2) for M,

M = [ \begin{matrix}r cos ϕ cos γ & s sin θ sin ϕ cos γ − s cos θ sin γ & t cos θ sin ϕ cos γ + t sin θ sin γ \\ r cos ϕ sin γ & s sin θ sin ϕ sin γ + s cos θ cos γ & t cos θ sin ϕ sin γ − t sin θ cos γ \\ −r sin ϕ & s sin θ cos ϕ & t cos θ cos ϕ \end{matrix} ]

will always have a c value, as in Theorem 3.1, given by

c = max {α : α is a value on the main diagonal of M^{T} M}
= max {r^{2}, s^{2}, t^{2}} .

Here, r, s, t ∈ ℝ^{+} are the scaling coefficients for i⃗, j⃗,and k⃗ respectively and so cannot be negative. Therefore, the scaling factor for T is

\sqrt{c} = max \{r, s, t\} .

Thus by definition of contraction map, if max {r, s, t} < 1 then T is a contraction mapping and max {r, s, t} is the best possible contraction factor. □

4.2. 2d Projections of 3d Attractors. For the definition of P_{z} see definition 3.5 above. Theorem 3.6. Let w : ℝ^{3} → ℝ^{3} and v : ℝ^{2} → ℝ^{2} be affine maps. Then P_{z} ∘ w = v ∘ P_{z} if and only if there exists a, b, c, d, e, f, g, h, i ∈ ℝ such that

w (x, y, z) = [ \begin{matrix}a & b & 0 \\ c & d & 0 \\ e & f & g \end{matrix} ] [ \begin{matrix}x \\ y \\ z \end{matrix} ] + [ \begin{matrix}h \\ i \\ j \end{matrix} ]

and

v (x, y) = [ \begin{matrix}a & b \\ c & d \end{matrix} ] [ \begin{matrix}x \\ y \end{matrix} ] + [ \begin{matrix}h \\ i \end{matrix} ]

for all (x, y, z) ∈ ℝ^{3}.

Proof. (=>)

Let w : ℝ^{3} → ℝ^{3} and v : ℝ^{2} → ℝ^{2} be affine maps such that

w (x, y, z) = [ \begin{matrix}w_{1} & w_{2} & w_{3} \\ w_{4} & w_{5} & w_{6} \\ w_{7} & w_{8} & w_{9} \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}w_{10} \\ w_{11} \\ w_{12} \end{matrix} ] and

v (x, y) = [ \begin{matrix}v_{1} & v_{2} \\ v_{3} & v_{4} \end{matrix} ] ( \begin{matrix}x \\ y \end{matrix} ) + [ \begin{matrix}v_{5} \\ v_{6} \end{matrix} ]

for some w_{i}, v_{i} ∈ ℝ.

Assume P_{z} ∘ w (x, y, z) = v ∘ P_{z} (x, y, z) for all (x, y, z) ∈ ℝ^{3}. Then

(w_{1}x + w_{2}y + w_{3}z + w_{10}, w_{4}x + w_{5}y + w_{6}z + w_{11})
= P_{z} ∘ (w_{1}x + w_{2}y + w_{3}z + w_{10}, w_{4}x + w_{5}y + w_{6}z + w_{11}, w_{7}x + w_{8}y + w_{9}z + w_{12})
= P_{z} ∘ w (x, y, z)
= v ∘ P_{z} (x, y, z)
= v (x, y)
= (v_{1}x + v_{2}y + v_{5}, v_{3}x + v_{4}y + v_{6})

By definition of equality for ordered pairs we have the following two equations,

w_{1}x + w_{2}y + w_{3}z + w_{10} = v_{1}x + v_{2}y + v_{5} and w_{4}x + w_{5}y + w_{6}z + w_{11} = v_{3}x + v_{4}y + v_{6}.

Since these equations must hold for all x, y, z ∈ ℝ^{3}, the first equation guarantees that

w_{1} = v_{1}
w_{2} = v_{2}
w_{3} = 0
w_{10} = v_{5}

and by the second equation we have

w_{4} = v_{3}
w_{5} = v_{4}
w_{6} = 0
w_{11} = v_{6}.

Therefore, w and v are of the form

w (x, y, z) = [ \begin{matrix}w_{1} & w_{2} & 0 \\ w_{4} & w_{5} & 0 \\ w_{7} & w_{8} & w_{9} \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}w_{10} \\ w_{11} \\ w_{12} \end{matrix} ] and

v (x, y) = [ \begin{matrix}w_{1} & w_{2} \\ w_{4} & w_{5} \end{matrix} ] ( \begin{matrix}x \\ y \end{matrix} ) + [ \begin{matrix}w_{10} \\ w_{11} \end{matrix} ] .

(<=)

Assume that w and v are defined such that

w (x, y, z) = [ \begin{matrix}w_{1} & w_{2} & 0 \\ w_{4} & w_{5} & 0 \\ w_{7} & w_{8} & w_{9} \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}w_{10} \\ w_{11} \\ w_{12} \end{matrix} ] and

v (x, y) = [ \begin{matrix}w_{1} & w_{2} \\ w_{4} & w_{5} \end{matrix} ] ( \begin{matrix}x \\ y \end{matrix} ) + [ \begin{matrix}w_{10} \\ w_{11} \end{matrix} ]

for some w_{i} ∈ ℝ. By the definition of P_{z} we have

P_{z} ∘ w (x, y, z) = P_{z} ∘ (w_{1}x + w_{2}y + w_{10}, w_{4}x + w_{5}y + w_{11}, w_{7}x + w_{8}y + w_{9}z + w_{12})
= (w_{1}x + w_{2}y + w_{10}, w_{4}x + w_{5}y + w_{11})
= v (x, y)
= v ∘ P_{z} (x, y, z)

thus P_{z} ∘ w (x, y, z) = v ∘ P_{z} (x, y, z). □

Theorem 3.7. Let W = w_{1} ∪ ⋯ ∪ w_{k} and V = v_{1} ∪ ⋯ ∪ v_{k} be IFS’s on ℝ^{3} and ℝ^{2} respectively such that P_{z} ∘ w_{i} = v_{i} ∘ P_{z} for all i ∈ {1, …, k} .Then

P_{z} ∘ Φ_{W} = Φ_{V} .

Proof. Let W = w_{1} ∪ ⋯ ∪ w_{k} and V = v_{1} ∪ ⋯ ∪ v_{k} be IFS’s on ℝ^{3} and ℝ^{2} respectively such that P_{z} ∘w_{i} = v_{i}∘P_{z} for all i ∈ {1, …, k} . Define the function Φ as in Section 2 but with the restrictions that Φ_{W} : Σ_{k} → ℝ^{3} and Φ_{V} : Σ_{k} → ℝ^{2}. Then P_{z} ∘ Φ_{W} : Σ_{k} → ℝ^{2}. Let s ∈ Σ_{k}, then s = s_{1}s_{2}... such that for all j, s_{j} ∈ {1, …, k} . By the definition of Φ_{W} we know that ∃a ∈ ℝ^{3} such that a = Φ_{W} (s) = lim_{n→∞} w_{s_{1}} ∘w_{s_{2}} ∘...∘w_{s_{n}} (\overrightarrow{0_{3}}) . Note that P_{z} is a continuous map. Since limits commute with continuous maps, lim_{i→∞} P_{z} (q_{i}) = P_{z} (lim_{i→∞} q_{i}) for any convergent sequence of points q_{i} ∈ ℝ^{3}. Therefore,

P_{z} ∘ Φ_{W} (s) = P_{z} ( lim_{n→∞} w_{s_{1}} ∘ w_{s_{2}} ∘ ... ∘ w_{s_{n}} (\overrightarrow{0_{3}}))
= lim_{n→∞} P_{z} (w_{s_{1}} ∘ w_{s_{2}} ∘ ... ∘ w_{s_{n}} (\overrightarrow{0_{3}}))
= lim_{n→∞} P_{z} ∘ w_{s_{1}} ∘ w_{s_{2}} ∘ ... ∘ w_{s_{n}} (\overrightarrow{0_{3}})
= lim_{n→∞} v_{s_{1}} ∘ P_{z} ∘ w_{s_{2}} ∘ ... ∘ w_{s_{n}} (\overrightarrow{0_{3}})
⋮
= lim_{n→∞} v_{s_{1}} ∘ v_{s_{2}} ∘ ... ∘ v_{s_{n}} ∘ P_{z} (\overrightarrow{0_{3}})
= lim_{n→∞} v_{s_{1}} ∘ v_{s_{2}} ∘ ... ∘ v_{s_{n}} (\overrightarrow{0_{2}})
= Φ_{V} (s) .

Thus P_{z} ∘ Φ_{W} (s) = Φ_{V} (s) for all s = s_{1}s_{2}... ∈ Σ_{k}. □

Theorem 3.8. Let W = w_{1} ∪ ⋯ ∪ w_{k} and V = v_{1} ∪ ⋯ ∪ v_{k} be IFS’s on ℝ^{3} and ℝ^{2} respectively such that P_{z} ∘ w_{i} = v_{i} ∘ P_{z} for all i ∈ {1, …, k} .Then

P_{z} (F_{W}) = F_{V} .

Proof. Let W = w_{1} ∪ ⋯ ∪ w_{k} and V = v_{1} ∪ ⋯ ∪ v_{k} be IFS’s on ℝ^{3} and ℝ^{2} respectively such that P_{z} ∘ w_{i} = v_{i} ∘ P_{z} for all i ∈ {1, …, k} . Let q ∈ P_{z} (F_{W}) . Then there must exist an a ∈ F_{W} such that q = P_{z} (a) . Since Φ_{W} (Σ_{k}) = F_{W} there exists s = s_{1}s_{2}... ∈ Σ_{k} such that a = Φ_{W} (s) . Therefore, by the previous theorem

q = P_{z} ∘ Φ_{W} (s)
= Φ_{V} (s) .

But Φ_{V} (s) ∈ Φ_{V} (Σ_{k}) = F_{V} . Thus q ∈ F_{V} and so P_{z} (F_{W}) ⊆ F_{V} . Now let q′ ∈ F_{V} . As above, we know that there must exist an r = r_{1}r_{2}... ∈ Σ_{k} such that q′ = Φ_{V} (r) . Thus

q′ = Φ_{V} (r)
= P_{z} ∘ Φ_{W} (r) .

Again we know, Φ_{W} (r) ∈ Φ_{W} (Σ_{k}) = F_{W} . So q′ ∈ P_{z} (F_{W}) and so F_{V} ⊆ P_{z} (F_{W}) . Thus P_{z} (F_{W}) = F_{V} . □

Similar proofs confirm these three theorems with respect to the alternate definitions for w, v,and P mentioned in this section above.

4.3. Sierpinski Symmetry groups. The definition for \overline{i} can be found in definition 3.11 above. For the following proof, note that αv\overline{i} = v_{i}α.

Theorem 3.12. Let α ∈ Δ and W = ⋃_{i=1}^{4} v_{i}d_{i} ∈ Ψ. Define W′ = ⋃_{i=1}^{4} v_{i}d′_{i} ∈ Ψ by

d′_{i} = αd\overline{i}α^{−1} for all i ∈ \{1, 2, 3, 4\} .

Then α (F_{W} ) = F_{W ′}.

Proof. Let α ∈ Δ and W = ⋃_{i=1}^{4} v_{i}d_{i} ∈ Ψ. F_{W} is the unique A ∈ K_{3} such that A = W (A).

Define W′ = ⋃_{i=1}^{4} v_{i}d′_{i} ∈ Ψ by d′_{i} = αd\overline{i}α^{−1} for all i ∈ \{1, 2, 3, 4\} .

αA = α (W (A))
= α (⋃_{i=1}^{4} v_{i}d_{i} (A))

= α (⋃_{i=1}^{4} v\overline{i}d\overline{i}A) by definition of \overline{i}

= ⋃_{i=1}^{4} αv\overline{i}d\overline{i}A

= ⋃_{i=1}^{4} αv\overline{i}d\overline{i}α^{−1}αA

= ⋃_{i=1}^{4} v_{i}αd\overline{i}α^{−1}αA

= ⋃_{i=1}^{4} v_{i}d′_{i} (αA)

So by definition and uniqueness of the attractor of an IFS we know that αA = W′ (αA) , and so α (F_{W}) = F_{W′}. □

Theorem 3.13. Let (Sym (t) , ∘) be the group of all symmetry operations on the Sierpinski right triangle and let (Sym (T) , ∘) be the group of all symmetry operations on the Sierpinski right tetrahedron. Then

a) Sym (t) is isomorphic to ℤ_{2} (and hence D_{2})

b) Sym (T) is isomorphic toS_{3} (and hence D_{3})

Proof. Part a)

On the Sierpinski Triangle, the only two points separated by a distance of \sqrt{2} are (1, 0) and (0, 1). Since a symmetry operation must preserve distances, these two points must map to each other. Furthermore, the only point contained in the Sierpinski right triangle that is exactly one unit away from each of these points is the origin. It therefore must map to itself. Since a symmetry map in ℝ^{2} is completely determined by where it sends three non-colinear points, the only possible symmetry operation on the Sierpinski Triangle other than the identity map is g : ℝ^{2} → ℝ^{2} such that g (x, y) = (y, x) for all x, y ∈ ℝ.

To see that g is indeed a symmetry operation on the Sierpinski right triangle, let W : ℝ^{2} → ℝ^{2} be the IFS whose attractor is the Sierpinski right triangle defined by V = v_{1} ∪ v_{2} ∪ v_{3} where v_{i} is defined as in (4) above. Using these definitions we can see that for any (x, y) ∈ ℝ^{2}

v_{1} (g (x, y)) = v_{1} (y, x) = (\frac{1}{2}y, \frac{1}{2}x) = g(\frac{1}{2}x, \frac{1}{2}y) = g (v_{1} (x, y))
v_{2} (g (x, y)) = v_{2} (y, x) = (\frac{1}{2}y + \frac{1}{2}, \frac{1}{2}x) = g(\frac{1}{2}x, \frac{1}{2}y + \frac{1}{2}) = g (v_{3} (x, y))
v_{3} (g (x, y)) = v_{1} (y, x) = (\frac{1}{2}y, \frac{1}{2}x + \frac{1}{2}) = g(\frac{1}{2}x + \frac{1}{2}, \frac{1}{2}y) = g (v_{2} (x, y)) .

And so we have

v_{1}g = gv_{1}
v_{2}g = gv_{3}
v_{3}g = gv_{2}.

Now let A = v_{1} (A) ∪ v_{2} (A) ∪ v_{3} (A) be the Sierpinski right triangle. From the above equations we see that

V (g (A)) = v_{1} (g (A)) ∪ v_{2} (g (A)) ∪ v_{3} (g (A))
= g (v_{1} (A)) ∪ g (v_{3} (A)) ∪ g (v_{2} (A))
= g (v_{1} (A) ∪ v_{2} (A) ∪ v_{3} (A))
= g (A) .

By uniqueness of the attractor of an IFS we can conclude that g (A) = A. So g is a symmetry operation on A and therefore, all possible symmetry operations on the Sierpinski right triangle are

g^{0} (x, y) = (x, y) = e
and g (x, y) = (y, z) .

Let Sym (t) = {e, g} then Sym (t) is isomorphic to ℤ_{2} (which is in turn isomorphic to D_{2}). Part b) On the Sierpinski right tetrahedron, the only three points separated by a distance of \sqrt{2}

are (1, 0, 0) and (0, 1, 0) and (0, 0, 1). Since a symmetry operation must preserve distances, these three points must map to one another. Furthermore, the only point contained in the Sierpinski right tetrahedron that is exactly one unit away from each of these points is the origin. It therefore must map to itself. Since a symmetry map in ℝ^{3} is completely determined by where it sends four non-coplanar points, the five possible symmetry operations on the Sierpinski right tetrahedron in addition to the identity map are defined by R_{i} : ℝ^{3} → ℝ^{3} for i ∈ {1, 2, 3, 4, 5} such that for all x, y, z ∈ ℝ

R_{1} (x, y, z) = (y, z, x)
R_{2} (x, y, z) = (z, x, y)
R_{3} (x, y, z) = (y, x, z)
R_{4} (x, y, z) = (x, z, y)
R_{5} (x, y, z) = (z, y, x)

We can show that all R_{i} can be generated by two of these functions. Therefore, rename R_{1} and R_{3} such that Q = R_{1} and L = R_{3}. Let W : ℝ^{3} → ℝ^{3} be the IFS producing the Sierpinski right tetrahedron defined by W = w_{1} ∪ w_{2} ∪ w_{3} ∪ w_{4} where w_{i} is defined as in (3) above. From the definitions for Q and L we see that

w_{1} (Q (x, y, z)) = w_{1} (y, z, x) = (\frac{1}{2}y, \frac{1}{2}z, \frac{1}{2}x) = Q (\frac{1}{2}x, \frac{1}{2}y, \frac{1}{2}z) = Q (w_{1} (x, y, z))
w_{2} (Q (x, y, z)) = w_{2} (y, z, x) = (\frac{1}{2}y + \frac{1}{2}, \frac{1}{2}z, \frac{1}{2}x) = Q (\frac{1}{2}x, \frac{1}{2}y + \frac{1}{2}, \frac{1}{2}z) = Q (w_{4} (x, y, z))
w_{3} (Q (x, y, z)) = w_{3} (y, z, x) = (\frac{1}{2}y, \frac{1}{2}z, \frac{1}{2}x + \frac{1}{2}) = Q (\frac{1}{2}x + \frac{1}{2}, \frac{1}{2}y, \frac{1}{2}z) = Q (w_{2} (x, y, z))
w_{4} (Q (x, y, z)) = w_{4} (y, z, x) = (\frac{1}{2}y, \frac{1}{2}z + \frac{1}{2}, \frac{1}{2}x) = Q (\frac{1}{2}x, \frac{1}{2}y, \frac{1}{2}z + \frac{1}{2}) = Q (w_{3} (x, y, z))

and

w_{1} (L (x, y, z)) = w_{1} (y, x, z) = (\frac{1}{2}y, \frac{1}{2}x, \frac{1}{2}z) = L (\frac{1}{2}x, \frac{1}{2}y, \frac{1}{2}z) = L (w_{1} (x, y, z))
w_{2} (L (x, y, z)) = w_{2} (y, x, z) = (\frac{1}{2}y + \frac{1}{2}, \frac{1}{2}x, \frac{1}{2}z) = L (\frac{1}{2}x, \frac{1}{2}y + \frac{1}{2}, \frac{1}{2}z) = L (w_{4} (x, y, z))
w_{3} (L (x, y, z)) = w_{3} (y, x, z) = (\frac{1}{2}y, \frac{1}{2}x, \frac{1}{2}z + \frac{1}{2}) = L (\frac{1}{2}x, \frac{1}{2}y, \frac{1}{2}z + \frac{1}{2}) = L (w_{3} (x, y, z))
w_{4} (L (x, y, z)) = w_{4} (y, x, z) = (\frac{1}{2}y, \frac{1}{2}x + \frac{1}{2}, \frac{1}{2}z) = L (\frac{1}{2}x + \frac{1}{2}, \frac{1}{2}y, \frac{1}{2}z) = L (w_{2} (x, y, z)) .

So now we have

\begin{matrix}w_{1}Q = Qw_{1} \\ w_{2}Q = Qw_{4} \\ w_{3}Q = Qw_{2} \\ w_{4}Q = Qw_{3} \end{matrix} and \begin{matrix}w_{1}L = Lw_{1} \\ w_{2}L = Lw_{4} \\ w_{3}L = Lw_{3} \\ w_{4}L = Lw_{2} \end{matrix} .

Let A = w_{1} (A) ∪ w_{2} (A) ∪ w_{3} (A) ∪ w_{4} (A) be the Sierpinski right tetrahedron. From the above equations we see that

W (Q (A)) = w_{1} (Q (A)) ∪ w_{2} (Q (A)) ∪ w_{3} (Q (A)) ∪ w_{4} (Q (A))
= Q (w_{1} (A)) ∪ Q (w_{4} (A)) ∪ Q (w_{2} (A)) ∪ Q(w_{3} (A))
= Q (w_{1} (A) ∪ w_{2} (A) ∪ w_{3} (A) ∪ w_{4} (A))
= Q (A)

and

W (L (A)) = w_{1} (L (A)) ∪ w_{2} (L (A)) ∪ w_{3} (L (A)) ∪ w_{4} (L (A))
= L (w_{1} (A)) ∪ L (w_{4} (A)) ∪ L (w_{3} (A)) ∪ L(w_{2} (A))
= L (w_{1} (A) ∪ w_{2} (A) ∪ w_{3} (A) ∪ w_{4} (A))
= L (A)

By uniqueness of the attractor of an IFS we can conclude that Q (A) = L (A) = A. So Q and L are symmetry operations on A. Furthermore, we can see that

R_{1} (x, y, z) = (y, z, x) = Q (x, y, z)
R_{2} (x, y, z) = (z, x, y) = Q^{2} (x, y, z)
R_{3} (x, y, z) = (y, x, z) = L (x, y, z)
R_{4} (x, y, z) = (x, z, y) = QL (x, y, x)
R_{5} (x, y, z) = (z, y, x) = Q^{2}L (x, y, z)

Therefore, all R_{i} are indeed symmetry operations on the Sierpinski right tetrahedron.

Let Sym(T) = {e, Q, Q^{2}, L, QL, Q^{2}L}. Then Sym (T) is isomorphic to S_{3} which is in turn isomorphic to D_{3}. □

4.4. Classification of Attractors According to Their Symmetry Group. The following theorems were devised and proven by Ken Monks. The definition for \overline{i} can be found in definition 3.11 above. Theorem 3.14. Let W ∈ Ψ such that W (A) = A. Let Sym (A) be the group of symmetry operations on A.

Case I : If | Sym (A) | = 1 then only 1 IFS in Ψ can produce A.
Case II : If | Sym (A) | = 2 then 2^{4} = 16 distinct IFS’s in Ψ can produce A.
Case III : If | Sym (A) | = 3 then 3^{4} = 81 distinct IFS’s in Ψ can produce A.
Case IV : If | Sym (A) | = 6 then 6^{4} = 1296 distinct IFS’s in Ψ can produce A.

Proof. Let W, W′ ∈ Ψ and W (A) = A. Assume W′ (A) = A. By definition of Ψ, W = v_{1}d_{1} ∪ v_{2}d_{2} ∪ v_{3}d_{3} ∪ v_{4}d_{4} and W′ = v_{1}d′_{1} ∪ v_{2}d′_{2} ∪ v_{3}d′_{3} ∪ v_{4}d′_{4}. By assumption,

v_{1}d_{1} (A) ∪ v_{2}d_{2} (A) ∪ v_{3}d_{3} (A) ∪ v_{4}d_{4} (A) = A = v_{1}d′_{1} (A) ∪ v_{2}d′_{2} (A) ∪ v_{3}d′_{3} (A) ∪ v_{4}d′_{4} (A)

Since the entire attractor must be contained within each v_{i}, we can reduce the problem to the separate equations v_{1}d_{1} (A) = v_{1}d′_{1} (A) , v_{2}d_{2} (A) = v_{2}d′_{2} (A) , v_{3}d_{3} (A) = v_{3}d′_{3} (A) , and v_{4}d_{4} (A) = v_{4}d′_{4} (A) . Furthermore, each v_{i} is invertible and therefore d_{1} (A) = d′_{1} (A) , d_{2} (A) = d′_{2} (A) , d_{3} (A) = d′_{3} (A) ,and d_{4} (A) = d′_{4} (A) . Thus

A = d_{1}^{−1}d′_{1} (A) , A = d_{2}^{−1}d′_{2} (A) , A = d_{3}^{−1}d′_{3} (A) and A = d_{4}^{−1}d′_{4} (A)

In order for this set of equations to be true, it must be the case that ∀i ∈ 1..4, d_{i}^{−1}d′_{i} ∈ Sym (A) . Since we only need to count the possible number of different combinations of d′_{i} and not their specific values, we have 4 cases to examine.

Case 1: | Sym (A) | = 1

Sym (A) = {e} . Therefore d_{1}^{−1}d′_{1} = d_{2}^{−1}d′_{2} = d_{3}^{−1}d′_{3} = d_{4}^{−1}d′_{4} = e and so ∀i, d′_{i} = d_{i}. Thus there is only one IFS that produces each asymmetric attractor.

Case 2: | Sym (A) | = 2

Sym (A) = {e, α} where α ∈ {α_{1}, α_{2}, α_{3}} . Therefore d_{1}^{−1}d′_{1} ∈ {e, α} , d_{2}^{−1}d′_{2} ∈ {e, α}, d_{3}^{−1}d′_{3} ∈ {e, α} and d_{4}^{−1}d′_{4} ∈ {e, α}. Thus ∀i, d′_{i} = d_{i} or d′_{i} = αd_{i}. And so there are 2^{4} = 16 different IFS’s that produce each attractor that is symmetric with respect to reflection across only one of the 3 diagonal planes.

Case 3: | Sym (A) | = 3

Sym (A) = {e, α_{4}, α_{5}} . Therefore as above ∀i, d′_{i} = d_{i} or d′_{i} = α_{4}d_{i} or d′_{i} = α_{5}d_{i}. Thus there are 3^{4} = 81 different IFS’s that produce each attractor that is symmetric only with respect to rotation about x = y = z.

Case 4: | Sym (A) | = 6

Sym (A) = {e, α_{1}, α_{2}, α_{3}, α_{4}, α_{5}}. As above ∀i, d′_{i} = d_{i} or d′_{i} = α_{1}d_{i} or d′_{i} = α_{2}d_{i} or d′_{i} = α_{3}d_{i} or d′_{i} = α_{4}d_{i} or d′_{i} = α_{5}d_{i}. Therefore there are 6^{4} = 1296 different IFS’s that produce each attractor having all six symmetries. □

Theorem 3.15. Let α ∈ Δ and let W ∈ Ψ such that A = W (A) = ⋃_{i=1}^{4} v_{i}d_{i} (A).

α ∈ Sym (A) ⇔ ∀i, d_{i}^{−1}αd\overline{i} ∈ Sym (A)

Proof. Let α ∈ Δ and let W ∈ Ψ such that A = W (A) = ⋃_{i=1}^{4} v_{i}d_{i} (A).

(=>) Assume α ∈ Sym (A) .

α (A) = A
α (W (A)) = W (A)
α (⋃_{i=1}^{4} v_{i}d_{i}A) = ⋃_{i=1}^{4} v_{i}d_{i}A

α (⋃_{i=1}^{4} v\overline{i}d\overline{i}A) = ⋃_{i=1}^{4} v_{i}d_{i}A By definition of \overline{i}

⋃_{i=1}^{4} αv\overline{i}d\overline{i}A = ⋃_{i=1}^{4} v_{i}d_{i}A

⋃_{i=1}^{4} v_{i}αd\overline{i}A = ⋃_{i=1}^{4} v_{i}d_{i}A

since the attractor is contained completely within each v_{i} we have

∀i, v_{i}αd\overline{i}A = v_{i}d_{i}A
∀i, αd\overline{i}A = d_{i}A
∀i, d_{i}^{−1}αd\overline{i}A = A

Thus ∀i, d_{i}^{−1}αd\overline{i} ∈ Sym (A) .

22 CARLA M. RIGGI

(<=) Assume ∀i, d_{i}^{−1}αd\overline{i} ∈ Sym (A) . Then ∀i, d_{i}^{−1}αd\overline{i} = β_{i} for some β_{i} ∈ Sym (A) and so ∀i, αd\overline{i} = d_{i}β_{i}.

α (A) = α (W (A))
= α (⋃_{i=1}^{4} v_{i}d_{i}A)
= α (⋃_{i=1}^{4} v\overline{i}d\overline{i}A)
= ⋃_{i=1}^{4} αv\overline{i}d\overline{i}A
= ⋃_{i=1}^{4} v_{i}αd\overline{i}A
= ⋃_{i=1}^{4} v_{i}d_{i}β_{i}A
= ⋃_{i=1}^{4} v_{i}d_{i}A
= W (A)
= A

Thus α ∈ Sym (A) . □ Theorem 3.16. Let α ∈ Δ and let W ∈ Ψ such that A = W (A) = ⋃_{i=1}^{4} v_{i}d_{i} (A).

If ∀i, αd\overline{i} = d_{i}α then α ∈ Sym (A)

Proof. Let α ∈ Δ and let W ∈ Ψ such that A = W (A) = ⋃_{i=1}^{4} v_{i}d_{i} (A). Assume ∀i, αd\overline{i} = d_{i}α then

α (A) = α(W (A))
= α ⋃_{i=1}^{4} v_{i}d_{i} (A)
= ⋃_{i=1}^{4} αv\overline{i}d\overline{i} (A)
= ⋃_{i=1}^{4} v_{i}αd\overline{i} (A)
= ⋃_{i=1}^{4} v_{i}d_{i}α (A)
= W (α (A))

Since the attractor of an IFS is unique, it follows that αA = A and so α ∈ Sym (A) . □

HUTCHINSON OPERATORS IN ℝ^{3} 23

Theorem 3.17. Let α ∈ Δ and let W, W′ ∈ Ψ such that W = ⋃_{i=1}^{4} v_{i}d_{i}, W′ (A) =

⋃_{i=1}^{4} v_{i}di^{′} (A) and A = W′ (A) .

A = W (A) ⇔ ∀i, d_{i}^{−1}d′_{i} ∈ Sym (A)

Proof. Let α ∈ Δ and let W, W′ ∈ Ψ such that W = ⋃_{i=1}^{4} v_{i}d_{i} and A = W′ (A) = ⋃_{i=1}^{4} v_{i}di^{′} (A)

(=>) Assume W (A) = A

W (A) = W′ (A)
⋃_{i=1}^{4} v_{i}d_{i}A = ⋃_{i=1}^{4} v_{i}d′_{i}A
∀i, v_{i}d_{i}A = v_{i}d′_{i}A
∀i, d_{i}A = d′_{i}A
∀i, A = d_{i}^{−1}d′_{i}A

therefore, ∀i, d_{i}^{−1}d′_{i} ∈ Sym (A) .

(<=) Assume ∀i, d_{i}^{−1}d′_{i} ∈ Sym (A). Then ∀i, A = d_{i}^{−1}d′_{i}A and so by substitution we have

W (A) = ⋃_{i=1}^{4} v_{i}d_{i} (A)
= ⋃_{i=1}^{4} v_{i}d_{i} (d_{i}^{−1}d′_{i}A)
= ⋃_{i=1}^{4} v_{i}d′_{i}A
= W′ (A)
= A

Thus W (A) = A. □

24 CARLA M. RIGGI

5. Appendix A

5.1. Symmetry Group Notation and CubeIFS Input Table. The group of 48 symmetries of the unit cube can be generated by the five following elements.

e = starting position, no rotation or reflection − = reflection through the plane parallel to the yz plane that cuts the cube in half

r  s

t

r = 90 deg rotation about the axis through the center of the cube from (\frac{1}{2}, 0, \frac{1}{2})

to (\frac{1}{2}, 1, \frac{1}{2}) CCW when looking down the y axis from − to +.

s = 90 deg rotation about the axis through the center of the cube from (0, \frac{1}{2}, \frac{1}{2})

to (1, \frac{1}{2}, \frac{1}{2}) CCW when looking down the x axis from − to +.

t = 90 deg rotation about the axis through the center of the cube from (\frac{1}{2}, \frac{1}{2}, 0)

to (\frac{1}{2}, \frac{1}{2}, 1) CCW when looking down the z axis from − to +.

Symmetry group elements are defined by composing the generator operations from left to right.

For example: −r^{2}s means first reflection followed by 180 degrees rotation in the r direction (literally

it means to perform r twice) followed by 90 degrees rotation in the s direction.

The following chart lists all 48 elements of the group. The labels along the top row and left column are equivalent notations that may be used in the CubeIFS routine that is explained in Appendix D.

F H  N  R  L  B
Up  e  s  s^{3}  t  t^{3}  r^{2}s^{2}
-Up −e −s −s^{3} −t −t^{3} −r^{2}s^{2}
Lt r rs rs^{3} rt rt^{3} r^{3}s^{2}
-Lt −r −rs −rs^{3} −rt −rt^{3} −r^{3}s^{2}
Dn r^{2} r^{2}s r^{2}s^{3} r^{2}t r^{2}t^{3}  s^{2}
-Dn −r^{2} −r^{2}s −r^{2}s^{3} −r^{2}t −r^{2}t^{3} −s^{2}
Rt r^{3} r^{3}s r^{3}s^{3} r^{3}t r^{3}t^{3}  rs^{2}
-Rt −r^{3} −r^{3}s −r^{3}s^{3} −r^{3}t −r^{3}t^{3} −rs^{2}

Definition 5.1. Let G be a group. The center of G is the subgroup

Z (G) = {z ∈ G : zd = dz for every d ∈ G} .

The center of our group of symmetry operations on the cube is {e, −s^{2}} . Definition 5.2. Let G be a group and α ∈ G. The centralizer of α is the subgroup

C (α) = {d ∈ G : dα = αd} .

As it turns out, the groups of centralizers of α ∈ Δ (and sometimes a coset of one of these groups) contain key elements for defining IFS’s W ∈ Ψ having symmetric attractors. Some particular definitions that we know of so far are given in Appendix E. Here we list the groups of centralizers for each α ∈ Δ (7).

C (α_{1}) = C (−t) = {e, r^{2}s^{2}, r^{2}t^{3}, r^{2}t, −r^{2}, −s^{2}, −t, −t^{3}}
C (α_{2}) = C (−r^{2}s^{3}) = {e, −e, r^{2}s^{3}, −r^{2}s^{3}, r^{2}s, −r^{2}s, s^{2}, −s^{2}}
C (α_{3}) = C (−r^{3}) = {e, r^{3}s^{2}, r^{2}, rs^{2}, −r, −r^{3}, −s^{2}, −r^{2}s^{2}}
C (α_{4}) = C (rt) = C (α_{5}) = C (r^{3}s^{3}) = {e, r^{3}s^{3}, rt, −rs, −s^{2}, −rt^{3}}

As you can see, for i ∈ {1, 2, 3} , |C (α_{i}) | = 8 and for j ∈ {4, 5} , |C (α_{j}) | = 6. It is also helpful to note that the only elements common among the centralizer groups are the elements of the center, {e, −s^{2}} .

Definition 5.3. Let H be a subgroup of a group G. The normalizer of H is the set N_{H} = {d ∈ G : d^{−1}Hd = H} .

It is well known that N_{H} is a subgroup of G that contains H.

Lemma 1. Let W = ⋃_{i=1}^{4} v_{i}d_{i} ∈ Ψ such that A = W (A) . If | Sym(A)| = 6, and hence

Sym (A) = Δ, then d_{1} ∈ N_{Δ}.

Here we list the normalizer of Δ and the normalizer of C (α_{i}) for all i ∈ {1, 2, 3, 4, 5} since they have been useful in determining the definitions of symmetric attractors in Ψ.

N_{Δ} = N_{C(α_{4})} = N_{C(α_{5})} = {e, rs^{2}, r^{2}s, r^{3}s^{3}, rt, r^{2}t^{3}, −r^{3}, −s^{2}, −t, −rs, −r^{2}s^{3}, −rt^{3}}
N_{C(α_{1})} = {e, r^{2}, s^{2}, t, t^{3}, r^{2}s^{2}, r^{2}t, r^{2}t^{3}, −e, −r^{2}, −s^{2}, −t, −t^{3}, −r^{2}s^{2}, −r^{2}t, −r^{2}t^{3}}
N_{C(α_{2})} = {e, r^{2}, s, s^{2}, s^{3}, r^{2}s, r^{2}s^{2}, r^{2}s^{3}, −e, −r^{2}, −s, −s^{2}, −s^{3}, −r^{2}s, −r^{2}s^{2}, −r^{2}s^{3}}
N_{C(α_{3})} = {e, r, r^{2}, r^{3}, s^{2}, rs^{2}, r^{2}s^{2}, r^{3}s^{2}, −e, −r, −r^{2}, −r^{3}, −s^{2}, −rs^{2}, −r^{2}s^{2}, −r^{3}s^{2}}

Notice here that |N_{Δ}| = 12 and for i ∈ {1, 2, 3} , |N_{C(α_{i})}| = 16. Here also, the only elements common to all normalizer groups are {e, −s^{2}} . However, the elements common to N_{C(α_{i})} for i ∈ {1, 2, 3} are {e, r^{2}, s^{2}, r^{2}s^{2}, −e, −r^{2}, −s^{2}, −r^{2}s^{2}} . Additionally notice that N_{Δ} and N_{C(α_{1})} share {e, −s^{2}, r^{2}t^{3}, −t} , N_{Δ} and N_{C(α_{2})} share {e, −s^{2}, r^{2}s, −r^{2}s^{3}} , and N_{Δ} and N_{C(α_{3})} share {e, −s^{2}, rs^{2}, −r^{3}} .

6. Appendix B

6.1. Affine Forms in ℝ^{3} for Use in the Chaos3d Maple Module. The following representations of affine maps are presented for use in the Maple Chaos3d module as structures for defining an affine map in ℝ^{3}. The Vector, Spherical, and Geometric forms are the three core methods for defining an affine map, T (x) = Mx + B. For the convenience of maple users I additionally define the Matrix form and Standard form data structures to allow greater flexibility in the syntax used to create an affine map using known values for M and B.

An affine map in ℝ^{3} is completely determined by where it sends four non-coplanar points. The simplest points to consider are (1, 0, 0), (0, 1, 0), (0, 0, 1), and the origin. Therefore, it is possible to define an affine map by its effect on the three base vectors, i⃗, j⃗, k⃗, and the translation of the origin. If we know that the origin is mapped to point t, and we know that point (1, 0, 0) is mapped to point p_{i}, then we can calculate p_{i} −t which we call v⃗_{i}, the image of i⃗ under the affine map before translation. A similar computation works for determining v⃗_{j} and v⃗_{k}whose values are used in the affineV3d data structure. For convenience in the Chaos3d Maple module, these values may also be given by spherical coordinates using the affineS3d data structure.

Definition 6.1. Let T : ℝ^{3} → ℝ^{3} given by T (x) = Mx + B. where M ∈ M_{3}(ℝ) defines the images of the basis vectors, i⃗, j⃗, and k⃗ under the affine map before translation, and B ∈ ℝ^{3} is the translation map on the origin. The following data structures define equivalent affine maps. Vector Form Let

v⃗_{i} = [ \begin{matrix}a_{i} \\ b_{i} \\ c_{i} \end{matrix} ], v⃗_{j} = [ \begin{matrix}a_{j} \\ b_{j} \\ c_{j} \end{matrix} ], v⃗_{k} = [ \begin{matrix}a_{k} \\ b_{k} \\ c_{k} \end{matrix} ], t⃗ = [ \begin{matrix}l \\ m \\ n \end{matrix} ]

affineV3d(v⃗_{i}, v⃗_{j}, v⃗_{k}, t⃗) defines

T ( \begin{matrix}x \\ y \\ z \end{matrix} ) = [ \begin{matrix}a_{i} & a_{j} & a_{k} \\ b_{i} & b_{j} & b_{k} \\ c_{i} & c_{j} & c_{k} \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}l \\ m \\ n \end{matrix} ]  (8)

Spherical Form Let r, s, and t ∈ ℝ be scaling coefficients for i⃗, j⃗, and k⃗ respectively. Let θ_{i}, θ_{j}, θ_{k}, ϕ_{i}, ϕ_{j},and ϕ_{k} be angle measurements in degrees. Then according to standard spherical coordinates, affineS3d (r, s, t, θ_{i}, ϕ_{i}, θ_{j}, ϕ_{j}, θ_{k}, ϕ_{k}, l, m, n) defines

T ( \begin{matrix}x \\ y \\ z \end{matrix} ) = [ \begin{matrix}cos θ_{i} sin ϕ_{i} & cos θ_{j} sin ϕ_{j} & cos θ_{k} sin ϕ_{k} \\ sin θ_{i} sin ϕ_{i} & sin θ_{j} sin ϕ_{j} & sin θ_{k} sin ϕ_{k} \\ cos ϕ_{i} & cos ϕ_{j} & cos ϕ_{k} \end{matrix} ] [ \begin{matrix}r & 0 & 0 \\ 0 & s & 0 \\ 0 & 0 & t \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}l \\ m \\ n \end{matrix} ]

= [ \begin{matrix}r cos θ_{i} sin ϕ_{i} & s cos θ_{j} sin ϕ_{j} & t cos θ_{k} sin ϕ_{k} \\ r sin θ_{i} sin ϕ_{i} & s sin θ_{j} sin ϕ_{j} & t sin θ_{k} sin ϕ_{k} \\ r cos ϕ_{i} & s cos ϕ_{j} & t cos ϕ_{k} \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}l \\ m \\ n \end{matrix} ]  (9)

Standard Form affine3d (a_{1}, a_{2}, a_{3}, b_{1}, b_{2}, b_{3}, c_{1}, c_{2}, c_{3}, t_{1}, t_{2}, t_{3}) defines

T ( \begin{matrix}x \\ y \\ z \end{matrix} ) = [ \begin{matrix}a_{1} & a_{2} & a_{3} \\ b_{1} & b_{2} & b_{3} \\ c_{1} & c_{2} & c_{3} \end{matrix} ] ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}t_{1} \\ t_{2} \\ t_{3} \end{matrix} ]  (10)

Matrix Form Let M be a Maple-defined “Matrix” or “matrix”. Let B be a Maple-defined “Vector” affineM3d (M, B) defines

T (x) = Mx + B  (11)

Any three-dimensional affine map will stretch and twist the unit cube into any desired parallelepiped fixing the origin and then translate that shape to any coordinate in the ℝ^{3} plane. However, for many applications it may be useful to rotate a cube-like figure uniformly around one axis without destroying its orthogonality. For this purpose, the affineG3d data structure produces an affine map that first scales i⃗, j⃗, and k⃗ by real-valued constants r, s,and t then consecutively rotates this orthogonal parallelepiped about the x, y, and z axes by angles θ, ϕ,and γ respectively before translating the origin to B. Definition 6.2. Let T : ℝ^{3} → ℝ^{3} given by T (x) = Mx + B. where M ∈ M_{3}(ℝ) defines the images of the basis vectors, i⃗, j⃗, and k⃗ under the affine map before translation, and B ∈ ℝ^{3} is the translation map on the origin. Let θ, ϕ, γ be angle measurements in degrees and let

M = [ \begin{matrix}cos γ & − sin γ & 0 \\ sin γ & cos γ & 0 \\ 0 & 0 & 1 \end{matrix} ] [ \begin{matrix}cos ϕ & 0 & sin ϕ \\ 0 & 1 & 0 \\ − sin ϕ & 0 & cos ϕ \end{matrix} ] [ \begin{matrix}1 & 0 & 0 \\ 0 & cos θ & − sin θ \\ 0 & sin θ & cos θ \end{matrix} ] [ \begin{matrix}r & 0 & 0 \\ 0 & s & 0 \\ 0 & 0 & t \end{matrix} ] (12)

= [ \begin{matrix}r cos ϕ cos γ & s sin θ sin ϕ cos γ − s cos θ sin γ & t cos θ sin ϕ cos γ + t sin θ sin γ \\ r cos ϕ sin γ & s sin θ sin ϕ sin γ + s cos θ cos γ & t cos θ sin ϕ sin γ − t sin θ cos γ \\ −r sin ϕ & s sin θ cos ϕ & t cos θ cos ϕ \end{matrix} ]

The data structure that defines an affine map for any orthogonal parallelepiped is given in Geometric form by affineG3d(r, s, t, θ, ϕ, γ, l, m, n) which defines

T ( \begin{matrix}x \\ y \\ z \end{matrix} ) = M ( \begin{matrix}x \\ y \\ z \end{matrix} ) + [ \begin{matrix}l \\ m \\ n \end{matrix} ]  (13)

7. Appendix C

7.1. Chaos3d Maple Module. The following Maple programming was written as a threedimensional extension to the Chaos module written by Ken Monks for working with IFS’s and their attractors in ℝ^{2}. The Chaos3d module uses the five new affine data structures defined in Appendix B in order to define IFS’s on ℝ^{3} and then render their attractors either randomly or deterministically. A procedure is included for calculating a scaling/contraction factor for any affine map. The CubeIFS procedure is very useful for studying the family, Ψ, of relatives to the Sierpinski right tetrahedron as defined by 3.10, and other cube-based IFS’s. The basic IFS definitions for the Menger Sponge and the standard Sierpinski Tetrahedron are stored within the Chaos3d module and can be used within any procedure accepting IFS input. Multiple seed figures are also stored for use with the DrawDetIFS3d command for rendering an attractor. One procedure is included for writing a POV-Ray (Persistence Of Vision Ray tracing) file to store data for small spheres corresponding to points in the attractor of an IFS. However, this procedure uses a random iteration method. To produce a more accurate fractal image, the procedure should be altered to use a deterministic method.

Detailed explanations of the procedure calls and how to use them are given in Appendix D.

# Define data types
#
# 3d Vector form       -- affineV3d(vi,vj,vk,T)
# 3d Spherical form    -- affineS3d(r,s,t,thetaY,thetaZ,phiX,phiZ,gammaX,gammaY,l,m,n)
# 3d Geometric form    -- affineS3d(r,s,t,theta,phi,gamma,l,m,n)
# 3d Matrix form       -- affineM3d(M,B)
# 3d Standard form     -- affine3d(a1,a2,a3,b1,b2,b3,c1,c2,c3,t1,t2,t3)
# IFS = List of affines (of any type)
chaos3d:=module()
export
   affineV3dFromPoints,ContractionFactor3d,IFSFromList,CubeIFS,MengerSponge,
   SierpinskiTetrahedron,MrBlockHead,MrBlockFrame,MrCorner,MrBlock,MrPyramid,MrPoint3d,
   DrawDetIFS3d,DrawIFS3d,WriteIFSPOV;
global
   `convert/affineV3d`,`convert/affineS3d`,`convert/affineM3d`,
   `convert/affine3d`,`convert/affine3df`:
local
   affineV3d2affineS3d,affineS3d2affineV3d,affineV3d2affineM3d,affineM3d2affineV3d,
   affineG3d2affineV3d,affine3d2affineM3d,affine3d2affineV3d,MapOne3d,
   Map3d,Transform3d,ApplyIFS3d;
with(plots):
with(plottools):
protect(affineS3d,affineG3d,affineV3d,affineM3d,affine3d,affine3df):
# internal affine form conversion
affineV3d2affineS3d:=proc(vi,vj,vk,T)
   local r,s,t,thetaI,phiI,thetaJ,phiJ,thetaK,phiK,l,m,n:

r:= sqrt(op(1,vi)^2+op(2,vi)^2+op(3,vi)^2);
phiI:= arccos(op(3,vi)/r);
thetaI:= arctan(op(2,vi),op(1,vi));
s:= sqrt(op(1,vj)^2+op(2,vj)^2+op(3,vj)^2);
phiJ:= arccos(op(3,vj)/s);
thetaJ:= arctan(op(2,vj),op(1,vj));
t:= sqrt(op(1,vk)^2+op(2,vk)^2+op(3,vk)^2);
phiK:= arccos(op(3,vk)/t);
thetaK:= arctan(op(2,vk),op(1,vk));
l:=op(1,T);
m:=op(2,T);
n:=op(3,T);
RETURN(r,s,t,evalf(thetaI*180/Pi),evalf(phiI*180/Pi),
         evalf(thetaJ*180/Pi),evalf(phiJ*180/Pi),
         evalf(thetaK*180/Pi),evalf(phiK*180/Pi),l,m,n);
end:
affineS3d2affineV3d:=proc(r,s,t,thetaI,phiI,thetaJ,phiJ,thetaK,phiK,l,m,n)
local vi,vj,vk,T:
vi:=[evalf(r*cos(thetaI*(Pi/180))*sin(phiI*(Pi/180))),
   evalf(r*sin(thetaI*(Pi/180))*sin(phiI*(Pi/180))),
   evalf(r*cos(phiI*(Pi/180)))];
vj:=[evalf(s*cos(thetaJ*(Pi/180))*sin(phiJ*(Pi/180))),
   evalf(s*sin(thetaJ*(Pi/180))*sin(phiJ*(Pi/180))),
   evalf(s*cos(phiJ*(Pi/180)))];
vk:=[evalf(t*cos(thetaK*(Pi/180))*sin(phiK*(Pi/180))),
   evalf(t*sin(thetaK*(Pi/180))*sin(phiK*(Pi/180))),
   evalf(t*cos(phiK*(Pi/180)))];
T:=[l,m,n];
RETURN(vi,vj,vk,T);
end:
affineG3d2affineV3d:=proc(r,s,t,theta,phi,gamma,l,m,n)
local vi,vj,vk,T:
vi:=[evalf(r*cos(phi*Pi/180)*cos(gamma*Pi/180)),
   evalf(r*cos(phi*Pi/180)*sin(gamma*Pi/180)),
   evalf(-r*sin(phi*Pi/180))];
vj:=[evalf(s*sin(theta*Pi/180)*sin(phi*Pi/180)*cos(gamma*Pi/180)
      -s*cos(theta*Pi/180)*sin(gamma*Pi/180)),
   evalf(s*sin(theta*Pi/180)*sin(phi*Pi/180)*sin(gamma*Pi/180)
      +s*cos(theta*Pi/180)*cos(gamma*Pi/180)),
   evalf(s*sin(theta*Pi/180)*cos(phi*Pi/180))];
vk:=[evalf(t*cos(theta*Pi/180)*sin(phi*Pi/180)*cos(gamma*Pi/180)
      +t*sin(theta*Pi/180)*sin(gamma*Pi/180)),
   evalf(t*cos(theta*Pi/180)*sin(phi*Pi/180)*sin(gamma*Pi/180)
      -t*sin(theta*Pi/180)*cos(gamma*Pi/180)),
   evalf(t*cos(theta*Pi/180)*cos(phi*Pi/180))];
T:=[l,m,n];

RETURN(vi,vj,vk,T);
end:
affineV3d2affineM3d:=proc(vi,vj,vk,T)
local M,B:
M:=Matrix([[op(1,vi),op(2,vi),op(3,vi)],[op(1,vj),
   op(2,vj),op(3,vj)],[op(1,vk),op(2,vk),op(3,vk)]],scan=columns);
B:=Vector([op(1,T),op(2,T),op(3,T)]);
RETURN(M,B);
end:
affineM3d2affineV3d:=proc(M,B)
local vi,vj,vk,T:
vi:=[M[1,1],M[2,1],M[3,1]];
vj:=[M[1,2],M[2,2],M[3,2]];
vk:=[M[1,3],M[2,3],M[3,3]];
T :=[B[1],B[2],B[3]];
RETURN(vi,vj,vk,T);
end:
affine3d2affineM3d:=proc(a1,a2,a3,b1,b2,b3,c1,c2,c3,t1,t2,t3)
local M,B:
M:=Matrix([[a1,a2,a3],[b1,b2,b3],[c1,c2,c3]]):
B:=Vector([t1,t2,t3]):
RETURN(M,B);
end:
affine3d2affineV3d:=proc(a1,a2,a3,b1,b2,b3,c1,c2,c3,t1,t2,t3)
local vi,vj,vk,T:
vi:=[a1,b1,c1]:
vj:=[a2,b2,c2]:
vk:=[a3,b3,c3]:
T:=[t1,t2,t3]:
RETURN(vi,vj,vk,T);
end:
# user affine form conversion
`convert/affineM3d`:=proc(A)
local nm:
if type(A,list) then RETURN(map(`convert/affineM3d`,A)) fi;
nm:=op(0,A):
if nm=affineS3d then
   RETURN(affineM3d(affineV3d2affineM3d(affineS3d2affineV3d(op(A)))))
elif nm=affineG3d then
   RETURN(affineM3d(affineV3d2affineM3d(affineG3d2affineV3d(op(A)))))
elif nm=affineV3d then
   RETURN(affineM3d(affineV3d2affineM3d(op(A))))

elif nm=affine3d then
   RETURN(affineM3d(affine3d2affineM3d(op(A))))
else RETURN(A)
fi;
end:
`convert/affineV3d`:=proc(A)
local nm;
if type(A,list) then RETURN(map(`convert/affineV3d`,A)) fi;
nm:=op(0,A):
if nm=affineS3d then
   RETURN(affineV3d(affineS3d2affineV3d(op(A))))
elif nm=affineG3d then
   RETURN(affineV3d(affineG3d2affineV3d(op(A))))
elif nm=affineM3d then
   RETURN(affineV3d(affineM3d2affineV3d(op(A))))
elif nm=affine3d then
   RETURN(affineV3d(affine3d2affineV3d(op(A))))
else RETURN(A)
fi;
end:
`convert/affineS3d`:=proc(A)
local nm:
if type(A,list) then RETURN(map(`convert/affineS3d`,A)) fi;
nm:=op(0,A):
if nm=affineV3d then
   RETURN(affineS3d(affineV3d2affineS3d(op(A))))
elif nm=affine3d then
   RETURN(affineS3d(affineV3d2affineS3d(affine3d2affineV3d(op(A)))))
elif nm=affineM3d then
   RETURN(affineS3d(affineV3d2affineS3d(affineM3d2affineV3d(op(A)))))
elif nm=affineG3d then
   RETURN(affineS3d(affineV3d2affineS3d(affineG3d2affineV3d(op(A)))))
else RETURN(A)
fi;
end:
`convert/affine3d`:=proc(A)
local nm,aff;
if type(A,list) then RETURN(map(`convert/affine3d`,A)) fi;
nm:=op(0,A):
if nm=affineM3d then
   RETURN(affine3d(op(1,A)[1,1],op(1,A)[1,2],op(1,A)[1,3],
      op(1,A)[2,1],op(1,A)[2,2],op(1,A)[2,3],op(1,A)[3,1],
      op(1,A)[3,2],op(1,A)[3,3],op(2,A)[1],op(2,A)[2],op(2,A)[3]))
elif nm=affineV3d then
   RETURN(affine3d(op(1,A)[1],op(2,A)[1],op(3,A)[1],op(1,A)[2],

      op(2,A)[2],op(3,A)[2],op(1,A)[3],op(2,A)[3],
      op(3,A)[3],op(4,A)[1],op(4,A)[2],op(4,A)[3]))
elif nm=affineS3d then
   aff:=affineV3d(affineS3d2affineV3d(op(A))):
   RETURN(affine3d(op(1,aff)[1],op(2,aff)[1],op(3,aff)[1],
      op(1,aff)[2],op(2,aff)[2],op(3,aff)[2],op(1,aff)[3],
      op(2,aff)[3],op(3,aff)[3],op(4,aff)[1],op(4,aff)[2],op(4,aff)[3]))
elif nm=affineG3d then
   aff:=affineV3d(affineG3d2affineV3d(op(A))):
   RETURN(affine3d(op(1,aff)[1],op(2,aff)[1],op(3,aff)[1],
      op(1,aff)[2],op(2,aff)[2],op(3,aff)[2],op(1,aff)[3],
      op(2,aff)[3],op(3,aff)[3],op(4,aff)[1],op(4,aff)[2],op(4,aff)[3]))
else RETURN(A)
fi;
end:
`convert/affine3df`:=proc(A)
if type(A,list) then RETURN(map(`convert/affine3df`,A)) fi;
map(evalf,convert(A,affine3d)):
end:
# alternate affine map definition
affineV3dFromPoints:=proc(p1,p2,p3,p4,Tp1,Tp2,Tp3,Tp4)
local x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,
   ai,bi,ci,aj,bj,cj,ak,bk,ck,l,m,n,
   e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12:
x1:=p1[1];   y1:=p1[2];   z1:=p1[3];
x2:=p2[1];   y2:=p2[2];   z2:=p2[3];
x3:=p3[1];   y3:=p3[2];   z3:=p3[3];
x4:=p4[1];    y4:=p4[2];    z4:=p4[3];
e1:=ai*x1+aj*y1+ak*z1+l=Tp1[1];
e2:=bi*x1+bj*y1+bk*z1+m=Tp1[2];
e3:=ci*x1+cj*y1+ck*z1+n=Tp1[3];
e4:=ai*x2+aj*y2+ak*z2+l=Tp2[1];
e5:=bi*x2+bj*y2+bk*z2+m=Tp2[2];
e6:=ci*x2+cj*y2+ck*z2+n=Tp2[3];
e7:=ai*x3+aj*y3+ak*z3+l=Tp3[1];
e8:=bi*x3+bj*y3+bk*z3+m=Tp3[2];
e9:=ci*x3+cj*y3+ck*z3+n=Tp3[3];
e10:=ai*x4+aj*y4+ak*z4+l=Tp4[1];
e11:=bi*x4+bj*y4+bk*z4+m=Tp4[2];
e12:=ci*x4+cj*y4+ck*z4+n=Tp4[3];
assign(solve({e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12}));
affineV3d([ai,bi,ci],[aj,bj,cj],[ak,bk,ck],[l,m,n]);
end:
# calculate scaling factor of an affine map

ContractionFactor3d:=proc(aff)
local MV,Mt,E,A:
if type(aff,list) then
   RETURN(map(`ContractionFactor3d`,aff))
fi:
MV:=op(convert(aff,affineM3d)):
Mt:=LinearAlgebra[Transpose](MV[1]):
A:=Matrix(evalm(Mt*MV[1])):
E:= linalg['eigenvals'](A):
RETURN(Re(evalf(sqrt(E[-1]))));
end:
# methods for defining an IFS
IFSFromList:=proc(List)    # assumed to be a list of affineV3d coefficients
map(x->affineV3d(op(x)),List);
end:
CubeIFS:=proc()
local check,tcheck,n,m,p,i,j,k,Lst,T,data,h,A,lab,d,names:
check:=evalf(nargs^1/3));   tcheck:=trunc(check);
if not check=tcheck then
   ERROR(''The number of arguments must be a perfect cube.'')
fi:
n:=tcheck:   m:=1:   Lst:=NULL;   p:=1/n;
A:={args};   data:=[args];
if not hastype(A,list) then
   A:=[args];   data:=[];
   lab:=[e,-e,r,-r,r2,-r2,r3,-r3,s,-s,rs,-rs,r2s,-r2s,r3s,-r3s,
      s2,-s2,rs2,-rs2,r2s2,-r2s2,r3s2,-r3s2,s3,-s3,rs3,-rs3,
      r2s3,-r2s3,r3s3,-r3s3,t,-t,rt,-rt,r2t,-r2t,r3t,-r3t,t3,-t3,rt3,
      -rt3,r2t3,-r2t3,r3t3,-r3t3];
   names:=[
      [e,[F,Up]],[-e,[F,-Up]],[r,[F,Lt]],[-r,[F,-Lt]],
      [r2,[F,Dn]],[-r2,[F,-Dn]],[r3,[F,Rt]],[-r3,[F,-Rt]],
      [s,[H,Up]],[-s,[H,-Up]],[rs,[H,Lt]],[-rs,[H,-Lt]],
      [r2s,[H,Dn]],[-r2s,[H,-Dn]],[r3s,[H,Rt]],[-r3s,[H,-Rt]],
      [s2,[B,Dn]],[-s2,[B,-Dn]],[rs2,[B,Rt]],[-rs2,[B,-Rt]],
      [r2s2,[B,Up]],[-r2s2,[B,-Up]],[r3s2,[B,Lt]],[-r3s2,[B,-Lt]],
      [s3,[N,Up]],[-s3,[N,-Up]],[rs3,[N,Lt]],[-rs3,[N,-Lt]],
      [r2s3,[N,Dn]],[-r2s3,[N,-Dn]],[r3s3,[N,Rt]],[-r3s3,[N,-Rt]],
      [t,[R,Up]],[-t,[R,-Up]],[rt,[R,Lt]],[-rt,[R,-Lt]],
      [r2t,[R,Dn]],[-r2t,[R,-Dn]],[r3t,[R,Rt]],[-r3t,[R,-Rt]],
      [t3,[L,Up]],[-t3,[L,-Up]],[rt3,[L,Lt]],[-rt3,[L,-Lt]],
      [r2t3,[L,Dn]],[-r2t3,[L,-Dn]],[r3t3,[L,Rt]],[-r3t3,[L,-Rt]]];
   for h to nargs do
      if member(A[h],lab) then

for d to 48 do
if A[h]=names[d][1] then
   data:=[op(data),names[d][2]];
   break;
fi;
od;
else
data:=[op(data),A[h]];
fi;
od;
fi;
for j from 0 to n-1 do
for k from 0 to n-1 do
for i from 0 to n-1 do
if op(1,data[m])=F then
if op(2,data[m])=Up then
   T:=affineM3d(Matrix(Vector([p,p,p]),
      shape=diagonal),Vector([i*p,j*p,k*p])):
elif op(2,data[m])=-Up then
   T:=affineM3d(Matrix(Vector([-p,p,p]),
      shape=diagonal),Vector([(i+1)*p,j*p,k*p])):
elif op(2,data[m])=Dn then
   T:=affineM3d(Matrix(Vector([-p,p,-p]),
      shape=diagonal),Vector([(i+1)*p,j*p,(k+1)*p])):
elif op(2,data[m])=-Dn then
   T:=affineM3d(Matrix(Vector([p,p,-p]),
      shape=diagonal),Vector([i*p,j*p,(k+1)*p])):
elif op(2,data[m])=Rt then
   T:=affineM3d(Matrix([[0,0,-p],[0,p,0],[p,0,0]],
      scan=columns),Vector([i*p,j*p,(k+1)*p])):
elif op(2,data[m])=-Rt then
   T:=affineM3d(Matrix([[0,0,p],[0,p,0],[p,0,0]],
      scan=columns),Vector([i*p,j*p,k*p])):
elif op(2,data[m])=Lt then
   T:=affineM3d(Matrix([[0,0,p],[0,p,0],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,j*p,k*p])):
elif op(2,data[m])=-Lt then
   T:=affineM3d(Matrix([[0,0,-p],[0,p,0],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,j*p,(k+1)*p])):
else T:=NULL;
fi:
elif op(1,data[m])=B then
if op(2,data[m])=Up then
   T:=affineM3d(Matrix(Vector([-p,-p,p]),
      shape=diagonal),Vector([(i+1)*p,(j+1)*p,k*p])):
elif op(2,data[m])=-Up then
   T:=affineM3d(Matrix(Vector([p,-p,p]),
      shape=diagonal),Vector([i*p,(j+1)*p,k*p])):

elif op(2,data[m])=Dn then
   T:=affineM3d(Matrix(Vector([p,-p,-p]),
      shape=diagonal),Vector([i*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=-Dn then
   T:=affineM3d(Matrix(Vector([-p,-p,-p]),
      shape=diagonal),Vector([(i+1)*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=Rt then
   T:=affineM3d(Matrix([[0,0,-p],[0,-p,0],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=-Rt then
   T:=affineM3d(Matrix([[0,0,p],[0,-p,0],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,k*p])):
elif op(2,data[m])=Lt then
   T:=affineM3d(Matrix([[0,0,p],[0,-p,0],[p,0,0]],
      scan=columns),Vector([i*p,(j+1)*p,k*p])):
elif op(2,data[m])=-Lt then
   T:=affineM3d(Matrix([[0,0,-p],[0,-p,0],[p,0,0]],
      scan=columns),Vector([i*p,(j+1)*p,(k+1)*p])):
else T:=NULL;
fi:
elif op(1,data[m])=L then
if op(2,data[m])=Up then
   T:=affineM3d(Matrix([[0,p,0],[-p,0,0],[0,0,p]],
      scan=columns),Vector([(i+1)*p,j*p,k*p])):
elif op(2,data[m])=-Up then
   T:=affineM3d(Matrix([[0,-p,0],[-p,0,0],[0,0,p]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,k*p])):
elif op(2,data[m])=Dn then
   T:=affineM3d(Matrix([[0,-p,0],[-p,0,0],[0,0,-p]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=-Dn then
   T:=affineM3d(Matrix([[0,p,0],[-p,0,0],[0,0,-p]],
      scan=columns),Vector([(i+1)*p,j*p,(k+1)*p])):
elif op(2,data[m])=Rt then
   T:=affineM3d(Matrix([[0,0,-p],[-p,0,0],[0,p,0]],
      scan=columns),Vector([(i+1)*p,j*p,(k+1)*p])):
elif op(2,data[m])=-Rt then
   T:=affineM3d(Matrix([[0,0,p],[-p,0,0],[0,p,0]],
      scan=columns),Vector([(i+1)*p,j*p,k*p])):
elif op(2,data[m])=Lt then
   T:=affineM3d(Matrix([[0,0,p],[-p,0,0],[0,-p,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,k*p])):
elif op(2,data[m])=-Lt then
   T:=affineM3d(Matrix([[0,0,-p],[-p,0,0],[0,-p,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,(k+1)*p])):
else T:=NULL;
fi:
elif op(1,data[m])=R then

if op(2,data[m])=Up then
   T:=affineM3d(Matrix([[0,-p,0],[p,0,0],[0,0,p]],
      scan=columns),Vector([i*p,(j+1)*p,k*p])):
elif op(2,data[m])=-Up then
   T:=affineM3d(Matrix([[0,p,0],[p,0,0],[0,0,p]],
      scan=columns),Vector([i*p,j*p,k*p])):
elif op(2,data[m])=Dn then
   T:=affineM3d(Matrix([[0,p,0],[p,0,0],[0,0,-p]],
      scan=columns),Vector([i*p,j*p,(k+1)*p])):
elif op(2,data[m])=-Dn then
   T:=affineM3d(Matrix([[0,-p,0],[p,0,0],[0,0,-p]],
      scan=columns),Vector([i*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=Rt then
   T:=affineM3d(Matrix([[0,0,-p],[p,0,0],[0,-p,0]],
      scan=columns),Vector([i*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=-Rt then
   T:=affineM3d(Matrix([[0,0,p],[p,0,0],[0,-p,0]],
      scan=columns),Vector([i*p,(j+1)*p,k*p])):
elif op(2,data[m])=Lt then
   T:=affineM3d(Matrix([[0,0,p],[p,0,0],[0,p,0]],
      scan=columns),Vector([i*p,j*p,k*p])):
elif op(2,data[m])=-Lt then
   T:=affineM3d(Matrix([[0,0,-p],[p,0,0],[0,p,0]],
      scan=columns),Vector([i*p,j*p,(k+1)*p])):
else T:=NULL;
fi:
elif op(1,data[m])=H then
if op(2,data[m])=Up then
   T:=affineM3d(Matrix([[p,0,0],[0,0,-p],[0,p,0]],
      scan=columns),Vector([i*p,j*p,(k+1)*p])):
elif op(2,data[m])=-Up then
   T:=affineM3d(Matrix([[-p,0,0],[0,0,-p],[0,p,0]],
      scan=columns),Vector([(i+1)*p,j*p,(k+1)*p])):
elif op(2,data[m])=Dn then
   T:=affineM3d(Matrix([[-p,0,0],[0,0,-p],[0,-p,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=-Dn then
   T:=affineM3d(Matrix([[p,0,0],[0,0,-p],[0,-p,0]],
      scan=columns),Vector([i*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=Rt then
   T:=affineM3d(Matrix([[0,-p,0],[0,0,-p],[p,0,0]],
      scan=columns),Vector([i*p,(j+1)*p,(k+1)*p])):
elif op(2,data[m])=-Rt then
   T:=affineM3d(Matrix([[0,p,0],[0,0,-p],[p,0,0]],
      scan=columns),Vector([i*p,j*p,(k+1)*p])):
elif op(2,data[m])=Lt then
   T:=affineM3d(Matrix([[0,p,0],[0,0,-p],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,j*p,(k+1)*p])):

elif op(2,data[m])=-Lt then
   T:=affineM3d(Matrix([[0,-p,0],[0,0,-p],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,(k+1)*p])):
else T:=NULL;
fi:
elif op(1,data[m])=N then
if op(2,data[m])=Up then
   T:=affineM3d(Matrix([[p,0,0],[0,0,p],[0,-p,0]],
      scan=columns),Vector([i*p,(j+1)*p,k*p])):
elif op(2,data[m])=-Up then
   T:=affineM3d(Matrix([[-p,0,0],[0,0,p],[0,-p,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,k*p])):
elif op(2,data[m])=Dn then
   T:=affineM3d(Matrix([[-p,0,0],[0,0,p],[0,p,0]],
      scan=columns),Vector([(i+1)*p,j*p,k*p])):
elif op(2,data[m])=-Dn then
   T:=affineM3d(Matrix([[p,0,0],[0,0,p],[0,p,0]],
      scan=columns),Vector([i*p,j*p,k*p])):
elif op(2,data[m])=Rt then
   T:=affineM3d(Matrix([[0,p,0],[0,0,p],[p,0,0]],
      scan=columns),Vector([i*p,j*p,k*p])):
elif op(2,data[m])=-Rt then
   T:=affineM3d(Matrix([[0,-p,0],[0,0,p],[p,0,0]],
      scan=columns),Vector([i*p,(j+1)*p,k*p])):
elif op(2,data[m])=Lt then
   T:=affineM3d(Matrix([[0,-p,0],[0,0,p],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,(j+1)*p,k*p])):
elif op(2,data[m])=-Lt then
   T:=affineM3d(Matrix([[0,p,0],[0,0,p],[-p,0,0]],
      scan=columns),Vector([(i+1)*p,j*p,k*p])):
else T:=NULL;
fi:
else T:=NULL;
fi:
Lst:=Lst,T:
m:=m+1;
od:
od:
od:
[  Lst];
end:
# Procedures for rendering IFS attractors
# Random method
DrawIFS3d:=proc(IFS,n)
local i,m,a,b,c,t,ifs,dice,r,pts,FastIFSData,s,x,y,z;
FastIFSData:=proc(a,b,c,t,r,n,s)

local i,j,k,pts;
pts:=hfarray(1..n,1..3);
pts[1,1]:=s[0]: pts[1,2]:=s[1]: pts[1,3]:=s[2]:
for i from 2 to n do
   k:=r[i];   j:=i-1;
   pts[i,1]:=a[0,k]*pts[j,1]+a[1,k]*pts[j,2]+a[2,k]*pts[j,3]+t[0,k];
   pts[i,2]:=b[0,k]*pts[j,1]+b[1,k]*pts[j,2]+b[2,k]*pts[j,3]+t[1,k];
   pts[i,3]:=c[0,k]*pts[j,1]+c[1,k]*pts[j,2]+c[2,k]*pts[j,3]+t[2,k];
od:
pts;
end:
ifs:=convert(IFS,affine3df);
pts:=hfarray(1..n,1..3);
m:=nops(ifs)-1;
dice:=rand(0..m);
r:=hfarray(1..n,[seq(dice(),i=1..n)]);
a:=hfarray(0..2,0..m,[[seq(op(1,ifs[i+1]),i=0..m)],
[  seq(op(2,ifs[i+1]),i=0..m)],[seq(op(3,ifs[i+1]),i=0..m)]]);
b:=hfarray(0..2,0..m,[[seq(op(4,ifs[i+1]),i=0..m)],
[  seq(op(5,ifs[i+1]),i=0..m)],[seq(op(6,ifs[i+1]),i=0..m)]]);
c:=hfarray(0..2,0..m,[[seq(op(7,ifs[i+1]),i=0..m)],
[  seq(op(8,ifs[i+1]),i=0..m)],[seq(op(9,ifs[i+1]),i=0..m)]]);
t:=hfarray(0..2,0..m,[[seq(op(10,ifs[i+1]),i=0..m)],
[  seq(op(11,ifs[i+1]),i=0..m)],[seq(op(12,ifs[i+1]),i=0..m)]]);
# compute a fixed point to use as a seed
s:=hfarray(0..2);
assign(solve({x=a[0,0]*x+a[1,0]*y+a[2,0]*z+t[0,0],
y=b[0,0]*x+b[1,0]*y+b[2,0]*z+t[1,0],
z=c[0,0]*x+c[1,0]*y+c[2,0]*z+t[2,0]})):
s[0]:=x:   s[1]:=y:   s[2]:=z:
plots[display3d](PLOT3D(POINTS( evalhf(FastIFSData(a,b,c,t,r,n,s))) ),
scaling=constrained,axes=none,symbol=POINT,args[3..nargs]);
end:
# DrawDetIFS3d internal use only
MapOne3d:=proc(aff)
local x,y,z,a1,a2,a3,b1,b2,b3,c1,c2,c3,t1,t2,t3,A;
A:=convert(aff,affine3d);
a1:=op(1,A);   b1:=op(4,A);   c1:=op(7,A);
a2:=op(2,A);   b2:=op(5,A);    c2:=op(8,A);
a3:=op(3,A);   b3:=op(6,A);    c3:=op(9,A);
t1:=op(10,A);   t2:=op(11,A);   t3:=op(12,A);
unapply([a1*x+a2*y+a3*z+t1,b1*x+b2*y+b3*z+t2,c1*x+c2*y+c3*z+t3],x,y,z);
end:
Map3d:=proc(aff)
if type(aff,list) then map(MapOne3d,aff);
else MapOne3d(aff);

   fi:
end:
Transform3d:=proc(aff)
   if type(aff,list) then map(x->transform(Map3d(x)),aff);
   else transform(Map3d(aff));
   fi:
end:
ApplyIFS3d:=proc(image,IFSTrans)
   local i;
   plots[display3d]({seq(IFSTrans[i](image),i=1..nops(IFSTrans))});
end:
# Deterministic method
DrawDetIFS3d:=proc(image,IFS,n)
local P,T,i;
P:=image;
T:=Transform3d(IFS);
for i to n do
   P:=ApplyIFS3d(P,T);
od:
plots[display3d](P,args[4..nargs],scaling=CONSTRAINED);
end:
# POV-Ray Interface
# Random method of attractor rendering analogous to DrawIFS3d
WriteIFSPOV:=proc(NewFile,IFS,n)
local fileId,ifs,m,dice,r,a,b,c,t,pts,di,i,j,k:
fileId := fopen(NewFile,WRITE,TEXT):
ifs:=convert(IFS,affine3df);
m:=nops(ifs)-1;
dice:=rand(0..m);
r:=array(1..n,[seq(dice(),i=1..n)]);
a:=hfarray(0..2,0..m,[[seq(op(1,ifs[i+1]),i=0..m)],
      [  seq(op(2,ifs[i+1]),i=0..m)],[seq(op(3,ifs[i+1]),i=0..m)]]);
b:=hfarray(0..2,0..m,[[seq(op(4,ifs[i+1]),i=0..m)],
      [  seq(op(5,ifs[i+1]),i=0..m)],[seq(op(6,ifs[i+1]),i=0..m)]]);
c:=hfarray(0..2,0..m,[[seq(op(7,ifs[i+1]),i=0..m)],
      [  seq(op(8,ifs[i+1]),i=0..m)],[seq(op(9,ifs[i+1]),i=0..m)]]);
t:=hfarray(0..2,0..m,[[seq(op(10,ifs[i+1]),i=0..m)],
      [  seq(op(11,ifs[i+1]),i=0..m)],[seq(op(12,ifs[i+1]),i=0..m)]]);
pts:=hfarray(1..n,1..3);
# compute a fixed point to use as a seed
assign(solve({x=a[0,0]*x+a[1,0]*y+a[2,0]*z+t[0,0],
         y=b[0,0]*x+b[1,0]*y+b[2,0]*z+t[1,0],
         z=c[0,0]*x+c[1,0]*y+c[2,0]*z+t[2,0]})):
pts[1,1]:=x:    pts[1,2]:=y:   pts[1,3]:=z:

   fprintf(fileId,''sphere { <%13.10f,%13.10f,%13.10f> IFSRadius texture{
                  SphereTexIFS } }n'',pts[1,1],pts[1,2],pts[1,3]);
   for k from 2 to n do
      di:=r[k];   j:=k-1;
      pts[k,1]:=a[0,di]*pts[j,1]+a[1,di]*pts[j,2]+a[2,di]*pts[j,3]+t[0,di];
      pts[k,2]:=b[0,di]*pts[j,1]+b[1,di]*pts[j,2]+b[2,di]*pts[j,3]+t[1,di];
      pts[k,3]:=c[0,di]*pts[j,1]+c[1,di]*pts[j,2]+c[2,di]*pts[j,3]+t[2,di];
      fprintf(fileId,''sphere { <%13.10f,%13.10f,%13.10f> IFSRadius texture{
                     SphereTexIFS } }n'',pts[k,1],pts[k,2],pts[k,3]);
   od:
   fclose(fileId);
end:
# famous 3D IFSs
MengerSponge:=CubeIFS([F,Up],[F,Up],[F,Up],[F,Up],
               none,[F,Up],[F,Up],[F,Up],[F,Up],[F,Up],
               none,[F,Up],none,none,none,[F,Up],
               none,[F,Up],[F,Up],[F,Up],[F,Up],[F,Up],
               none,[F,Up],[F,Up],[F,Up],[F,Up]):
SierpinskiTetrahedron:=[affineG3d(.5,.5,.5,0,0,0,0,0,0),
                  affineG3d(.5,.5,.5,0,0,0,0,-.5,0),
                  affineG3d(.5,.5,.5,0,0,0,-0.5,-.25,0),
                  affineG3d(.5,.5,.5,0,0,0,-.25,-.25,.5)];
# Common starting figures for Deterministic IFS -- DrawDetIFS3d
MrPoint3d:=pointplot3d([0,0,0]):
MrCorner:=display3d(
         {line([0,0,0],[1,0,0]),
         line([0,0,0],[0,1,0]),
         line([0,0,0],[0,0,1])},
         scaling=CONSTRAINED):
MrPyramid:=display3d(
         plots[polygonplot3d]([
            [[.5,-evalf(sqrt(3)/6),0],[0,0,evalf(sqrt(2/3))],
             [0,evalf(sqrt(3)/3),0]],
            [[.5,-evalf(sqrt(3)/6),0],[0,0,evalf(sqrt(2/3))],
             [-.5,-evalf(sqrt(3)/6),0]],
            [[-.5,-evalf(sqrt(3)/6),0],[0,0,evalf(sqrt(2/3))],
             [0,evalf(sqrt(3)/3),0]],
            [[.5,-evalf(sqrt(3)/6),0],[0,evalf(sqrt(3)/3),0],
             [-.5,-evalf(sqrt(3)/6),0]]]),
         scaling=CONSTRAINED,
         style=PATCHNOGRID):

MrBlock:=display3d(
cuboid([0,0,0],[1,1,1]),
scaling=CONSTRAINED,
style=PATCHNOGRID):
MrBlockHead:=display3d(
{cuboid([0,0,0],[1,1,1]),
spacecurve({[.3+0.075*cos(t),0,.7+0.075*sin(t),t=0..2*Pi],
[t,0,.7,t=0.6..0.8],[.5+.25*cos(t),0,.5+.25*sin(t),t=-3*Pi/4..-Pi/4]}
,color=black),
spacecurve({[.5+.4*cos(t),.5+.4*sin(t),1,t=0..2*Pi],
[.5+.3*cos(t),.5+.3*sin(t),1,t=0..2*Pi]}
,color=black),
spacecurve({[1,.25+.25*cos(t),.5+.25*sin(t),t=-Pi/2..Pi/2],
[1,(.25+.25*cos(Pi/3))+.07*cos(t),(.5+.25*sin(-Pi/3))+.07*sin(t),
   t=-5*Pi/4..0]}
,color=black),
spacecurve({[0,.25+.25*sin(t),.5+.25*cos(t),t=0..Pi]}
,color=black),
spacecurve({[t,.5*t-.15,1,t=.4..0.6],[t,-.5*t+.45,1,t=.4..0.6],
[.4,s,1,s=.05..0.25]}
,color=black),
spacecurve({[(t+.54)/1.6,1,t,t=.41..0.9],[(t+.3)/1.6,1,t,t=.1..0.3],
[(t-1.06)/(-1.6),1,t,t=.5..0.9],[(t-1.3)/(-1.6),1,t,t=.1..0.3],
[(t+1.7)/4,1,t,t=.44..0.9],[(t+1.5)/4,1,t,t=.1..0.3],
[(t-2.3)/(-4),1,t,t=.47..0.9],[(t-2.5)/(-4),1,t,t=.1..0.3],
[t,1,-.4*t+.6,t=.3..0.7],
[s,1,.32,s=.3..0.7],
[.3,1,s,s=.32..0.48]}
,color=black),
spacecurve({[t,.5*t-.05,0,t=.3..0.7],[t,(-.5)*t+.45,0,t=.3..0.7],
[t,5*t-2.3,0,t=.5..0.6],[t,(-5)*t+2.7,0,t=.4..0.5],
[t,t+.3,0,t=.4..0.6],[t,1.5-t,0,t=.6..0.8],
[.3,t,0,t=.1..0.3],[.7,t,0,t=.1..0.3],[t,.7,0,t=.4..0.8]}
,color=black)
},
scaling=CONSTRAINED,
style=PATCHNOGRID):
MrBlockFrame:=display3d(
{line([0,0,0],[0,0,1]), line([0,0,1],[1,0,1]), line([1,0,0],[1,0,1]),
line([0,0,0],[1,0,0]), line([0,0,0],[0,1,0]), line([0,1,0],[0,1,1]),
line([0,1,1],[1,1,1]), line([0,1,0],[1,1,0]), line([1,0,0],[1,1,0]),
line([1,1,0],[1,1,1]), line([1,0,1],[1,1,1]), line([0,0,1],[0,1,1]),
spacecurve({[.3+0.075*cos(t),0,.7+0.075*sin(t),t=0..2*Pi],
[t,0,.7,t=0.6..0.8],[.5+.25*cos(t),0,.5+.25*sin(t),t=-3*Pi/4..-Pi/4]}
,color=red),

spacecurve({[.5+.4*cos(t),.5+.4*sin(t),1,t=0..2*Pi],
[.5+.3*cos(t),.5+.3*sin(t),1,t=0..2*Pi]}
,color=red),
spacecurve({[1,.25+.25*cos(t),.5+.25*sin(t),t=-Pi/2..Pi/2],
[1,(.25+.25*cos(Pi/3))+.07*cos(t),(.5+.25*sin(-Pi/3))+.07*sin(t),
   t=-5*Pi/4..0]}
,color=red),
spacecurve({[0,.25+.25*sin(t),.5+.25*cos(t),t=0..Pi]}
,color=red),
spacecurve({[t,.5*t-.15,1,t=.4..0.6],[t,-.5*t+.45,1,t=.4..0.6],
[.4,s,1,s=.05..0.25]}
,color=red),
spacecurve({[(t+.54)/1.6,1,t,t=.41..0.9],[(t+.3)/1.6,1,t,t=.1..0.3],
[(t-1.06)/(-1.6),1,t,t=.5..0.9],[(t-1.3)/(-1.6),1,t,t=.1..0.3],
[(t+1.7)/4,1,t,t=.44..0.9],[(t+1.5)/4,1,t,t=.1..0.3],
[(t-2.3)/(-4),1,t,t=.47..0.9],[(t-2.5)/(-4),1,t,t=.1..0.3],
[t,1,-.4*t+.6,t=.3..0.7],
[s,1,.32,s=.3..0.7],
[.3,1,s,s=.32..0.48]}
,color=red),
spacecurve({[t,.5*t-.05,0,t=.3..0.7],[t,(-.5)*t+.45,0,t=.3..0.7],
[t,5*t-2.3,0,t=.5..0.6],[t,(-5)*t+2.7,0,t=.4..0.5],
[t,t+.3,0,t=.4..0.6],[t,1.5-t,0,t=.6..0.8],
[.3,t,0,t=.1..0.3],[.7,t,0,t=.1..0.3],[t,.7,0,t=.4..0.8]}
,color=red)
},
scaling=CONSTRAINED,
style=PATCHNOGRID):
end module:

8. Appendix D

8.1. Introduction to the Chaos3d module. Calling Sequence:

readlib(chaos3d) function ( args )

Description:

– To use a chaos3d function, you must first execute the readlib(chaos3d) command.

The functions available are: IFS’s and AFFINE maps In the following routines, we define two data types: AFFINE and IFS.

– An IFS is a single AFFINE or a list of one or more AFFINEs.

– An AFFINE represents an affine map of 3-dimensional space and comes in five forms: affine3d(a_{1}, a_{2}, a_{3}, b_{1}, b_{2}, b_{3}, c_{1}, c_{2}, c_{3}, t_{1}, t_{2}, t_{3}) - standard 3d form affineS3d (r, s, t, θ_{i}, ϕ_{i}, θ_{j}, ϕ_{j}, θ_{k}, ϕ_{k}, l, m, n) - spherical 3d form affineG3d(r, s, t, θ, ϕ, γ, l, m, n) - geometric 3d form affineV3d(v⃗_{i}, v⃗_{j}, v⃗_{k}, t⃗) - vector 3d form (where each of v⃗_{i}, v⃗_{j}, v⃗_{k}, t⃗ are Vectors or three-

element lists)

affineM3d(M, B) - matrix 3d form These are inert forms representing the maps defined in Appendix B.

Built-in IFS’s: The following IFS data types are built-in:

MengerSponge SierpinskiTetrahedron

Built-in starting figures: The following plot structures are built-in, and can be used as the starting figure for

plotting the iterations of the deterministic method for an IFS (see DrawDetIFS3d() below).

MrPoint3d MrBlockHead MrBlockFrame MrBlock MrCorner MrPyramid

IFS routines:

‘convert/affineV3d‘(A::AFFINE) - converts an AFFINE or IFS object to its affineV3d() form

‘convert/affineS3d‘(A::AFFINE) - converts an AFFINE or IFS object to its affineS3d() form

‘convert/affineM3d‘(A::AFFINE) - converts an AFFINE or IFS object to its affineM3d() form

‘convert/affine3d‘ (A::AFFINE) - converts an AFFINE or IFS object to its affine3d() form

‘convert/affine3df‘(A::AFFINE) - converts an AFFINE or IFS object to its affine3d() form

having all floating point arguments

affineV3dFromPoints(p_{1}, p_{2}, p_{3}, p_{4}, T_{p1}, T_{p2}, T_{p3}, T_{p4})

– Computes the AFFINE object affineV3d(v⃗_{i}, v⃗_{j}, v⃗_{k}, t⃗) for the unique affine map T (u⃗) = M(v⃗_{i}, v⃗_{j}, v⃗_{k})(u⃗) + t⃗ that maps points p_{1} to T_{p1}, p_{2} to T_{p2}, p_{3} to T_{p3}, and p_{4} to T_{p4}. These arguments are each lists containing three elements. A valid procedure call contains 8 lists of the form [x, y, z].

IFSFromList3d(List)

-creates an IFS data structure from a list of m lists of the form [v⃗_{i}, v⃗_{j}, v⃗_{k}, t⃗]where each of v⃗_{i}, v⃗_{j}, v⃗_{k}, t⃗ is a Vector or a list of the form [a, b, c] as would be found in an affineV3d call.

Map3d(IFS)

- converts the AFFINE or IFS object to the map T (x, y, z) = (a_{1} ∗ x + a_{2} ∗ y + a_{3} ∗ z + t_{1}, b_{1} ∗ x + b_{2} ∗ y + b_{3} ∗ z + t_{2}, c_{1} ∗ x + c_{2} ∗ y + c_{3} ∗ z + t_{3}).

Transform3d(IFS)

- converts the AFFINE or IFS object to a transform which can be applied to plot objects. DrawDetIFS3d(figure,IFS,n)

- plots the n^{th} iteration of the attractor of IFS using the deterministic method starting with figure as the seed. figure can be a predefined starting figure (see above) or another well-defined plot structure.

DrawIFS3d(IFS,n)

- plots the attractor of IFS by the random iteration method using n points. (Note: An image rendered by maple using this command will be extremely fuzzy unless MANY more points are rendered than you have time to wait for. Use of a value for n larger than about 3000 is not recommended. In general, DrawDetIFS3d() using MrPoint3d as figure will render a more accurate image more quickly.)

ContractionFactor3d(Aff)

- where Aff is an AFFINE object or a list of such objects. Returns the scaling/contraction factor(s) for the affine transformation(s) defined by the AFFINE object.

WriteIFSPOV(NewFile,IFS,n)

- creates a POV-Ray include file at the location specified by NewFile containing a list of n POV-Ray spheres each centered at a point in the attractor of IFS. NewFile is defined in quotes in this form: “c:\\documents\\NewFileName.inc”. The .inc file includes a list of these spheres with placeholders for user-defined radius, texture, etc.

In addition, classic IFS’s can be quickly produced with the following routine:

CubeIFS(seq)

- This routine takes n^{3} arguments. Each argument defines one affine map in an IFS based on the following structure: Begin with the unit cube [0..1] × [0..1] × [0..1] situated in x, y, z space such that the front lower left corner is located at the origin, the rear lower left corner is at (0, 1, 0), and the front upper left corner is located at (0, 0, 1). Thereby, the cube is resting with its base in quadrant I of the x, y plane and its body entirely in the +z half plane. This is standard position. Divide the unit cube into n by n smaller cubes all contained within the area of the original unit cube. Number the smaller cubes from 1 to n^{3} starting in the front lower left corner and moving from x = 0 to x = 1 (left to right) and then from z = 0 to z = 1 (bottom to top) and then from y = 0 to y = 1 (front to back). For example: If the unit cube is divided into 2 × 2 sub-cubes, view the structure by looking in the positive y-direction with the xz plane in front of you and the +x axis pointing to your right. Then cube #1 has its front lower left corner at the origin, cube #4 has its front upper right corner at (1, 0, 1), and cube #6 has its rear lower right corner at (1, 1, 0). The i^{th} argument of the CubeIFS procedure call describes an affine transformation which maps the unit cube into the ith smaller cube in a manner corresponding one of the 48 symmetry operations on the unit cube as defined in Appendix A.

- There are 2 acceptable forms for entering arguments (affine maps) into this IFS structure; symmetry notation form and MrBlockHead form. Differing notational forms may not be mixed in one CubeIFS call.

Symmetry Notation Form

Acceptable compositions of the generating elements for the group of symmetry operations on the cube are listed in the table in Appendix A. The terms along the top row and left column correspond to MrBlockHead form and may not be combined with symmetry notation in the same CubeIFS call.

When the i^{th} argument of a CubeIFS call contains none (or any unrecognized symbol), no affine map is assigned to that cube cell. If the i^{th} argument contains a symmetry operation from the table above, then the corresponding affine map performs the action of first positioning the cube according to the symmetry operation stated, then contracting the cube from its usual 1 × 1 × 1 size and finally relocating it to the i^{th} cell position.

MrBlockHead Form

Each argument in MrBlockHead form describes the appearance of the pre-defined starting figure, MrBlockHead, in one of his 48 possible orientations corresponding to the symmetries of the cube. An argument comes from the set

{[P, Q], none : P ∈ {F, B, R, L, H, N} and Q ∈ {Lt, Rt, Up, Dn, −Lt, −Rt, −Up, −Dn}}.

A precise definition of these arguments is given by the symmetry group table in Appendix A. The P value is along the top row and the Q value descends the left column. The ordered pair [P, Q] represents the affine map corresponding to one of the 48 elements in the symmetry group as defined by the said composition of the generating elements. When the i^{th} argument contains none, no affine map is assigned to that cube cell.

The visual interpretation of [P, Q] notation is as follows. If the i^{th} argument contains an ordered pair [P, Q], then the corresponding affine map performs the action of shrinking MrBlockHead from his usual 1 × 1 × 1 size and upright position into the i^{th} cell such that when viewed from standard position his face is now on the side of his head normally occupied by P and tipped in the direction defined by Q. The abbreviations in the set P stand for Face, Back, his Right ear, his Left ear, Hat, Necktie. The abbreviations in the set Q correspond to the direction his face is tipped (assuming you are now viewing him head-on) and mean: Up, eyes near top with open eye on your left; Dn, eyes near bottom with open eye on your right; Rt, eyes near your right with open eye on top; Lt, eyes near your left with open eye on bottom; −Up, eyes near top with open eye on your right; −Dn, eyes near bottom with open eye on your left; −Rt, eyes near your right with open eye on bottom; −Lt, eyes near your left with open eye on top.

Example: The pair [R, −Lt] is equivalent to the symmetry group element rt and it means that MrBlockHead’s face is now on the side of his head normally occupied by HIS Right ear, and if we now look at that side of his head, we notice that his eyes are toward OUR left with the open eye on top. (Note: when using this notation it is important to understand the correct position for [H, Up] and [N, Up] especially.)

[F, Up]  [R, −Lt]  [H, Up]  [N, Up]

Symmetry Notation form and MrBlockHead form produce identical IFS’s.

For example, one IFS that produces the Menger Sponge can be entered:

CubeIFS([F,Up],[F,Up],[F,Up],[F,Up],none,[F,Up],[F,Up],[F,Up],[F,Up],
[F,Up],none,[F,Up],none,none,none,[F,Up],none,[F,Up],
[F,Up],[F,Up],[F,Up],[F,Up],none,[F,Up],[F,Up],[F,Up],[F,Up]);

or using the symmetry notation:

CubeIFS(e,e,e,e,none,e,e,e,e,
e,none,e,none,none,none,e,none,e,
e,e,e,e,none,e,e,e,e);

However, it is INCORRECT to enter:

CubeIFS([F,Up],[F,Up],[F,Up],[F,Up],none,[F,Up],[F,Up],[F,Up],[F,Up],
 e,none,e,none,none,none,e,none,e,
 e,e,[F,Up],e,none,e,e,[F,Up],e);

9. Appendix E

9.1. Known Attractors of W ∈ Ψ having |S (A) | = 6. Let W ∈ Ψ be given by W (A) = v_{1}d_{1} (A) ∪ v_{2}d_{2} (A) ∪ v_{3}d_{3} (A) ∪ v_{4}d_{4} (A) . Then | Sym (A) | = 6 when d_{i} is chosen according to one of the eight possibilities below. Each attractor can be rendered by 6^{4} distinct IFS’s as determined by Theorem 3.14. The notation below lists all possible values for d_{i} in order to produce the corresponding attractor. These values are elements of the group of symmetry operations on the cube as defined in Appendix A.

Sierpinski Right Tetrahedron

∀i ∈ 1..4, d_{i} ∈ {e, rt, r^{3}s^{3}, −r^{3}, −r^{2}s^{3}, −t}

The Jax

∀i ∈ 1..4, d_{i} ∈ {−s^{2}, −rt^{3}, −rs, rs^{2}, r^{2}t^{3}, r^{2}s}

Triangular Plane (of area \sqrt{3})

d_{1} ∈ {−s^{2}, −rt^{3}, −rs, rs^{2}, r^{2}t^{3}, r^{2}s}
∀i ∈ 2..4, d_{i} ∈ {e, rt, r^{3}s^{3}, −r^{3}, −r^{2}s^{3}, −t}

The Bloom

d_{1} ∈ {e, rt, r^{3}s^{3}, −r^{3}, −r^{2}s^{3}, −t}
∀i ∈ 2..4, d_{i} ∈ {−s^{2}, −rt^{3}, −rs, rs^{2}, r^{2}t^{3}, r^{2}s}

Deadly Butterfly

d_{1} ∈ {e, rt, r^{3}s^{3}, −r^{3}, −r^{2}s^{3}, −t}
d_{2} ∈ {r^{3}t, −r^{2}t, s^{2}, −r^{2}s, r^{3}s, −r^{3}s^{2}}
d_{3} ∈ α_{4}d_{2} = {−t^{3}, rt^{3}, rs^{3}, r^{2}s^{2}, −rs^{2}, −s^{3}}
d_{4} ∈ α_{4}d_{3} = {r^{2}, −r, rs, −r^{2}t^{3}, r^{3}t^{3}, −s}

Galactica Rocket Ship

d_{1} ∈ {−s^{2}, −rt^{3}, −rs, rs^{2}, r^{2}t^{3}, r^{2}s}
d_{2} ∈ {r^{3}t, −r^{2}t, s^{2}, −r^{2}s, r^{3}s, −r^{3}s^{2}}
d_{3} ∈ α_{4}d_{2} = {−t^{3}, rt^{3}, rs^{3}, r^{2}s^{2}, −rs^{2}, −s^{3}}
d_{4} ∈ α_{4}d_{3} = {r^{2}, −r, rs, −r^{2}t^{3}, r^{3}t^{3}, −s}

“Mr. T”

d_{1} ∈ {e, rt, r^{3}s^{3}, −r^{3}, −r^{2}s^{3}, −t}
d_{2} ∈ {r, r^{2}s^{3}, t^{3}, −e, −rs^{3}, −r^{3}t^{3}}
d_{3} ∈ α_{4}d_{2} = {r^{2}t, s, −r^{3}s, −r^{2}, r^{3}, −rt}
d_{4} ∈ α_{4}d_{3} = {s^{3}, r^{3}s^{2}, −r^{2}s^{2}, −r^{3}t, t, −r^{3}s^{3}}

Klingon Battle Cruiser

d_{1} ∈ {−s^{2}, −rt^{3}, −rs, rs^{2}, r^{2}t^{3}, r^{2}s}
d_{2} ∈ {r, r^{2}s^{3}, t^{3}, −e, −rs^{3}, −r^{3}t^{3}}
d_{3} ∈ α_{4}d_{2} = {r^{2}t, s, −r^{3}s, −r^{2}, r^{3}, −rt}
d_{4} ∈ α_{4}d_{3} = {s^{3}, r^{3}s^{2}, −r^{2}s^{2}, −r^{3}t, t, −r^{3}s^{3}}

10. Appendix F

10.1. Extension to Maple’s group package; for using notation for the elements in the symmetry group of the cube as defined in Appendix A. The following collection of procedures was written by Ken Monks.

with(group):
# The generators of the Cube Group
M:=[[1,2],[3,4],[5,6],[7,8]]:
R:=[[1,2,4,3],[5,6,8,7]]:
S:=[[1,3,7,5],[2,4,8,6]]:
T:=[[1,5,6,2],[3,7,8,4]]:
E:=[]:
# This routine allows mulperms to work with a list of more than two
# elements at time composing left to right.
Mult:=proc()
local a,p:
a:=[];
for p in args do a:=mulperms(a,p) od:
a;
end:
# lookup table to name the elements according to the notation in Appendix A.
Names:=[
[E,e],[M,-e],
[R,r],[Mult(M,R),-r],
[Mult(R,R),r2],[Mult(M,R,R),-r2],
[Mult(R,R,R),r3],[Mult(M,R,R,R),-r3],
[S,s],[Mult(M,S),-s],
[Mult(R,S),rs],[Mult(M,R,S),-rs],
[Mult(R,R,S),r2s],[Mult(M,R,R,S),-r2s],
[Mult(R,R,R,S),r3s],[Mult(M,R,R,R,S),-r3s],
[Mult(S,S),s2],[Mult(M,S,S),-s2],
[Mult(R,S,S),rs2],[Mult(M,R,S,S),-rs2],
[Mult(R,R,S,S),r2s2],[Mult(M,R,R,S,S),-r2s2],
[Mult(R,R,R,S,S),r3s2],[Mult(M,R,R,R,S,S),-r3s2],
[Mult(S,S,S),s3],[Mult(M,S,S,S),-s3],
[Mult(R,S,S,S),rs3],[Mult(M,R,S,S,S),-rs3],
[Mult(R,R,S,S,S),r2s3],[Mult(M,R,R,S,S,S),-r2s3],
[Mult(R,R,R,S,S,S),r3s3],[Mult(M,R,R,R,S,S,S),-r3s3],
[T,t],[Mult(M,T),-t],
[Mult(R,T),rt],[Mult(M,R,T),-rt],
[Mult(R,R,T),r2t],[Mult(M,R,R,T),-r2t],
[Mult(R,R,R,T),r3t],[Mult(M,R,R,R,T),-r3t],

[Mult(T,T,T),t3],[Mult(M,T,T,T),-t3],
[Mult(R,T,T,T),rt3],[Mult(M,R,T,T,T),-rt3],
[Mult(R,R,T,T,T),r2t3],[Mult(M,R,R,T,T,T),-r2t3],
[Mult(R,R,R,T,T,T),r3t3],[Mult(M,R,R,R,T,T,T),-r3t3],
[''Error'',ERROR] ]:
# To use a permutation as a function
Apply:=proc(p,n)
local a,i,j:
a:=n;
for i from 1 to nops(p) do
   if has(p[i],a) then
      for j to nops(p[i]) while p[i][j]<>a do
      od:
      if j=nops(p[i]) then
         a:=p[i][1]
      else a:=p[i][j+1]
      fi;
   fi:
od:
a
end:
# to determine if two permutations are the same element of the group
Equal:=proc(n,p,q)
local i;
for i to n while Apply(p,i)=Apply(q,i) do od:
if i>n then true else false fi:
end:
# to find the name of a permutation in the group
PermToName:=proc(a)
local i; global Names:
for i in Names while not Equal(8,i[1],a) do od;
i[2];
end:
# to find the permutation associated with a name in the group
NameToPerm:=proc(n)
global Names:
select(has,Names,n)[1][1];
end:
# To multiply the elements using the names defined above instead

# of permutation notation
Times:=proc(a,b)
   PermToName(mulperms(NameToPerm(a),NameToPerm(b)))
end:

References

[1] Monks, K. Chaos: A Maple package for Chaos and Fractals. version 3.16. 2000. on-line, www.scranton.edu/~monks

[2] Lay, David C. Linear Algebra and Its Applications. New York: Addison-Wesley Publishing Company, 1994.

[3] Barnsley, M. Fractals Everywhere. San Diego: Academic Press, 1988.

[4] Peitgen, Heinz-Otto, Hartmut Jürgens, and Dietmar Saupe. Chaos and Fractals: New Frontiers of Science New York: Springer-Verlag, 1992.

152 E Elm St. Apt 4, Dunmore, PA 18512

   E-mail address: riggic2@scranton.edu
   URL: http://facweb.uofs.edu/~monks/fsrp.html