MATHEMATICA IN ACTION

The Traveling Salesman and the Turtle

Combining a two-dimensional turtle with Mathematica’s three-dimensional capabilities gives us a new way to look at space-filling curves.

Turtle.m package and code for figures to be added to the Electronic Supplement

The Turtle.m Package

Off[General::spell1];

Remove[FractalTurtle, StartDirection, TerminalSubstitution, TurtleStep, 
StartPosition,  TurtlePath]

  
BeginPackage["Turtle`", {"Utilities`FilterOptions`"}];

CurveDepth::usage = "CurveDepth is an option to TravelingSalesmanTour that specifies the depth of the space-filling curve approximation to use when estimating the position of points on the curve.";

FractalTurtle::usage = "FractalTurtle[r, s, d, a] creates (in \"TurtlePath\") and shows the list of points that a turtle would make if it starts from s and then applies the rewrite rule d times. The size of the turns is set by a (assumed to be in degrees). Strings can involve F, B, +, -, and i, which instruct the turtle to move forward, move back, turn left, turn right, and roll over, respectively. The characters [ (push) and ] (pop) can also be used to store a location in a stack and pop back to it.";

ShowPointLabels::usage = "ShowPointLabels is an option to TravelingSalesmanTour that causes the points to be large, with labels inside them.";

LineStyle::usage = "LineStyle is an option to TravelingSalesmanTour that specifies the style of the line forming the tour.";

PointSwitching::usage = "PointSwitching is an option to TravelingSalesmanTour that follows up the space-filling heuristic with a three passes of a check on switches of consecutive pairs and triples to find improvements.";

PointStyle::usage = "PointStyle is an option to TravelingSalesmanTour that specifies the style of the points.";

RandomPoints::usage = "RandomPoints[n] gives n random points in the unit square.";

ShowDistance::usage = "ShowDistance is an option to TravelingSalesmanTour that asks the distance of the tour to be used as a plot label.";

ShowLabeledPoints::usage = "ShowLabeledPoints[p] shows the points in p in a large size with random soft colors and the labels inside the points.";

ShowIndices::usage = "ShowIndices is an option to TravelingSalesmanTour that causes the index list to be returned.";

StartPosition::usage = "StartPosition is an option to FractalTurtle that specifies the point at which the turtle starts its journey.";

StartDirection::usage = "StartDirection is an option to FractalTurtle that specifies the initial direction the turtle faces. It must be a unit vector.";

TerminalSubstitution::usage = "TerminalSubstitution is an option to FractalTurtle that causes one additional substitution (or many if it is set to be a list) to be made after the specified depth is reached.";

TravelingSalesmanTour::usage = "TravelingSalesmanTour[p] uses a Hilbert space-filling curve to come up with a traveling salesman tour through the given points that is, on average, not too far from optimal.";

TSPDistance::usage = "TSPDistance[data, tour] gives the sum of the straight-line lengths corresponding to a trip through the points in data in the order determined by tour. The last entry in tour must be the same as the first.";

TurtleStep::usage = "TurtleStep is an option to FractalTurtle that specifies the length of each turtle step.";

TurtlePath::usage = "TurtlePath is where the points of the turtle's path are stored by FractalTurtle.";

TSPDistance::bddata = "A traveling salesman tour must have its last point being the same as its first point.";

SecondaryStepSize::usage = "SecondaryStepSize is an option to FractalTurtle that specifies the step size for a secondary set of moves, corrsponding to the instructions \"f\" and \"b\".";

Begin["`Private`"];

Options[FractalTurtle] = {StartPosition -> {0, 0},
  StartDirection -> {1, 0}, 
  TurtleStep -> 1, TerminalSubstitution -> None,
  PlotStyle -> AbsoluteThickness[0.4]};

FractalTurtle[r_, s_String, depths_List, angle_, 
  opts___] := Show[GraphicsArray[
   Map[FractalTurtle[r,s,#,angle,opts,
     DisplayFunction->Identity]&, depths, {Depth[depths]-1}]]]

FractalTurtle[rewrite_, start_String, depth_Integer, 
  angle_, opts___] := Module[{
   turtle, sdir, scale, i = 1, 
   rotateleft, rotateright, stack,gg},
 inds = {1};
 {startingPt, dir, ts, term, sty, secsize} = 
   {StartPosition, StartDirection, TurtleStep, 
     TerminalSubstitution, PlotStyle, SecondaryStepSize} /.
     {opts} /. Options[FractalTurtle];
 TurtlePath = {last = startingPt};
 stack = {{startingPt, dir}};
 rotateleft = N[
   {{Cos[angle Degree], -Sin[angle Degree]}, 
    {Sin[angle Degree], Cos[angle Degree]}}];
 rotateright = Transpose[rotateleft];
 If[!ListQ[sty], sty = {sty} /. {Automatic} ->{}];

 s = Nest[StringReplace[#, rewrite]&, start, depth];
If[term =!= None, s = StringReplace[s, term]];




 Characters[s] /. {
  "+" :> (dir = rotateleft . dir;),
  "-" :> (dir = rotateright . dir;), 
  "F" :> (i++;TurtlePath = 
    {TurtlePath, last += ts*dir};),
  "B" :> (i++;TurtlePath = 
     {TurtlePath,last -= ts*dir};),
  "f" :> (i++;TurtlePath = 
    {TurtlePath, last += secsize*dir};),
  "b" :> (i++;TurtlePath = 
     {TurtlePath,last -= secsize*dir};),
  "i" :> ({rotateleft, rotateright} = {rotateright, rotateleft};),
  "[" :> (stack = Flatten[{{{last, dir}},stack},1];),
  "]" :> (i++; dir = stack[[1,2]]; TurtlePath = 
     {TurtlePath, last=stack[[1,1]]}; stack = Rest[stack];
      inds = Flatten[{inds,i}];)};
      
TurtlePath = Partition[Flatten[TurtlePath], 2];

   gg = Graphics[If[inds == {1},
    {Append[sty, Line[TurtlePath]]}, 
TurtlePath = (TurtlePath[[Range[#[[1]], #[[2]]-1]]] & /@ 
  Partition[inds, 2, 1]);
Join[sty, Line /@ TurtlePath]],

   FilterOptions[Graphics, opts], 
   AspectRatio -> Automatic,
   PlotRange -> All];
   
{{x0,x1}, {y0,y1}} = PlotRange[gg];

 Show[gg,
  Evaluate@FilterOptions[Graphics, opts],
    PlotRange -> {{x0,x1}+{-.1, .1} (x1-x0),
                  {y0 ,y1} + {-.1, .1} (y1-y0)},
    Frame->True, FrameTicks->None]]

dC = Compile[{x0,y0,x1,y1},Sqrt[(x0-x1)^2+(y0-y1)^2]];

distance[points_] := Plus @@ 
  Apply[dC, Flatten /@ Partition[points,2,1],{1}]

TSPDistance[data_, tour_] := (If[Last[tour] != First[tour],
Message[TSPDistance::bddata]; Return[HoldForm[TSPDistance[tour]]]];
distance[data[[tour]]])

position[{x_, y_}, k_Integer] := Which[
  x < 0.5 && y < 0.5, 
	Mod[0.875 + position[{1-2 x, 1- 2y}, k-1]/4, 1],
  x < 0.5 && y >= 0.5, 
      0.125 + position[{2y-1,1-2x},k-1]/4,
  x >= 0.5 && y >= 0.5, 
    0.375 + position[ {2x-1,2y-1},k-1]/4,
  True, 0.625 + position[{1-2y,2x-1},k-1]/4 ]

position[{x_, y_}, 0] := Which[
  x < 0.5 && y < 0.5, 	 0,
  x < 0.5  && y >= 0.5,  0.25,
  x >= 0.5 && y > 0.5,   0.5,
  True, 				 0.75]

tsp[pointList_, k_Integer] := First@Thread@Sort[
 Map[{#, position[#, k]}&, pointList],
  #1[[2]] < #2[[2]] &]
			
improve[path_] := Module[{
 path1 = path, len = Length[path]},
 adj[s_] := Mod[s,len] /. {0->len}; 
 Do[If[distance[path1[[adj@{i-1,i,i+1,i+2}]]]> 
       distance[path1[[adj@{i-1,i+1,i,i+2}]]],
  path1 = ReplacePart[ReplacePart[path1, 
     path1[[adj[i+1]]],i], path1[[i]],adj[i+1]]],
    {i, len}];
 Do[If[distance[path1[[adj@Range[i-1,i+3]]]]> 
    distance[path1[[adj@{i-1,i+2,i+1,i,i+3}]]], 
  path1 = ReplacePart[ReplacePart[path1, 
    path1[[adj[i+2]]],i], path1[[i]],adj[i+2]]], 
    {i, len}];
  path1]
			
Options[TravelingSalesmanTour] = {
 CurveDepth -> 10, PointSwitching->True, 
 ShowDistance->True, PointStyle->PointSize[0.015], 
 LineStyle -> {}, ShowIndices -> False, 
 ShowPointLabels->False};

TravelingSalesmanTour[pts_, opts___] := Module[
 {dep, impQ, disQ, tour, tourQ},
 {dep, impQ, disQ, lsty, psty, tourQ, labQ} = 
   {CurveDepth, PointSwitching, ShowDistance,
   LineStyle, PointStyle, ShowIndices,ShowPointLabels} /. 
      {opts} /. Options[TravelingSalesmanTour];
 If[!ListQ[lsty], lsty = {lsty}];
 If[!ListQ[psty], psty = {psty}];
 tour = tsp[pts,dep];
 If[impQ, tour = Nest[improve,tour,3]];

 ss=Show[Graphics[{Thickness[0.0005],
Append[lsty, Line[realtour=Append[tour, First[tour]]]],
	If[labQ, {{AbsolutePointSize[16], 
  {RGBColor @@ Array[Random[Real, {0.4,1}]&,3], 
     Point[#]}& /@ pts},
     MapIndexed[Text[#2[[1]], #1+{.005,0}]&, pts]},
	Append[psty, Point /@ pts]]}],

  Evaluate[FilterOptions[Graphics, opts]],
  PlotLabel->If[disQ, 
  StringForm["Length is ``",
 NumberForm[distance[realtour],6]], None],
  PlotRange->{{-0.1,1.1},{-0.1,1.1}}, 
  AspectRatio->Automatic,
  DefaultFont->{"Times-Bold",12}];
  If[!tourQ, ss, {ss, Flatten[
   Position[pts, #]& /@ Append[tour, First[tour]]]}]]

RandomPoints[n_] := Array[{Random[], Random[]}&, n]

ShowLabeledPoints[pts_] := Show[Graphics[{
{AbsolutePointSize[16], 
  {RGBColor @@ Array[Random[Real, {0.4,1}]&,3], Point[#]}& /@ pts},
MapIndexed[Text[#2[[1]], #1+{.005,0}]&, pts]}],
AspectRatio->Automatic,
DefaultFont -> {"Times-Bold", 12},
PlotRange->{{-.03,1.03}, {-0.03,1.03}}]

End[]; EndPackage[];

Code For Figures

Code for Figure 1

Needs["Graphics`Arrow`"];
Show[Graphics[{Arrow[{-1.4, -0}, {-0.7, -0}],
AbsoluteThickness[2],
Line[{{-3, 0}, {-2, 0}}], 
Line[{{0, 0}, {1, 0}, {1.5, -Sqrt[3]/2}, {2, 0}, 
   {3, 0}}]}], 
   PlotRange -> {{-3, 3.2}, {-1, 0.3}}, AspectRatio -> Automatic];

Code for Figure 4

Block[{$DisplayFunction = Identity},
 Do[im[d] = FractalTurtle[
 {"U" -> "+VF-UFU-FV+", "V" -> "-UF+VFV+FU-"}, "U", d, 90,
 Frame -> True, TurtleStep -> 2^-(d), 
 Axes->Automatic, PlotRange -> {{0, 1}, {0, 1}},
  StartPosition -> {2^-(d+1), 2^-(d+1)},
  FrameTicks -> {Range[0, 1, 0.5], Range[0, 1, 0.5], None, None}], {d, 6}]]
  
Show[GraphicsArray[{
  {im[1], im[2]}, {im[3], im[4]}, {im[5], im[6]}  } ]];

Code for Figure 5

Block[{$DisplayFunction=Identity}, Do[
FractalTurtle[{"U" -> "+VF-UFU-FV+", "V" -> "-UF+VFV+FU-"}, "U", 
 d, 90, Frame -> True, Turtle -> 2^-d, 
 Axes->Automatic, PlotRange -> {{0, 1}, {0, 1}},
StartPosition -> {2^-(d+1), 2^-(d+1)},
FrameTicks -> {Range[0, 1, 0.25], Range[0, 1, 0.25]}];
pathTemp[d] = TurtlePath, {d, 2}]]

Needs["Graphics`Arrow`"];

Show[Graphics[{
Line[pathTemp[1] /. {x_, y_} :> {x - 2.9, y}],
Line[pathTemp[2]],
Arrow[{-1.3, 0.8}, {-.4, .8}]

}], AspectRatio -> Automatic];

Code for Figure 7

f[s_]  :=  FontForm[s,  {"Times-Bold",  12}]\(im1\  = \ 
      Show[Graphics[{

            Line[{{0, 0}, \ {1, 0}, \ {1, 1}, \ {0, 1}, \ {0, 0}}], 

            Line[{{0.5, \ 0}, \ {0.5, \ 1}}], 

            Line[{{0, \ 0.5}, \ {1, \ 0.5}}], 

            Text[f@"0", \ {0.25, \ 0.15}], 

            Text[f@"0.3", \ {3/4, \ 1/4 -  .1}], 

            Text[f@"0.1", \ {1/4, \ 0.85}], 

            Text[f@"0.2", \ {3/4, \ 3/4 +  .1}], 
            
{AbsoluteDashing[{1, 2.5}], 
\ \ 
              Line[ll = {{1/4, \ 1/4}, {1/4, \ 3/4}, \ {3/4, \ 3/4}, 
                    \ {3/4, \ 1/4}}\ ]}, 
\ \ AbsolutePointSize[3], \ 
            Point\ /@\ ll

}], \ 
        Ticks\  -> \ Array[Range[0.25, \ 1, \ 0.25]&, \ 2], 

        DisplayFunction\  -> \ Identity, 
\ Axes -> Automatic, 
        AspectRatio\  -> \ Automatic]; 
\ 
d\  = \ 2; 

    FractalTurtle[

      \ {"U"\  -> \ "+VF-UFU-FV+", \ 
        "V"\  -> \ "-UF+VFV+FU-"}, \ "U", \ d, \ 90, 
\ 
      TurtleStep\  -> \ 2^\(-d\), \ \ 
      StartPosition\  -> \ {2^\(-\((d + 1)\)\), \ 2^\(-\((d + 1)\)\)}, 
\ 
      Axes\  -> \ Automatic, \ Frame\  -> \ True, \ 
\ 
      PlotRange\  -> \ {{0, \ 1}, \ {0, \ 1}}, 
      DisplayFunction\  -> \ Identity, 
\ 
      FrameTicks\  -> \ {Range[0, \ 1, \ 0.25], \ Range[0, \ 1, \ 0.25]}]; 

    
im2\  = \ \ 
      Show[Graphics[{
\ {AbsoluteDashing[{1, 2.5}], \ Line[TurtlePath]}, 

            Line[{{0, 0}, \ {1, 0}, \ {1, 1}, \ {0, 1}, \ {0, 0}}], 

            Line[{{0.5, \ 0}, \ {0.5, \ 1}}], 

            Line[{{0, \ 0.5}, \ {1, \ 0.5}}], 

            l\  = \ ToString\ /@\ 
                ToExpression[
                  StringTake[
\ \ 
                    First[ToString\ /@\ 
                        InputForm[BaseForm[Range[0, \ 1, \ 1. /16], \ 4]]], 

                    \ \ 111]]; \ 

            j\  = \ {{0, \ \(- .065\)}, {0, \ \(- .065\)}, {0, \ 
                  \(+ .045\)}, {0, \ \(- .065\)}, { .07, \ 0.013}, 
                
{0, \  .05}, 
                
{0, \  .05}, {0, \ \(- .065\)}, {0, \ \(- .065\)}, {0, 
                  \  .05}, {0, \  .05}, 
                
{\(-0.078\); \(-0.065\), \ 0.013}, {0, \ \(- .065\)}, {0, 
                  \  .065}, {0, \ \(- .065\)}, {0, \ \(- .065\)}\ \ }; 

            MapIndexed[

              Text[\ FontForm[
\ \ l[\([First@#2]\)], 
                    \ {"Times-Roman", \ 8}], 
                  \ \ #1\  + \ j[\([First@#2]\)]]&\ , \ 
\ \ \ TurtlePath], 
            
\ \ \ 
\ \ \ {AbsolutePointSize[3], \ Point\ /@\ TurtlePath}, 
            
\ \ Line[{{ .25, \ 0}, \ { .25, \ 1}}], 
\ \ 
            Line[{{ .75, \ 0}, \ { .75, \ 1}}], 
\ \ 
            Line[{{0,  .25}, \ {1,  .25}}], 
\ \ \ 
            Line[{{0,  .75}, \ {1,  .75}}]


}], \ 
        Ticks\  -> \ Array[Range[1/4, \ 1, \ 1/4]&, \ 2], 

\ 
        Axes -> Automatic, \ AspectRatio\  -> \ Automatic, 
\ 
        DisplayFunction -> Identity]; 
\ 

    Show[GraphicsArray[{{im1}, \ {im2}}]]; \)

Code for Figure 8 (* WARNING: You must have System 3.5 for this code to work, as a bug in 3.0 kills it in 3.0 *)

\(d = 6; 

    FractalTurtle[{"U" -> "+VF-UFU-FV+", 
        "V" -> "-UF+VFV+FU-"}, "U", d, 90, 
      TurtleStep -> 2\^\(-d\), 
      StartPosition -> {2\^\(-\((d + 1)\)\), 2\^\(-\((d + 1)\)\)}, 
      DisplayFunction -> Identity]; \)

\(\(Show[
      Graphics3D[{AbsoluteThickness[0.4], 
			
          Line[MapIndexed[
              Flatten[{\((#2 - 1)\)/Length[TurtlePath], \ #1}, \ 1]&, \ 

              \ \ \ \ TurtlePath]]}], Axes -> Automatic, 
      BoxRatios -> {1.5, 1, 1}, 
	
      AxesLabel -> 
        \((FontForm[#, {"Courier-Bold", 12}]&)\)/@{"t", "x", 
            "y"}]; \)\)

Code for Figure 10

image[d_?OddQ] := FractalTurtle["U" -> "iU+iU+Ui+++Ui++", 
 "+U", (d-1)/2, 45,
 TerminalSubstitution -> ("U" -> "F--F+"),
 TurtleStep -> 2^(-(d-1)/2-0.5),
 PathStyle -> AbsoluteThickness[2],
 Axes -> True,DefaultFont -> {"Courier", 8},
 Ticks -> {Range[0, 1, 1/4], Range[0, 1, 1/4]},
 PlotRange -> {{0, 1}, {0, 1}}]
 
image[d_?EvenQ] := FractalTurtle["U" -> "+iUi--U+UiUi-", 
 "U", d/2 - 1, 90,
 TerminalSubstitution -> ("U" -> "+F-FF-F"),
 TurtleStep -> 2^(-(d/2)), 
 PathStyle -> AbsoluteThickness[2],
 Axes->True,
 DefaultFont -> {"Courier", 8},
 Ticks -> {Range[0, 1, 1/4], Range[0, 1, 1/4]},
 PlotRange -> {{0, 1}, {0, 1}}]
 
Show[GraphicsArray[Block[{$DisplayFunction = Identity},
Map[image, {{1, 2, 3}, {4, 5, 6}}, {2}]], 
GraphicsSpacing -> {.15, .1} ]];

Code for Figure 11

\(\(
Clear[im]; 

    Do[im[i] = 

        FractalTurtle[{"F" -> "F+f---F++F", 
            "f" -> "+F+f---F+f"}, 
		"F++F++F++F", 
	i, 45, 
          
	StartPosition -> {2, 2}\^\(-\((i + 2)\)\), 
          TurtleStep -> 2\^\(-\((i + 1)\)\), 
	
          SecondaryStepSize -> 2\^\(-\((i + 0.5)\)\), 
          DisplayFunction -> Identity, 
			
          PlotStyle -> AbsolutePointSize[0.65], 
	
          PlotRange -> \(({{0, 1}, {0, 1}})\)], \ {i, 1, 3}]; 


    Show[GraphicsArray[{im[1], \ im[2], \ im[3]}]]; \)\)

by Stan Wagon

Readers familiar with space-filling curves may well be intrigued by the problem of using Mathematica to generate them. These objects are complicated, and it is not immediately clear what the best approach is. In this column I will show how to simulate a turtle and use it to generate approximations to a space-filling curve. As often happens, Mathematica’s superb 2- and 3-dimensional graphics allow us to generate some new ways of visualizing these monstrous objects. Of course, there is ample opportunity to learn some tricks to save time and memory.

1. A Literate Turtle

The turtle, an idealized traveler on a computer screen, will be familiar to readers who have experience with the LOGO language: it is an abstract object that wanders about in the plane. We think of it as occupying a point and able to accomplish certain simple tasks: a step forward or backward, or a right or left turn through a fixed angle. Because Mathematica generates entire images at once, we cannot easily show the real-time motion of the turtle. Instead, we will generate a complete sequence of moves, and then show the turtle's path as a sequence of line segments.

Our approach is based on the idea of string rewriting (also known as “L-systems”, after A. Lindenmayer, who used them to generate images that simulate plant growth; see [Prusinkiewicz and Lindenmeyer 1990] and chapter 7 of [Wagon 1998]). The idea is to attach motions to certain characters: "F" means forward, "B" means backward, and "+" and "-" denote a right and left turn, respectively. We can start with a simple string, such as "F", and then repeatedly apply a rule such as "F" -> "F+F--F+F". Such a substitution is simple to implement using StringReplace; we very quickly generate long strings.

Nest[StringReplace[#, "F" -> "F+F--F+F"] &, "F", 3]

The next step would be to have our turtle scan the characters of the string and perform the indicated moves. On the face of it, this seems simple. Break down the string into its characers and have the turtle, as it scans each character, either add a point to a list of points it has visited (if it sees an "F" or "B") or change its heading if it sees "+" or "-". But there are a couple of complications.

Recall that AppendTo is a natural way to add a point to a list. This command has the advantage that we need not know in advance how long the list will ultimately be. This is indeed convenient, but AppendTo turns out to be quite slow, so much so that it should never be used to build long lists. It is much better to build the list as a massively nested object and then flatten it appropriately. Here is an example. I leave it the reader to compare this technique to others on examples using thousands of points. I believe this is the fastest way to do it.

l = {}; Do[l = {l, {Random[], Random[]}}, {10}]; Partition[Flatten[l], 2]

With this preamble, the details are straightforward. We maintain the current direction in dir and the current path-end in last; path starts out as the empty set and grows appropriately. The key pieces of the code are the definition of turtle, which takes a string, looks at the characters, and carries out the proper instructions. In this example we assume that the turtle starts at the origin, turns through 60˚, and has steps having length 1.

\(turtle[s_String]\  := \ 
      \((θ = 60  °; 
\ \ \ dir = {1, 0}; 
\ \ \ 
        rotateleft
           = {{Cos[θ], \(-Sin[θ]\)}, {Sin[θ], 
              Cos[θ]}}; 
\ \ \ rotateright = Transpose[rotateleft]; 

        \ \ \ last = {0., 0.}; 
\ \ \ path = {last}; 
	
        Characters[s] /. {
		\ 
            "+" :> \((dir = rotateleft.dir; )\), 
		\ 
            "-" :> \((dir = rotateright.dir; )\), 
		\ 
            "F" :> \((path = {path, last += dir}; )\), 

            		\ "B" :> \((path = {path, last -= dir}; )\)}; 
        						
	path = Partition[Flatten[path], 2])\)\)

Now we start with the string "F++F++F", which instructs the turtle to trace out a counterclockwise equilateral triangle. When we then use the rewrite rule "F" -> "F-F++F-F", we are asking the turtle to replace each single step by a four-step walk, as illustrated in Figure 1. Of course, repeated application of this recursive idea leads to the Koch snowflake, perhaps the most famous of the family of curves known as fractals.

FIGURE 1:  A graphic view of the rewriting step that produces a Koch snowflake.

To continue, we define the starting strong and the rewrite rule, and then iterate to get longer strings to feed to the turtle.

rewrite = "F" -> "F-F++F-F"; \(start = "F++F++F"; 

    Show[Graphics[
        Line[turtle[Nest[StringReplace[#, \ rewrite]\ &, start, 1]]], 
		
        AspectRatio -> Automatic]]; 
\)

FIGURE 2:  A single application of the rewriting rule leads to instructions that cause the turtle to trace out a six-pointed star.

The electronic supplement contains a fleshier version of this code, called FractalTurtle. It takes four arguments: the rewriting rule(s), the starting string, the depth, and the turning angle in degrees. It has several options, including the ability to set the starting position and direction and the length of each step.

FractalTurtle["F" -> "F-F++F-F", "F++F++F", 3, 60];

FIGURE 3:  After three applications of a rewriting rule, the turtle traces out the familiar Koch snowflake, in this case a 192-sided polygon.

Note that after a use of FractalTurtle, the path, called TurtlePath in the package, is available for use in other graphics.

Short[TurtlePath, 2]

Of course, the reader might wish to play with variations obtained by changing some of the Fs to Bs, changing the type or number of the turns, and so on.

2. A Space-Filling Turtle

We will go into some of the mathematical details in section 3, but first let us consider the problem of how to get our turtle to trace out the shapes in Figure 4.

FIGURE 4:  Six approximations to a space-filling curve. These images can be generated by judicious use of two string rewriting rules.

Evidently, we need a recursive scheme that turns each U-shape, in any of four orientations, into a complicated shape made up of four linked U-shapes, as shown in Figure 5.

FIGURE 5:  A graphic view of the main recursive step to construct a space-filling curve.

The following rewrite rules do the job:

  "+F-F-F+" -> "+-F+F+F-F-+F-F-F+F+F-F-F+-F-F+F+F-+"
  "-F+F+F-" -> "-+F-F-F+F+-F+F+F-F-F+F+F-+F+F-F-F+-"

where the starting string is "+F-F-F+" . The code that follows generates the depth-3 approximation, with the TurtleStep and StartPosition options used to scale and center the path. Note that ordinary graphics options, such as FrameTicks, are just passed along.

d = 3;
FractalTurtle[
	\({"+F-F-F+"\  -> \ 
                  "+-F+F+F-F-+F-F-F+F+F-F-F+-F-F+F+F-+", 
\ \ \ 
                "-F+F+F-"\  -> \ 
                  "-+F-F-F+F+-F+F+F-F-F+F+F-+F+F-F-F+-"}\), "+F-F-F+", d, 90,TurtleStep -> (2^(-(d + 1))),StartPosition 
              -> {(2^(-(d + 2))), (2^(-(d + 2)))},Axes -> None,Frame -> True,PlotRange -> ({{0, 1}, {0, 1}}),FrameTicks 
              -> {Range[0, 1, 1/4], Range[0, 1, 1/4],  None, None}];

FIGURE 6:  An approximation to a space-filling curve can be traced out by a turtle following two rewriting rules.

Of course, there are many, many geometric shapes that can be generated by our turtle. We present just three examples. The reader might enjoy trying to guess what shapes will result from these rules!

FractalTurtle["F" -> "F-F++F-F", "F-F-F-F-F-F", 2, 
      60];

FractalTurtle["B" -> "+FB--FB+B", "FB", 5, 60];

In the next example, "Y" is just a placeholder. It plays a role in the rewriting process but, at the end, nothing happens at a "Y" because turtle["Y"] is undefined.

FractalTurtle[{"F" -> "FF", 
        "Y" -> "F-Y+Y+Y-F"}, "Y-F-F", 4, 120];

3. But Have We Really Filled Space?

So far we have not been too precise about what a space-filling curve is and what the shapes in Figure 4 have to do with a space-filling curve. The definition is simple: a space-filling curve is a continuous function that maps the unit interval onto the unit square. G. Peano, in 1890, was the first to construct such a beast. We present here an example due to David Hilbert. Here is how one can rigorously construct Hilbert’s space-filling curve H. Given a real number t in the unit interval, we can define H(t) by writing t in base 4. First divide the square into four equal squares and label them, with 0, 0.1, 0.2, 0.3 as in Figure 7(a). Then subdivide each of these squares and label them with 0, 0.01, 0.02, . . . , 0.33 as in Figure 7(b). Note how a path through the centers of the 16 labeled squares in order of the labels corresponds to our 4-U object. Continue subdividing and labeling in this way forever. Now, every real number t has a base-4 representation and considering initial segments leads to a nested sequence of squares. Since nested sequences of closed sets whose sizes approach zero converge to a unique point, we have defined the point that will be H(t). Thus, for example, 0’s representation consists of all 0s, and so its nested sequence converges to the point (0, 0).We will omit the details of the proof that H is continuous. It is fairly clear that H covers the entire square, since every point P in the square lies in the intersection of a unique sequence of squares. That sequence, via our convoluted labeling, defines the real t such that P = H(t).

FIGURE 7:  Tracing through squares in ever-finer grids according to precise rules leads to the formal definition of Hilbert’s space-filling curve.

So now we can see exactly how the piecewise linear paths of Figure 4 approximate the space-filling curve. Given a number t, use initial segments of its base-4 representation to locate points on the approximations. These points converge to H(t).

There are a few misleading aspects of this point of view, however. First, all the approximations are one-to-one. But the ideal function H is definitely not one-to-one: there is no continuous one-to-one function from an ineterval onto a square. And second, note that there is little point in buying more memory and generating deeper approximations to H; one would soon end up with an uninteresting square that is completely black. We will address these two issues in the next section.

4. The True Graph of a Space-Filler

One shortcoming of viewing the space-filling curve in the target square is that deeper images of the approximation are not visually enlightening. But remember that a function from an interval to the square has a graph that lives in 3-space. It is very simple for us to view an approximation of this graph: we simply generate the planar paths as before, and then attach the proper t-coordinate to each point. Set d to a smaller value if you have insufficient memory, because depth 6 leads to a path with (4^6) or, 4096 points.

\(d = 6; 

    FractalTurtle[{"U" -> "+VF-UFU-FV+", 
        "V" -> "-UF+VFV+FU-"}, "U", d, 90, 
      TurtleStep -> 2\^\(-d\), 
      StartPosition -> {2\^\(-\((d + 1)\)\), 2\^\(-\((d + 1)\)\)}, 
      DisplayFunction -> Identity]; \)

\(\(Show[
      Graphics3D[{AbsoluteThickness[0.4], 
			
          Line[MapIndexed[
              Flatten[{\((#2 - 1)\)/Length[TurtlePath], \ #1}, \ 1]&, \ 

              \ \ \ \ TurtlePath]]}], Axes -> Automatic, 
      AxesLabel -> {"t", "x", "y"}]; \)\)

FIGURE 8:  A view of the true 3-dimensional graph of Hilbert’s space-filling curve.

This view brings home the remarkable property of a space-filling curve: remember, it is not simply that it comes close to every point in the square; the ideal curve actually passes through every single point in the square.

An added bonus of having all the points in hand is that we can view x as a function of t. It turns out that, using the ideal curve H, this first projection function is a continuous function that fails to be differentiable at any point in its domain. We can see an approximation to such a continuous, nowhere differenetiable function by just plotting the first coordinates of the path.

a = 100; ListPlot[
      Transpose[{Range[0, 1, (1./(a - 1))], Take[First/@TurtlePath, a]}], 
      PlotJoined -> True, PlotStyle -> AbsoluteThickness[0.4]];

FIGURE 9:  Looking at only the x-coordinates of a space-filling curve leads to a continuous function that fails to be differentiable at every point of its domain.

5. Another Viewpoint

There is another way of looking at Hilbert’s curve. We could look at the exact values of, say,  H(0), H(1/2), and H(1) and connect them linearly. Then we could do the same with H(0), H(1/4), H(1/2), H(3/4), H(1); and so on. 

First we need to improve the turtle. There are many ways to enhance a turtle, some of which are discussed in chapter 5 of [Wagon, 1991]. One simple enhancement is to add a flip command (we use "i" for inversion). Upon seeing an "i" the turtle flips onto its back, thus switching left and right. FractalTurtle has this ability. One can then generate the Hilbert paths using only one substitution rule as follows.

FractalTurtle["U" -> "+iUiF-UFU-FiUi+", "U", 3, 90];

Another enhancement is the ability to generate long words using meaningless letters, and then, at the end, replace some of the letters with a string of letters corresponding to motions (i.e., a string such as "F++F++"). We call this a terminal substitution, and it is implemented as an option to FractalTurtle. We leave it to the reader to work out why the rules in the following code do the job; they do generate the approximations to Hilbert’s curve that are perfectly correct on the appropriate binary fractions as discussed above. Although they look quite different than the traditional approximations of Figure 4, a comparison of the 3-dimensional view derived from this approach with Figure 8 will confirm that thay approach the same limit.

image[d_?OddQ] := FractalTurtle["U" -> "iU+iU+Ui+++Ui++", 
 "+U", (d-1)/2, 45,
 TerminalSubstitution -> ("U" -> "F--F+"),
 TurtleStep -> 2^(-(d-1)/2-0.5),
 PathStyle -> AbsoluteThickness[2],
 Axes -> True,
 Ticks -> {Range[0, 1, 1/4], Range[0, 1, 1/4]},
 PlotRange -> {{0, 1}, {0, 1}}]
 
image[d_?EvenQ] := FractalTurtle["U" -> "+iUi--U+UiUi-", 
 "U", d/2 - 1, 90,
 TerminalSubstitution -> ("U" -> "+F-FF-F"),
 TurtleStep -> 2^(-(d/2)), 
 PathStyle -> AbsoluteThickness[2],
 Axes->True,
 Ticks -> {Range[0, 1, 1/4], Range[0, 1, 1/4]},
 PlotRange -> {{0, 1}, {0, 1}}]
 
Show[GraphicsArray[Block[{$DisplayFunction = Identity},
Map[image, {{1, 2, 3}, {4, 5, 6}}, {2}]], 
GraphicsSpacing -> {.15, .1} ]];

FIGURE 10:  These piecewise linear approximations to the Hilbert space-filling curve have two advantages over the traditional view in Figure 4: (1) they are not one-one (just as the ideal curve is not one-one), and (2) the points at all the bends are perfectly correct: the first curve connects H(0), H(1/2), and H(1), the second connects H(0), H(1/4), H(1/2), H(3/4), and H(1), and so on.

5. The Obsessive Traveling Salesman

If curves that fill space are not mind-blowing enough, here is something that is totally outrageous. These abstract objects, among the most complicated constructions of 19th century mathematics, have been used to guide truck drivers on their delivery routes in Atlanta! J. J. Bartholdi and his coworkers at Georgia Tech have realized that space-filling curves can provide an important and practical heuristic to the traveling salesman problem. The idea is simple: given a bunch of points in a square, check out the order in which an obsessive traveling salesman — one who visits every single point in the square — would visit the points. Of course, one simply uses an approximation to a space-filling curve, though the curve should be chosen to be a loop, starting and finishing at the same point.. Figure 11 shows approximations to such a loop. Generating these curves requires some advanced features of the Turtle package, since there are two step sizes in each image. Here is the code for these images, where "f" is used for a secondary step, whose length is set by the SecondaryStepSize option.

\(i\  = \ 2; 

    FractalTurtle[{"F" -> "F+f---F++F", \ 
        "f" -> "+F+f---F+f"}, 
	\ "F++F++F++F", 	i, 45, 

      	StartPosition -> {2, 2}\^\(-\((i + 2)\)\), 
      TurtleStep -> 2\^\(-\((i + 1)\)\), 
	
      SecondaryStepSize -> 2\^\(-\((i + 0.5)\)\), 
	
      PlotRange -> \(({{0, 1}, {0, 1}})\)]; \)

FIGURE 11:  Some approximations to a space-filling loop. Such loops can guide a traveling salesman through a region so that he visits points in an order that is close to optimal as far as total length is concerned.

For more details see [Bartholdi and Platzman 1988] or [Wagon 1998]; the latter contains a Mathematica implementation. In fact, their method has been used to deliver meals-on-wheels in Atlanta! A great advantage is that a computer is not needed; it is sufficient to have a space-filling curve superimposed on a city map! And of course it is very easy to add or remove addresses from the list. I am not a great believer in supporting all kinds of abstract mathematical research on the grounds that one never knows what will be useful in the future. That said, it is indeed remarkable that one of the most abstract constructions of the 19th century has found its way to a real-world problem 100 years later.

REFERENCES

Bartholdi, J. J., and L. K. Platzman, Heuristics based on space-filling curves for combinatorial problems in Euclidean space, Management Science 34 (1988) 291–305.
Prusinkiewicz, P., and A. Lindenmayer, The Algorithmic Beauty of Plants, Springer, New York, 1990.
Wagon, S., Mathematica in Action, Springer/TELOS, New York, 1988.

ABOUT THE AUTHOR

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition, and is especially appreciative of the new ways Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He also teaches a high-altitude course, Rocky Mountain Mathematica, every summer in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu