Clear[Ktile]

Off[General::spell1]

A Translation of
"The Construction of Self-Similar Tiles",
Richard Kenyon, into Mathematica
Submitted to Math Source by Roger L. Bagula 30 April 2005

http://www.math.ubc.ca/~kenyon/papers/index.html
The construction of self-similar tilings , Geom. and Func. Analysis 6,(1996):417-488. Thurston showed that the expansion constant of a self-similar tiling of the plane must be a complex Perron number (algebraic integer strictly larger in modulus than its Galois conjugates except for its complex conjugate). Here we prove that, conversely, for every complex Perron number there exists a self-similar tiling. We also classify the expansion constants for self-similar tilings which have a rotational symmetry of order n.
Section 6: page 15-18

"A complex Perron number is an algebraic integer which is strictly larger in modulus than its Galois conjugates (the other roots of its minimal polynomial), except for its complex conjugate."

List of operations performed in the calculations to pick best roots:
1) Turn polynomial Rule into List,
2) Get rid of λ symbols,
3) Get rid of Real solutions,
4) Discard roots with negative Arg, they just give mirror image tiles of those with positive Arg
5) Get rid of empty lists
6) Get rid of cases with no complex root remaining
7) Get rid of roots with Modulus less than the root with greatest Modulus, and with multiple roots having same greatest modulus

Kenyon' s paper on the web.

You are advised to evaluate each subsection  separately as the amount of output is large for each section.

compute a tile from a set of polynomial specifications

\(Ktile[n_, p_, q_, r_, λ_] := 
      Module[{}, 
Print[n, ",", p, ",", q, ",", 
          r, ",", λ]; 
 (*\ 
          create\ the\ rules\ that\ generate\ the\ edge\ of\ the\ tile\ *) \ 
phi[a] = b; phi[\(-a\)] = \ \(-b\); phi[b] = c; 
        phi[\(-b\)] = \ \(-c\); 
phi[c] = 
          Flatten[{Table[c, {p}], Table[\(-a\), {r}], 
              Table[\(-b\), {q}]}]; 
phi[\(-c\)] = 
          Flatten[{Table[b, {q}], Table[a, {r}], 
              Table[\(-c\), {p}]}]; 
phi[s_List] := 
          FixedPoint[
            Function[z, 
              Flatten[Split[
                      z, #1 === \(-#2\) &] /. 
{_, _} :> Sequence[]] /. {{\(-\((x_)\)\), m___, 
                      x_} :> {m}, {x_, 
                      m___, \(-\((x_)\)\)} :> {m}}
], 
            Flatten[Map[phi, Flatten[s]]]]; 
commutator[
            x_Symbol, y_Symbol] := 
          Flatten[{x, y, \(-x\), \(-y\)}]; 
commutator[
            x_Symbol, y_List] := 
          Flatten[{x, 
              y, \(-x\), \(-Reverse[y]\)}]; 
commutator[
            x_List, y_Symbol] := 
          Flatten[{x, 
              y, \(-Reverse[x]\), \(-y\)}]; 
commutator[
            x_List, y_List] := 
          Flatten[{x, 
              y, \(-Reverse[x]\), \(-Reverse[y]\)}]; 
 (*\ 
          generate\ the\ list\ of\ vectors\ a, b, 
          c, \(-a\), \(-b\), \(-c\)\ that\ make\ up\ the\ tile\ *) \ 
edges[i_] := 
          Nest[phi, commutator[b, c], i]; 
 (*\ 
          turn\ edges\ into\ verticies\ *) 
verticies[i_] := 
          FoldList[Plus, 0, edges[i]]; 
 (*\ turn\ a, b, 
          c\ vectors\ into\ 1, λ, λ^2\ vectors\ and\ then\ substitute\ numeric\ value\ for\ λ\ *) 
lambdaMap[
            i_] := \(\(Map[#/λ^\((i - 1)\) &, verticies[i]] /. 
                a -> 1\) /. b -> λ\) /. 
            c -> λ^2; 
 (*\ 
          split\ a\ complex\ value\ into\ a\ list\ of\ the\ real\ and\ complex\ coefficients, \ 
          for\ ListPlot\ *) 
splitComplex[z_] := {Re[z], 
            Im[z]}; 
Timing[
          Table[ListPlot[Map[splitComplex, lambdaMap[i]], 
              PlotJoined -> True, Axes -> False], {i, 11, 
              11}]]
]\)

pick best roots

CherryPicker[{n_, p_, q_, r_, roots_}] := 
      Module[{s}, 
s = 
          Sort[roots, 
            OrderedQ[{Abs[#1], Abs[#2]}] &]; 
While[Abs[
                s[[1]]][[1]] < Abs[s[[-1]]][[1]], 
          s = Drop[s, 1]]; 
If[
          Length[s] > 1, {}, {n, p, q, 
            r, s[[1]][[1]]}]]

Kenyonarray = 
        Flatten[Table[{n, p, q, r, 
              NSolve[λ^n - p λ^(n - 1) + q λ + 
                    r == 0, λ]}, {n, 3, 6}, {p, 0, 3}, {q, 0, 
              3}, {r, 1, 4}], 3];

refine the array

Kenyonarray2 = Kenyonarray /. 
                    Rule -> List /. 
List[λ, 
                      x_] -> 
                    x /. 
List[_Real] -> 
                  Sequence[] /. 
List[x_Complex] -> 
                If[Arg[x] < 0, List[], 
                  List[x]] /. 
List[] -> 
              Sequence[] /. 
List[_Integer, _Integer, _Integer, _Integer] -> Sequence[];

Kenyonmap = Map[CherryPicker, Kenyonarray2];

apply the tile function to the best complex Perron numbers found

Apply[Ktile, Kenyonmap, 2];

new substitution tile computed from a set of polynomial specifications

\(Poly2tile[n_, p_, q_, r_, λ_] := 
      Module[{}, 
Print[n, ",", p, ",", q, ",", 
          r, ",", λ]; 
 (*\ 
          create\ the\ rules\ that\ generate\ the\ edge\ of\ the\ tile\ *) \ 
phi[a] = b; phi[\(-a\)] = \ \(-b\); phi[b] = c; 
        phi[\(-b\)] = \ \(-c\); 
phi[c] = d; 
        phi[\(-c\)] = \ \(-d\); phi[d] = e; 
        phi[\(-d\)] = \ \(-e\); 
phi[e] = 
          Flatten[{Table[e, {p}], Table[\(-a\), {r}], Table[\(-b\), {1}], 
              Table[\(-d\), {q}]}]; 
phi[\(-e\)] = 
          Flatten[{Table[d, {q}], Table[b, {1}], Table[a, {r}], 
              Table[\(-e\), {p}]}]; 
phi[s_List] := 
          FixedPoint[
            Function[z, 
              Flatten[Split[
                      z, #1 === \(-#2\) &] /. 
{_, _} :> Sequence[]] /. {{\(-\((x_)\)\), m___, 
                      x_} :> {m}, {x_, 
                      m___, \(-\((x_)\)\)} :> {m}}
], 
            Flatten[Map[phi, Flatten[s]]]]; 
commutator[
            x_Symbol, y_Symbol] := 
          Flatten[{x, y, \(-x\), \(-y\)}]; 
commutator[
            x_Symbol, y_List] := 
          Flatten[{x, 
              y, \(-x\), \(-Reverse[y]\)}]; 
commutator[
            x_List, y_Symbol] := 
          Flatten[{x, 
              y, \(-Reverse[x]\), \(-y\)}]; 
commutator[
            x_List, y_List] := 
          Flatten[{x, 
              y, \(-Reverse[x]\), \(-Reverse[y]\)}]; 
 (*\ 
          generate\ the\ list\ of\ vectors\ a, b, 
          c, \(-a\), \(-b\), \(-c\)\ that\ make\ up\ the\ tile\ *) \ 
edges[i_] := 
          Nest[phi, commutator[d, e], i]; 
 (*\ 
          turn\ edges\ into\ verticies\ *) 
verticies[i_] := 
          FoldList[Plus, 0, edges[i]]; 
 (*\ turn\ a, b, 
          c\ vectors\ into\ 1, λ, λ^2\ vectors\ and\ then\ substitute\ numeric\ value\ for\ λ\ *) 
lambdaMap[
            i_] := \(\(\(\(Map[#/λ^\((i - 1)\) &, verticies[i]] /. 
                    a -> 1\) /. b -> λ\) /. 
                c -> λ^2\) /. d -> λ^3\) /. 
            e -> λ^4; 
 (*\ 
          split\ a\ complex\ value\ into\ a\ list\ of\ the\ real\ and\ complex\ coefficients, \ 
          for\ ListPlot\ *) 
splitComplex[z_] := {Re[z], 
            Im[z]}; 
Timing[
          Table[ListPlot[Map[splitComplex, lambdaMap[i]], 
              PlotJoined -> True, Axes -> False], {i, 11, 
              11}]]
]\)

pick best roots

CherryPicker[{n_, p_, q_, r_, roots_}] := 
      Module[{s}, 
s = 
          Sort[roots, 
            OrderedQ[{Abs[#1], Abs[#2]}] &]; 
While[Abs[
                s[[1]]][[1]] < Abs[s[[-1]]][[1]], 
          s = Drop[s, 1]]; 
If[
          Length[s] > 1, {}, {n, p, q, 
            r, s[[1]][[1]]}]]

Make array of roots of the new polynomial  based on  power 5

Poly2array = 
        Flatten[Table[{n, p, q, r, 
              NSolve[λ^n - p λ^(n - 2) + q λ + 
                    r == 0, λ]}, {n, 3, 6}, {p, 0, 3}, {q, 0, 
              3}, {r, 0, 4}], 3];

Poly2array2 = Poly2array /. 
                    Rule -> List /. 
List[λ, 
                      x_] -> 
                    x /. 
List[_Real] -> 
                  Sequence[] /. 
List[x_Complex] -> 
                If[Arg[x] < 0, List[], 
                  List[x]] /. 
List[] -> 
              Sequence[] /. 
List[_Integer, _Integer, _Integer, _Integer] -> Sequence[];

Poly2map = Map[CherryPicker, Poly2array2];

apply the tile function to the best complex Perron numbers found

Apply[Poly2tile, Poly2map, 2];

new substitution tile computed from a set of polynomial specifications

\(Poly3tile[n_, p_, q_, r_, λ_] := 
      Module[{}, 
Print[n, ",", p, ",", q, ",", 
          r, ",", λ]; 
 (*\ 
          create\ the\ rules\ that\ generate\ the\ edge\ of\ the\ tile\ *) \ 
phi[a] = b; phi[\(-a\)] = \ \(-b\); phi[b] = c; 
        phi[\(-b\)] = \ \(-c\); 
phi[c] = d; 
        phi[\(-c\)] = \ \(-d\); phi[d] = e; 
        phi[\(-d\)] = \ \(-e\); 
phi[e] = f; 
        phi[\(-e\)] = \ \(-f\); phi[f] = g; 
        phi[\(-f\)] = \ \(-g\); 
phi[g] = 
          Flatten[{Table[g, {p}], Table[\(-e\), {r}], 
              Table[\(-f\), {q}]}]; 
phi[\(-g\)] = 
          Flatten[{Table[f, {q}], Table[e, {r}], 
              Table[\(-g\), {p}]}]; 
phi[s_List] := 
          FixedPoint[
            Function[z, 
              Flatten[Split[
                      z, #1 === \(-#2\) &] /. 
{_, _} :> Sequence[]] /. {{\(-\((x_)\)\), m___, 
                      x_} :> {m}, {x_, 
                      m___, \(-\((x_)\)\)} :> {m}}
], 
            Flatten[Map[phi, Flatten[s]]]]; 
commutator[
            x_Symbol, y_Symbol] := 
          Flatten[{x, y, \(-x\), \(-y\)}]; 
commutator[
            x_Symbol, y_List] := 
          Flatten[{x, 
              y, \(-x\), \(-Reverse[y]\)}]; 
commutator[
            x_List, y_Symbol] := 
          Flatten[{x, 
              y, \(-Reverse[x]\), \(-y\)}]; 
commutator[
            x_List, y_List] := 
          Flatten[{x, 
              y, \(-Reverse[x]\), \(-Reverse[y]\)}]; 
 (*\ 
          generate\ the\ list\ of\ vectors\ a, b, 
          c, \(-a\), \(-b\), \(-c\)\ that\ make\ up\ the\ tile\ *) \ 
edges[i_] := 
          Nest[phi, commutator[f, g], i]; 
 (*\ 
          turn\ edges\ into\ verticies\ *) 
verticies[i_] := 
          FoldList[Plus, 0, edges[i]]; 
 (*\ turn\ a, b, 
          c\ vectors\ into\ 1, λ, λ^2\ vectors\ and\ then\ substitute\ numeric\ value\ for\ λ\ *) 
lambdaMap[
            i_] := \(\(\(\(\(\(Map[#/λ^\((i - 1)\) &, verticies[i]] /. 
                        a -> 1\) /. b -> λ\) /. 
                    c -> λ^2\) /. d -> λ^3\) /. 
                e -> λ^4\) /. f -> λ^5\) /. 
            g -> λ^6; 
 (*\ 
          split\ a\ complex\ value\ into\ a\ list\ of\ the\ real\ and\ complex\ coefficients, \ 
          for\ ListPlot\ *) 
splitComplex[z_] := {Re[z], 
            Im[z]}; 
Timing[
          Table[ListPlot[Map[splitComplex, lambdaMap[i]], 
              PlotJoined -> True, Axes -> False], {i, 11, 
              11}]]
]\)

pick best roots

CherryPicker[{n_, p_, q_, r_, roots_}] := 
      Module[{s}, 
s = 
          Sort[roots, 
            OrderedQ[{Abs[#1], Abs[#2]}] &]; 
While[Abs[
                s[[1]]][[1]] < Abs[s[[-1]]][[1]], 
          s = Drop[s, 1]]; 
If[
          Length[s] > 1, {}, {n, p, q, 
            r, s[[1]][[1]]}]]

Make array of roots of the new polynomial  based on  power 7

Poly3array = 
        Flatten[Table[{n, p, q, r, 
              NSolve[λ^n - p λ^(n - 4) + q λ + 
                    r == 0, λ]}, {n, 6, 8}, {p, 0, 3}, {q, 0, 
              3}, {r, 0, 4}], 3];

Poly3array2 = Poly3array /. 
                    Rule -> List /. 
List[λ, 
                      x_] -> 
                    x /. 
List[_Real] -> 
                  Sequence[] /. 
List[x_Complex] -> 
                If[Arg[x] < 0, List[], 
                  List[x]] /. 
List[] -> 
              Sequence[] /. 
List[_Integer, _Integer, _Integer, _Integer] -> Sequence[];

Poly3map = Map[CherryPicker, Poly3array2];

apply the tile function to the best complex Perron numbers found

Apply[Poly3tile, Poly3map, 2];

Clear[g1, g2, tiling, tile1, tile2, r, g, b, f]

Kenyon tiles calculated

tile1[n_, p_, q_, r_, λ_] := 
      Module[{}, 

 (* 
          create the rules that generate the edge of the tile *)  
phi[a] = b; phi[-a] =  -b; phi[b] = c; 
        phi[-b] =  -c; 
phi[c] = 
          Flatten[{Table[c, {p}], Table[-a, {r}], 
              Table[-b, {q}]}]; 
phi[-c] = 
          Flatten[{Table[b, {q}], Table[a, {r}], 
              Table[-c, {p}]}]; 
phi[s_List] := 
          FixedPoint[
            Function[z, 
              Flatten[Split[
                      z, #1 === -#2 &] /. 
{_, _} :> Sequence[]] /. {{-(x_), m___, 
                      x_} :> {m}, {x_, 
                      m___, -(x_)} :> {m}}
], 
            Flatten[Map[phi, Flatten[s]]]]; 
commutator[
            x_Symbol, y_Symbol] := 
          Flatten[{x, y, -x, -y}]; 
commutator[
            x_Symbol, y_List] := 
          Flatten[{x, 
              y, -x, -Reverse[y]}]; 
commutator[
            x_List, y_Symbol] := 
          Flatten[{x, 
              y, -Reverse[x], -y}]; 
commutator[
            x_List, y_List] := 
          Flatten[{x, 
              y, -Reverse[x], -Reverse[y]}]; 
 (* 
          generate the list of vectors a, b, 
          c, -a, -b, -c that make up the tile *)  
edges[i_] := 
          Nest[phi, commutator[b, c], i]; 
 (* 
          turn edges into verticies *) 
verticies[i_] := 
          FoldList[Plus, 0, edges[i]]; 
 (* turn a, b, 
          c vectors into 1, λ, λ^2 vectors and then substitute numeric value for λ *) 
lambdaMap[
            i_] := Map[#/λ^(i - 1) &, verticies[i]] /. 
                a -> 1 /. b -> λ /. 
            c -> λ^2; 
 (* 
          split a complex value into a list of the real and complex coefficients,  
          for ListPlot *) 
splitComplex[z_] := {Re[z], 
            Im[z]}; 
Timing[
          Table[ListPlot[Map[splitComplex, lambdaMap[i]], 
              PlotJoined -> True, Axes -> False, 
              PlotStyle -> {RGBColor[0, 0, 1]}], {i, 14, 
              14}]]
]

Kenyon tiles with first substitution set to zero: fractal subset, one of three possible

tile2[n_, p_, q_, r_, λ_] := 
      Module[{}, 

 (* 
          create the rules that generate the edge of the tile *)  
phi[a] = b; phi[-a] =  -b; phi[b] = c; 
        phi[-b] =  -c; 
phi[c] = 
          Flatten[{Table[c, {p}], Table[-a, {r}], 
              Table[-b, {q}]}]; 
phi[-c] = 
          Flatten[{Table[b, {q}], Table[a, {r}], 
              Table[-c, {p}]}]; 
phi[s_List] := 
          FixedPoint[
            Function[z, 
              Flatten[Split[
                      z, #1 === -#2 &] /. 
{_, _} :> Sequence[]] /. {{-(x_), m___, 
                      x_} :> {m}, {x_, 
                      m___, -(x_)} :> {m}}
], 
            Flatten[Map[phi, Flatten[s]]]]; 
commutator[
            x_Symbol, y_Symbol] := 
          Flatten[{x, y, -x, -y}]; 
commutator[
            x_Symbol, y_List] := 
          Flatten[{x, 
              y, -x, -Reverse[y]}]; 
commutator[
            x_List, y_Symbol] := 
          Flatten[{x, 
              y, -Reverse[x], -y}]; 
commutator[
            x_List, y_List] := 
          Flatten[{x, 
              y, -Reverse[x], -Reverse[y]}]; 
 (* 
          generate the list of vectors a, b, 
          c, -a, -b, -c that make up the tile *)  
edges[i_] := 
          Nest[phi, commutator[b, c], i]; 
 (* 
          turn edges into verticies *) 
verticies[i_] := 
          FoldList[Plus, 0, edges[i]]; 
 (* turn a, b, 
          c vectors into 1, λ, λ^2 vectors and then substitute numeric value for λ *) 
lambdaMap[
            i_] := Map[#/λ^(i - 1) &, verticies[i]] /. 
                a -> 0 /. b -> λ /. 
            c -> λ^2; 
 (* 
          split a complex value into a list of the real and complex coefficients,  
          for ListPlot *) 
splitComplex[z_] := {Re[z], 
            Im[z]}; 
Timing[
          Table[ListPlot[Map[splitComplex, lambdaMap[i]], 
              PlotJoined -> True, Axes -> False, 
              PlotStyle -> {RGBColor[1, 0, 0]}], {i, 14, 
              14}]]
]

CherryPicker[{n_, p_, q_, r_, roots_}] := 
      Module[{s}, 
s = 
          Sort[roots, 
            OrderedQ[{Abs[#1], Abs[#2]}] &]; 
While[Abs[
                s[[1]]][[1]] < Abs[s[[-1]]][[1]], 
          s = Drop[s, 1]]; 
If[
          Length[s] > 1, {}, {n, p, q, 
            r, s[[1]][[1]]}]]

A complex Perron number is an algebraic integer which is strictly larger in modulus than its Galois conjugates (the other roots of its
minimal polynomial), except for its complex conjugate.

aa = 
        Flatten[Table[{n, p, q, r, 
              NSolve[λ^n - p λ^(n - 1) + q λ + 
                    r == 0, λ]}, {n, 4, 4}, {p, 0, 2}, {q, 0, 
              2}, {r, 1, 2}], 3];

Turn Rule into List,
Get rid of λ symbols,
Get rid of Real solutions,
Discard roots with negative Arg, they just give mirror image tiles of those with positive Arg
Get rid of empty lists
Get rid of cases with no complex root remaining

bb = aa /. 
                    Rule -> List /. 
List[λ, 
                      x_] -> 
                    x /. 
List[_Real] -> 
                  Sequence[] /. 
List[x_Complex] -> 
                If[Arg[x] < 0, List[], 
                  List[x]] /. 
List[] -> 
              Sequence[] /. 
List[_Integer, _Integer, _Integer, _Integer] -> Sequence[];

Now get rid of roots with Modulus less than the root with greatest Modulus, and with multiple roots having same greatest modulus

cc = Map[CherryPicker, bb]

Discard empty lists

% /. List[] -> Sequence[] // TableForm

Length[cc]

plots showing tile and fractal superimposed on the same scale

Table[
      Show[{Flatten[Apply[tile1, {cc[[i]]}, 2]][[2]], Flatten[
              Apply[tile2, {cc[[i]]}, 2]][[2]]}], {i, 3, 
        Length[cc]}]

The way I got the ideas for new higher substitutions and polynomials  is hard to put together. 
    I had done work on Pisot Theta[0], Theta[1] and Theta[2] trying to get  IFS tiles.
    S.R. Hinsley claimed to have a  Kenyon tile IFS. I spent several months trying to duplicate that.
    I finally went to the original article and tried from  there. Up until lately I've had little luck with substitution fractals and tiles.

In probably totally unrelated work on minimal surfaces ,
CMC ( constant mean curvature surfaces),
and solitons, I went into the theory of elliptical equations pretty deeply.
If we make the Kenyon polynomial into:
(dx/dt)2=x^n-p*x^(n-1)+q*x+r
Then these algebraic numbers are associated with the cycles of the
lattice that the elliptic equation establishes on the plane.
That's why the degenerates are lines or quadrilaterals, I suppose.
It stills seems strange that this type of cycling behavior should give minimal substitutions
and tiles. I had thought of using Weierstrass elliptical equations in
Weierstrass fractal functions before as Besicovitch-Ursell like functions,
but I doubt that would give tiles. The two numbers {g2,g3}
map to the complex Perron numbers, I think , for n=3.
It gets harder for higher n, but seems to be related to a lattice structure of:
{a,b,c}-> -{a,b,c}
The Complex Perron number
z0=x0+I*y0
acts as the self-similar scaling/ affine scaling of the tiles.
In IFS terms that comes out like
S=scale number=x0/(x02+y02)-I*y0/(x02+y02)
z'=S*z+a(i)
a(i) are the digit set of the tile. I used that in the Reptiles I constructed.
This result isn't a "total" answer, but it makes  the tiles make more sense
and my polynomial which is mostly hyper elliptic with substitution symmetry:
{a,b,c,d,e}-->- {a,b,c,d,e}
Which is the next Prime Galois field level.
The Reptiles are the quadratic two symbol version of this:
{a,b}->- {a,b}
So we have levels: 2,3,5
The next polynomial of Kenyon in Galois field terms would/should  be:
{a,b,c,d,e,f,g}-> - {a,b,c,d,e,f,g}
w2=x^n-p*x^(n-4)+q*x+r
At higher levels it appears that negative q at least will give tiles ( I got one last night
from my elliptical invariant analysis of polynomial).
The higher level seem to roughly correspond to Theta(0),Theta(1) and Theta(2)
levels of minimal symmetry.

On the Fractal subsets found for Kenyon  tiles using Mauldin-Williams theory and null substitutions:
{1,λ, λ^2}->{1,λ,0},{1, 0, λ^2},{0,λ, λ^2}
Essentially this result is like replacing one symbol with a {} null substitution. 
Any substitution can be represented as a state machine digraph of the Mauldin-Williams type.
Removing one vertex  seems to result in tiles producing fractals.

Some References:
        Related paper link to Perron Numbers:
            http://www.google.com/url?sa=U&start=1&q=http://www.mrlonline.org/mrl/2004-011-003/2004-011-003-001.pdf&e=10187
Or
http://www.mrlonline.org/mrl/2004-011-003/2004-011-003-001.pdf
On Mauldin -Williamd Digraphs:
      "On the Hausdorff Dimension of Some Graphs," (with S.C.Williams),Transactions American Mathematical Society 298 (2) (1986),793-803.
http://www.math.unt.edu/~mauldin/papers/no60.pdf
"Hausdorff Dimension in Graph Directed Constructions," (with S.C.Williams),Transactions American Mathematical Society,309(1988),811-829.
http://www.math.unt.edu/~mauldin/papers/no67.pdf