**Kleinian and Quasi-Fuchsian Limit Sets: An Open Source Toolbox
Genotype:
1.1.33 PDF Overview**

Chris King

**Download** the complete **Matlab toolbox, C scripts and Mac app** with web resource

This page contains an in-depth article with a Mac app, a cross-platform C implementation and a Matlab toolbox that depict limit sets of Kleinian and quasi-Fuchsian groups.

**Contents**

- Overview of the Computational Package
- Diverse Mathematical Foundations
- The Depth First Search Algorithm
- Rational Maps and Period Counting
- Finding the Elusive Traces
- The Maskit Slice as a Fractal Solution Space
- Endless Varieties among Limits of Limits
- Lost in Four Dimensional Space
- Finding Conjugate Limit Sets
- Rational Periodicity becomes Irrational Degeneracy
- The Total Eclipse of Double Degeneracy
- The Spiral Nemesis
- The Holy Grail at the Edge of Chaos
- Discreteness Enters the Chaotic Milieu
- The Non-Free Case and Automatic Group Tables
- The Stochastic Algorithm
- Escape-time and Extra Dimensions
- Program Listings

**Overview of the Computational Package**

The package provides a computational toolbox in an OS-independent form which is freely available to those seeking to understand these systems and explore them, including the key examples in "**Indra's Pearls: The Vision of Felix Klein**" by David Mumford, Caroline Series and David Wright, Cambridge Univ Pr. 2002 ISBN 978-0-521-35253-6. The images seek to fully explore and elucidate the underlying mathematical systems, but anyone seeking to produce fractal art via 3D or ray-traced versons could easily extend the toolbox algorithms in Matlab to do this.

The depth first search algorithm, derived from Indra's Pearls, is fully articulated in the** ****Matlab function**, which is also compiled in a fully functional native gcc **C script**, which is a 'smoking' hundred times faster. The C version is also encapsulated in the Mac OS app **DFS Viewer**. The Matlab toolbox includes supporting functions, has the capacity to find all the coefficients of high degree trace formulae and easily address new research questions. However all features can also be explored without Matlab, using DFS Viewer on Mac and the key native C scripts on any gcc C compatible platform. This makes the algorithm widely and freely available to the general public, as well as maths students and researchers.

The Mac app **DFS Viewer** has a variety of controls to aid successful viewing, most of which are also included as parameters in the **C script** and **Matlab** versions:

The **C script** has been configured to be generic and universal – able on any platform to generate very high resolution colour or B&W images and output them as .ppm files which are easily converted to png or jpg formats, so for all users, this is the most accessible and time efficient way to generate high quality images. It is 93 times faster on the right-hand example in fig 7 (15.2 secs vs 1422 secs for levmax = 1000 and epsilon = 0.0001). It can also be easily adapted to generate movies, as batch files, and has a timing meter, based on the successive tree branches of the first four search levels, which gives a good indication of the progress of a high precision computation.

A list of key fractional and degenerate traces explored in this article is also provided, which can be input into any of the versions. Also included are Matlab and C script versions of **tracepq** and **tracepqr** which can be used to find traces of new examples.

The as well as a variety of suplementary functions aiding further research, the Matlab toolbox also includes a stochastic algorithm which can depict chaotic cases, a colour filling escape-time algorithm, and a 3-D example.

**An Intersection of Diverse Mathematical Areas**

Kleinian and quasi-Fuchsian limit sets form an intersection of diverse branches of mathematics, becoming a marriage of group theory, the topology of handles and pinched manifolds formed by the transformations on the Riemann sphere, hyperbolic geometry, fractal tilings, and the theory of chaos and fractal dimensions

Kleinian and Quasi-Fuschsian limits sets are generated by a pair of Mobius transformations in the complex plane, *a *and *b*, along with their inverses *A = a *^{–1}* *and *B* = *b*^{ –1}. Mobius transformations are the analog of affine Euclidean transformations (translations, scalings, rotations) on the Riemann sphere (1 in fig2), the complex plane with a single point at infinity added, forming the north pole. Mobius maps are fractional linear transformations *m*(*z*) = (*pz*+*q*)/(*rz*+*s*), which compose by matrix left-multiplication:

Mobius transformations, as noted, correspond to complex versions of the rotations, scalings and translations of affine Euclidian maps.They can be decomposed into similar elementary transforms and are classified into elliptic, parabolic, hyperbolic and loxodromic types. Parabolic transforms are like transalations and have only one fixed point on the Riemann sphere, like the point at infinity of a translation. Loxodromic maps are complex dilations. They have a repelling and an attracting fixed point and the map expands in a spiral form from the repeller to the attractor. Elliptic maps are conjgugate to pure rotations. Pairs of Mobius transformation induce a group structure of transformations under composition. The limt sets consist of the asymptotic points the rest of the Riemann sphere is eventually mapped into under the group transformations. Conjugacy powerfully simplifies the classification, because any map *M *=* C*^{-1}*NC* generates isomorphic behavior, so we can regard the groups as equivalent and reduce the analysis to conjugacy classes, or convert our mappings into convenient conjugate forms. The interplay of the group's generators *a*, *b*, *A*, *B *become defining words in the construction of the limit set algorithm, were the exact action of each generator on the limit set becomes apparent, as in fig 1(right) and fig 4. The traces of parabolic maps, defined by *p* + *s* above (the sum of the diagonal elements) lie on the real line between –2 and 2. Those of elliptic maps lie on the unit circle excepting 1 and –1. Traces of hyperbolic maps lie on the real line outside [–2,2] and loxodromic maps lie anywhere in the complex plane outside these ranges, although hyperbolic are sometimes included in loxodromic.

A brief summary of some of the mathematical ideas underlying these limit sets is illustrated in fig 2. Mobius mappings preserve circles and lines in the complex plane Lines are simply circles going through the north pole at infinity. When two Mobius transformations are combined, the actions of the matrices and their inverses generate a symmetry group, just as the 3 rotations and three reflections of an equilateral triangle form a symmetry group. These groups can take varying forms. If we take two disjoint pairs of circles and map the inside of one onto the outside of the other on the Riemann sphere, we also generate a topological transformation of the surface as shown in (e), generating a 2-handled torus – the pretzel. When we apply these two transformations to the surface, successive images of the pairs of Schottky circles become mapped fractally as shown in (b) forming a fractal dust or Cantor set – the limit set of the group. As the circles come together and kiss or touch as in (c) the limit set can become a circle. If the commutator is parabolic wiith value -2 the double torus becomes pinched as in (f) grandma's recipe can be used to generate a quasi-Fuchsian group consisting of a fractal quasi-circle forming a topologically closed (Jordan) curve. These transformations generate what is called a quasi-Fuchsian limit set, which is Fuchsian if the limit set is a geometrical circle as in (c). The topology of the surface can undergo other pinching and puncturing transformations as in (g), resulting in Kleinian groups such as the Apollonian gasket (partial Schottky computation in (d) and full limit set fig 1 (left), which is a limiting case of neighbouring quasi-Fuchsian sets illustrated in fig 1 left, as the traces* *tend from 2.2 down to 2. The groups generate tiling symmetries on the hyperbolic plane as illustrated by Maurice Escher (i) and the symmetries become reflected in the forms of the circle packings generated in many of the limit sets where the second matrix trace is 2. The limit sets are also intimately connected with hyperbolic geometry as examplified by the modular group of symmetries illustrated in the Escher print (i), with topological pinching (h) also echoed in the fractal structures of the circular limit sets were *TB* = 2.

In fig 3 one can see the effects of changing *Ta* and *Tb* from the 2, 2 of the gasket to complex values. While *Ta* alters the inner region to form spiral fronds, as in fig 1 right, *Tb* transforms the outer circle into exterior fronds. When both have complex values, the interactive effects can become unexpected, as in the case of + / + (left), which is not simply a partial reflection of the centre figure and is also about ten times larger. Changing the signs of real and/or imaginary parts yeild four symmetry versions of each figure, consistent with the symmetries of the Maskit plane and the double roots for finding *Tab*.

**The Depth First Search Algorithm – **Program listings in colour**: **(a) C script (b)** **Matlab

The Depth First Search algorithm is the quintessance of depiction methods of quasi-Fuchsian and Kleinian limit sets and is the key instrument for studying the diverse key examples. The DFS algorithm does a full depth search of all the key examples in the book, including the non-free cases, which require an automatic group schedule to take account of degenerate words. It includes intelligent line linking using the attracting fixed points of the commutators and special words to take account of the Apollonian Gasket and other examples with whorls with narrow necks. The essential inputs are the traces (*p*+*s*) of the two generating Mobius transforms *a* = [*p q*;*r s*], sending *z* to (*pz*+*q*)/(*rz*+*s*), plus other parameters to control the accuracy of the plots, noted in the discussion.

The gallery of high resolution 1100 x 1100 examples demonstrate the quality of the algorithm, with examples whch test it to the limits of degeneracy, 'accidental' coincidences of double cusps, and situations in which the group is not free due to relations between the generators.

The DFS algorithm has an ingenious strategy designed to work with transformations with parabolic commutators that produce quasi-Fuchsian groups which can be scanned in a circular manner by exploring the tree of words defined by compositions of the transformations. Firstly "grandma's algorithm" is used to generate Mobius transforms with parabolic commutators (*aba*^{-1}*b*^{-1} etc) from the two traces, *Ta* and *Tb *that are unique up to complex conjugacy. In the case of the Apollonian gasket, both *Ta* and *Tb *are 2 so *a* and *b* are also parabolic. One can see how grandma's algorithm is defined by reading the subfunction grandma() in the program listing.

The algorithm then does a depth-first exploration of the tree of addresses of all possible iterations of the transformations, proceeding around the addresses in a manner that ensures we steadily traverse the maze of address space and in particular, the fractal topological quasi-circle of the limit set, by ignoring cancelling paths such as *aa*^{-1}, storing the matrices of composed transformations going down the depth search, so these don't have to be re-calculated whenwe back up part way to explore another vacant branch. The cancellations are made by always going down the tree corkscrew-wise, making a right-hand turn after each level descent, and then exiting the level when we reach the inverse of the transform we entered the level by after three left turns to explore the other tunnels. It makes a branch terminal check by applying the composed transformation to the fixed point of the appropriate commutator defined by the last transform e.g. *t*(1)**t*(2)*...**t*(*n*-1)*a operates on the attracting fixed point of *bA*^{}*B*^{}*a*, which effectively defines the infinite path composed of the *n* terms and the infinite repeat of the right-hand commutator. It then compares this "new" point with the last "old" point found the same way (equivalent to using the fixed point of the opposite commutator *BAba*) and if they are less than *ε* apart draws a line between. This approach neatly overcomes a critical flaw in random iteration algorithms, that visit some parts of the limit set exponential rarely and thus effectively never get drawn in a finite time. This means that it doesn't need to explore further down the tree once it has a satisfactory linkage of the scale we are depicting, avoiding exponential runaway. Nevertheless on a dual core intel Mac at a resolution of *ε* = 0.0001, each image takes around an hour to produce, although moderately lower resolution images will display in a matter of seconds or minutes.

Special words are introduced to take account of two anomalies, Firstly in the Apollonian gasket example to close the top and bottom circles and their fractal offspring. Secondlythis process is also invoked to avoid pinching off the whorls of limit sets with large appendages with narrow necks as in fig 4. To do this, in addition to the attracting fixed points of the four sets of commutators we also include those of the original four generators, which, as can be seen from the first image in fig 1 and the image on the right from Indra's pearls, are at the very tips of the whorls. Fig 4 illustrates the problem for the second example in fig 3, where two "new" and "old" candidates at the neck are *ε-*close, but the "tip" of the whorl is further away and joining them would cut the whorl at the neck. This means we can then decide to fill in a line only if the distances from our "new" and "old" points to the local whorl tip we call "mid" are both less than *ε*. In the case of the gasket we fill "new" to "mid" to "old" as the 'tip' is right in the centre of the gap, but for other sets just "new" to "old" provided both *ε*'s are small enough.

Further examples aided by special words to get vey high reolution images efficiently arise when dealing with the cusp examples below, where other elements are chosen to be parabolic, such as *a*^{15}*B* in the 1/15 example. We then proceed by making a table of all the parabolic elements for each generator, including the commutators and joing "new" and "old" points only if all the elements arising in the tree exploration sequence between these two are within *ε*. These examples are not included as they would each apply to only a single case and the algorithm is already giving high reolution images of these in real time.

The many non-free examples in figs 19 – 22 go further and have generators connected by an algebraic relation (*aba*^{-1}*b*^{-1})^{2} = I, which couples them and causes many of the words in the search tree to become null, so we have to set up an automaton generating the table of an automatic group to correctly portray the set. This requires a new array "state" which iteratively generates successive word states by recursively entering the lookup table FSA, forbiding entry to tunnels involving any of the null word combinations.

I have so included a function specialwords() in the Matlab implementation so that anyone who wishes to set up a special words directory for other more complex examples can do so easily, as the plotting routine is already included to process the special word extensions, so all that is needed to be done is to replace the array entries in specialwords() with the ones one has defined from the permutations of the parabolic element. As an example in fig 6, I have included IPs 8.16 designed to show the method in the case of *aaB*. This ends up with the look up structure shown in fig 6 left, while the image on the right illustrates two levels of resolution of the limit set algorithm. The words in the tree maze can take multiple forms due to equivalences and cancellations of inverses. For example the commutator *ABab* is equivalent to *BAba* because its single parabolic fixed point is the same as its inverse. This means that a word address like *aaBaab*[*ABab*] = *aaBaab*[*BAba*] = *aaBa*[*abAB*]. One also has to bear in mind these examples apply only to one example each and there is a significant time cost in checking a large number of extra cases for higher denominator words, and their many permutations all having to be below threshold drives the computation further down the exponentiating levels of the search tree. The example above needs some final tuning in its parameters, indicated by irregulatiries in the low resolution image, but otherwise, the Matlab implementation gives a faithful rendition of a full spread of examples of limit sets in the book.

**Rational Maps and Period Counting**

Group theorists take a special interest in the spectrum of properties of limit sets where one trace factor is 2, both because they provide a space of examples which can be analysed and because many other limit sets can be shown to be conjugate to these, so below are several key examples.

In fig 7 you can see examples of rational cusps 1/15, 2/19 and 7/43, with the central chain of circles numbered to highlight the relationship with the denominator. In the left figure the fixed point of the parabolic element *a*^{15}*B* is in green and circles 1 and 2 are thus mapped to themselves. B carries C1 to C16 and then *a*^{15} carries it back to C1. In a similar way, B carries C2 to C17 and then *a*^{15} carries it back to C2. This means that *a*^{15}*B* maps both the disks C0 and C1 to themselves, the fixed points of *a*^{15}*B* must thus be on both these circles meaning there is one fixed point (parabolic) on the tangent point. b is parabolic with fixed point -*i* and pushes points out from its fixed point along clockwise circular trajectories. a is loxodromic and pushes the disks in the numbered spiral chain from right to left across the shrinking and spiralling into the sink on the left.The chain of 15 circles connecting to the opposite large circle 16 map to one another and the mauve circles are in turn emerging from the action of loxodromic *a* in a spiral from the repeller and proceeding in the orange numbered spiral to the attractor. In the 7/43 case, you can also see the way the chain is running in cycles of seven, corresponding to the numerator, and cycles of 2 in the 2/19 case.

These mappings give a good idea of the discrete action of the underlying group, where the isometries map discs to discs so that the images of a given point under the group transformations form a discrete set and the images of the discs effectively form a hyperbolic tiling of the Riemann sphere. In non-discrete cases, as in fig 25 and fig 26, where the stochastic algorithm is simply displaying the orbits of one or more points on the limit set under random sequences of transformations, one can see that the images of such points will become chaotically mixed, so that the isometries of a given point are no longer discrete sets. Any discrete group introduces a tiling of the plane, or higher dimensional space e.g. with the annular tiles exemplified by the hyperbolic tilings in fig 2(2,3).

More accurate plots of the small circles can be gained by including the parabolic elements such as* a*^{15}*B* for 1/15 or *a*^{3}*Ba*^{3}*Ba*^{2}*Ba*^{3}*Ba*^{2}*B * for 5/13 as special words, along with all their permutations, linking "new" and "old" only if all of the steps in the chain between, involving these parabolic critical points are within *ε*.

To explore a given type of limit set we need to be able to calculate thre traces which will have particular properties and periodicities in terms of the generators. We can work out the correct trace of parabolic words recursively from *T*(*p*/*q*) of lower denominator fractions, as they combine in the manner of mediants and Farey tree.

The traces combine in two key relations: (1) The "grandfather" identity *Tmn* + *Tm ^{-1}n* =

Key to this is the nature of any composition of the generators which is parabolic. As an exercise in finding the trace coefficients of some of these examples lets look at the cases 1/*n* where the parabolic element is *a*^{n}*B*. If we are not trying to use the special words algorithm at this point, we only need to find the correct value of *Ta.* Using the extended Markov and grandfather identities, we can immediatley see

*–2 = TabAB* = (*Ta*)^{2} + (*TB*)^{2 }+ (*TaB*)^{2}–* Ta.TB.TaB* – 2, or (*Ta*)^{2} + 4^{ }+ (*TaB*)^{2}–* Ta.TB.TaB *= 0, from which we can deduce that *TaB *= *Ta* + 2*i* . Then using the grandfather identity on *a* and *B *, with *TB* = 2, we get a series of equations *Ta *^{n }*B* = *Ta*.*Ta*^{ n-1}* B* –* Ta*^{ n -2}* B*, making it possible to write a recursive process to define *Ta* as an *n*th degree polynomial. We can do this readily in Matlab and at the same time solve the polynomial numerically by a Newton-type method as follows:

function s=recurscoeffs(stps,opt)

syms ta

a=cell(1,stps+1);

a{1}=2; % This is
TB

a{2}=ta+2i; %This is TaB

for i=3:stps+1

a{i}=expand(ta*a{i-1}-a{i-2}); % Ta^(i)B = Ta.Ta^(i-1)B-Ta^(i-2)B

end

f=symfun(a{stps+1},ta); % Make it a function

if opt==0

s=eval(vpasolve(f==2,ta)); %Solve to get the roots of f=2

else

s=eval(vpasolve(f==-2,ta));

end

This will give a list of *n *roots. Theory confirms that among all possible solutions to the trace equation (for a fixed p/q), exactly one will be discrete, so only one of these will have a credible root one can try. You can directly input the values such as s(7) into depthfirst as parameter, having run this function, giving full precision. Example values are 1/7 1.846276 – 0.09628*i*, 1/15 1.958591 – 0.01128*i*, 1/23 1.98179 – 0.00320*i*, 1/31 1.98987 – 0.00131*i*, and 1/60 1.99726820 – 0.00018236, tending towards *Ta = *2, (or *μ* = 2*i* ) in fig 11 centre.

Going beyond the 1/*n* fractions to *p*/*q*, we can also calculate the trace *Ta* for general fractions, because parabolic expressions like *a*^{3}*Ba*^{3}*Ba*^{2}*Ba*^{3}*Ba*^{2}*B* can also be generated recursively e.g. from *a*^{3}*B* and *a*^{2}*B* by similar relations. To be more precise, examination shows that the word *w*_{(p+r)/(q+s)} = *w*_{p/r}.*w*_{q/s}, in the manner of mediants and the Farey tree, so we can use the grandfather identity as follows:

*T**w*_{(p+r)/(q+s)} = *Tw*_{p/r}.T*w*_{q/s} – *T*(*w*_{p/r}^{-1}.*w*_{q/s}). These relations simplify in rather beautiful ways if we have generated a list of Farey tree predecessors. For example

*Tw*_{(3+5)/(8+13)} = *T**w*_{8/21} = *Tw*_{3/8}.T*w*_{5/13} – *T*(*w*_{3/8}^{-1}.*w*_{5/13}) = *Tw*_{3/8}.T*w*_{5/13} – *T**w*_{2/5}, because *w*_{3/8}^{-1}.*w*_{5/13} = (*a ^{3}Ba^{3}Ba^{2}B*)

The Fibonacci fractions used to illustrate this above are are all linked as immediate neighbours in the Farey tree (right). Now it turns out that any other pair of fractions that are Farey neighbours also have generator word sequences which will correctly cancel in the same way as above. Thus to find the correct trace polynomial for an arbitrary fraction, we simply need to determine the path down the Farey tree from the precursors to get to the fraction. So the function **tracepq()** firstly calculates the sequence of left and right handed moves to get from 0/1 and 1/0, (which represent *Ta* and *TB* = 2 respectively, with 1/1 corresponding to *TaB* = *Ta*+2*i *) to the particular fraction whose trace we are searching for, using the subfunction **fareyseq()**. It then builds up the trace of *p*/*q *iteratively using the grandfather identity. I previously had a much more cumbersome algorithm recursecoeffs() using recursion on arbitrary splits, which is still listed in the toolbox. The trouble with the top-down recursion is that you have to find a split which avoids degenerate fractions and is in Farey neighbours. Consequently some larger fractions failed to solve out correctly and the recursion is liable to exponential runaway.

A list of some fractional traces.

**Coefficient recursion** tracepq() program listing in colour

A sample listing of tracepq() for 2/7 iis shown below. The 7 roots are listed in increasing order of their real part, so that one can easily search the lowest few roots for the one which portrays correctly using depthfirst().

>> s=tracepq(2,7,0,1);

2 -1.53746746-0.21870978i

5 -1.40222578-0.58650847i

4 -0.39303611-0.25387457i

3 -0.38479895-1.27170824i

1 0.79526188-1.15961696i

6 1.25146576-0.06152298i

7 1.67080064-0.44805900i

We can then follow with a portrait using depthfirst()

depthfirst(s(7),2 ,0.0001,1000,1,3,0,1100);

The parameters in depthfirst select s(7) the 7th solution listed, 2 for *Tb,* followed by 0.001 accuracy scale 1000 maximum tree depth, 1/3 image scaling, sav = 3 both prints progress (≥2) and saves the image (mod(2)=1), option 0 is for a standard group without special words for the generators (see listing), and 1100 is the image size. Printing progress means you can see a series of stages of exploring the tree so you know progress is taking place for a detailed plot that may take an hour to complete. The report lists the time and olace of each time it falls back below level 4, and if the process reaches levmax it prints a dot. This can happen either because the limit set has issues or because it is larger than the scaling. Ideally one wants levmax large enough so that a nice curve on screen doesn't reach it for a given epsilon. One can also use sav=-1 as a test mode which will report whether the process reached levmax by setting outfl false. This will happen with a bad limit set and cause the process to strand, so this mode can be used within a function to test for a good limit set. The function scantraces() scans all the trace solutions *i*/*n* for a given *n* and automatically finds the one correct root which will produce a low resolution image and declare outfl=true to provide all the *i*/*n* solutions automatically. These can then be plotted in sequence using batchlimsets().

**The Maskit Slice as a Fractal Solution Space**

In figs 9 and fig10 below are shown the Maskit slice, consisting of all solutions to these polynomial equations for the trace of the commutator being both 2 and -2.

A group in which no simplifications other than cancelling out adjacent transformations and their inverses (*aA* and *bB*) is called free. Groups for which large products stay away from the identity I are called discrete. In the non-discrete, case *m*(*z*) can eventually come very close to *z*, resulting in a chaotic recurrent orbit which prevents the areas outside the limit set being consistently tiled by the transformations, resulting in the Cantor dust limit sets we can see using the stochastic algorithm. We are here looking for the double-cusp groups which hover on the boundary between discrete and non-discrete. The points in the white upper region all represent discretegroups, while 'most' points in the shaded region are not discrete. The upper (and lower) cuspy curve is called the Maskit boundary. *μ*-values on the boundary give groups right on the transition between order and chaos.

Above are shown all the fractional rotation states modulo 29, the Maskit slice up to *q* ≤ 51 generated using the script slicec, which collects all the roots found by tracepq(), saving them in mus.mat The conjugate roots constituting parabolic solutions can then be plotted using the function plotmus() . The "Maskit slice" is the plane of *μ*(*p*/*q*) highlighting all polynomial roots where the generating *p*/*q* word, such as a^{15}B for *μ*(1/15), is parabolic. In the Maskit parametrization, whose limit set corresponding to *Ta* = 1.91 + 0.1*i*, *Tb* = 2 is illustrated in fig 23, we have one parameter *μ* defining a Mobius transformation which is paired with translation by 2, which is parabolic with fixed point ∞. The two generators arise immediately from the parameters: *a*(*z*) = *μ* + 1/*z* and *b*(*z*) = *z* + 2. In the standard parametrization we get the same limit set in conjugate form from *Ta* = –*i**μ*(*p*/*q*), *Tb* = 2. We will thus refer to such *Ta* simply as *μ*(*p*/*q*)).

Theory shows that no value of *μ* with imaginary part between –1 and +1 could give a discrete group, while values with imaginary part more than +2 or less than −2 would always be discrete. The points in the white upper region all represent discrete groups, while 'most' points in the shaded region are not discrete. Some isolated *μ*-values in the shaded central region give discrete groups, but such groups are never free. For example *Ta *= 1 and *Ta* = 2^{1/2} , which emerge as solutions to *μ*(1/6) and *μ*(1/8) give non-free groups with *a *^{6} = I and *a* ^{8} = I respectively and isolated limit sets, fig 24, surrounded by non-discrete values.

**Endless Varieties of Limits of Limits**

In fig 11 are shown a sequence of limit sets of 1/*n* tending in limit to 1/i which should theortetically be the Apollonian gasket. Note that in the right superimposed figure for 1/31 in fig 11, the gasket coincides with the upper and lower circles and their fractal subcircles center. Note as well, the dark circles which have emerged in both 1/23 and 1/37, and are present as cicles interrupted by the central chain of circles in 1/15 (fig 8 left), forming a chain between the fractal loxodromic sources and sinks that lies in the limit along the large Apollonian circles. The limiting limit set results in a four-fold rotational symmetry, rather than the two-fold symmetry of the gasket.

If we examine the matrices of the gasket arising from grandma's recipe we have:

Hence we can see that (*a ^{n}B*)

One can take the limit of 2/(2*n*+1) and derive yet another distinct limit limit set which also has fourfold symmetry but has pairs of middle-sized circles on each diagonal rather than a single one, as shown left for *μ*(2/121). And we could do this for 3/(3*n*+1) and 3/(3*n*+2) and so on so we have an infinite number of different limits of limit sets, all tendi*ng to μ*(0/1).

**Lost in Four Dimensional Space**

When we move away from *Tb *= 2, and search for *Ta*, we no longer have a 2D Maskit slice of the parabolic generator complementing Tb = 2, but essentially a 4D space, in which we can't find all polynomial roots as before because the grandfather identities no longer produce polynomials. Fig 13a and discussion shows a way around this problem.

The two sets in fig 13a were generated by finding *Ta for Tb* = 2, using tracepq() for 1/7 and then letting *Tb* = *Ta *and using a tacepqs() for *Tb* ≠ 2 to get a new parabolic *Ta* for 1/5. The generalized equation is no longer a polynomial, so the root solving will give only one root, but you can enter a suggestion as the last parameter in vpasolve() to prompt it to search near the original value of *Ta*. As above, you can see hidden aspects of the periodicities in the limit set on the left in the five and seven crossing circles, which was the actual solution generated, and to a lesser extent in the congugated version right.

Chapter 9 of Indra's Pearls shows that limit sets with parabolic traces can be shown to be conjugate to a sister limit set with a parabolic *Ta* and *Tb* = 2, by a transformation of generators to a new pair (*u*, *v*) expressed as words in (*a*, *b*) so that *u*, *v* are both parabolic and (*a*, *b*) can also be expressed as words in (*u*,*v*), causing the new group to be a conjugate of the original. For example if *u* = *w*_{1/2} = *aaB* and *v* = *w*_{2/3} = *aaBaB* then
we have:

*aaB*(*aaBaB*)^{-1}*aaB* = *aaB*(*bAbAA*)*aaB* = *aaBbAbAAaaB* = *a*

(*aaB*)^{-1}*aaBaB*(*aaB*)^{-1}(*aaB*)^{-1}*aaBaB* = (*bAA*)(*aaBaB*)(*bAA*)(*bAA*)(*aaBaB*) = *B*.

If *a* and *B* are parabolic then so are *u* and *v* and finding their traces, which we can calculate by multiplying out the products of *a* and *B* and entering the traces into grandma's algorithm we can get matrices *u'* and *v'* with the same traces as the words* u* and *v*, with the commutator parabolic. Thus given any double cusp limit set defined by *Ta* and *Tb* = 2, we can find a double cusp limit set with *u'* and *v'* words in *a* and *B.* We can explore this numerically with the above example. If we apply tracepq(8,13), we get *Ta* = 1.62217905-1.28170820*i. *If we then put this *Ta *into grandma() with *Tb* = 2, we get two matrices *u' *and *v'. *If we now calculate the following traces for *w*_{1/2} and *w*_{2/3}, which are close to the middle of the Farey journey from 1/0 to 8/13, directly with Matlab, we find *TaaB = Tu'u'V'* = 1.55210536 – 0.91396228*i* and *TaaBaB *= *Tu'u'V'**u'V'* = 1.55210536 + 0.91396228*i,* conjugates to 8 decimal places! We can in turn plug this into either depthfirst algorithm and we get the results shown below in fig 13b.

We can see this process in action in the context of searching for doubly degenerate limit sets in the conjugacy between fig 15 (left) and fig 16 (centre), where the construction is key to undestanding how we can in fact get double degeneracy. That said, demonstrating the images in fig 13 are conjugate to ones with *Tb* = 2 is another matter. All we have proved here is a special case. In fact, if we took the alternative grandma's root to the original we would get a similar limit set to the one centre.

**Rational Periodicity becomes Irrational Degeneracy**

In fig 15 is an illustration of the singly-degenerate limit set defined by Troels' point, named after Troels Jorgensen –*μ*(*γ*)=*μ*((–1 + 5^{1/2})/2), of the golden mean limit of the Fibonnacci fractions representing the lowest point on the Maskit slice. The only noticeable difference between the right hand image in fig 14 and the limit set of Troels' point in fig 15 is right on the boundary of the horizontal circles, where *μ*(89/144) still has infinitessimal hammer heads defining the horizontal circles.

The Troels' point example above is deemed singly degenerate because half the region (in this case inside the bounding circle minus its circular replicates, has been engulfed by the limit set in a space-filling fractal of fractal dimension 2. Such points exist because the double cusp values on the boundary of the Maskit slice are countable and the (continuous) boundary has been proven to be uncountable so irrational vaues clearly exist. In fact, almost all values are irrational so degeneracy is the rule! Notice our manner of approach to the limit alternately moving left and right down the Farey tree 2/3, 3/5, 5/8, 8/13 etc. alternating on either side of the golden mean and tending to it in limit. Looking at the sequence in fig 14, one can see why this limit set is degenerate and englufs the free interior. Each successive Fibonacci fraction runs to a smaller pair of horizontal circles, which in limit engulf each of the horizontal circles and their fractal replicates. The fractally self-similar nature of the limit sets of Mobius maps of fractional linear transformations, means that the scale approaches a geometric space-filling limit.

From the fundamental transformations, one can that under L or R moves, the traces transform as follows:

(*Ta, TaB,* *Ta.TaB* – *TB*) <*– *L* –* (*Ta*, *TB*, *TaB*)* – *R*–> * (*TaB*,* TB,* *TBTaB* – *Ta*)

and these compose, so that for RL we get (*Ta*, *TB*, *TaB*)* – *RL* –> * (*Ta.TaB* – *TB**, TaB,** T ^{ 2}aB*.

The examples we have used of fractional traces *μ*(*p*/*q*) have been based on finding the roots the *q*-th degree polynomial using the polynomial solver vpasolve() in Matlab. This provides a complete list of all possible roots without prompting, but tends to break down as the degrees pass about 400, but a carefully crafted Newton iteration as used in Indra's Pearls, such as the one included as **tracepqr()**, which is provided in both Matlab and a native C version, so it can be used without Matlab on a Mac. When suitably tweaked, and given a suitable priming value, this can capture much higher fractions using an initial trace value such as a Farey-related lower fractional trace. One can input successive Fibonacci fractions as inputs to the next and easily reach *μ*(987/1597), *Ta* = 1.616888739503 -1.294321525461*i*, very close to the limiting value above, which was also generated this way in Indra's Pearls. The spirally degenerate examples in fig 17 have had their traces discovered by the Newton method.

**The Total Eclipse of Double Degeneracy**

One can take the degeneracy issue up a notch to doubly (or totally) degenerate if we can find a pair of traces engulfing both the complements e.g the outside and the inside sets. The trace parameters for a close approximation can be been derived using generator transformations, as described above.

We start with a high order Fibonacci fraction close to the Troels' point golden mean, generated from 8 left-right passes of Farey mediant addition from 0/1 and 1/0 down the Fibonacci fractions, for example using the function fareyadd() below, use trace recursion, to get the limit set of –*iμ*(987/1597) and apply a generator transformation to (*w*_{21/34}, *w*_{13/21}) the middle pair of terms in the left-right sweep, using grandma's algorithm as noted above, we end up with the two nearly conjugate traces: T (*w*_{21/34}) = 1.501347474086 – 0.865385203320*i* , Tr (*w*_{13/21}) = 1.501347474169 + 0.865385203218*i*. However this limit set (fig16 centre) is still singly degenerate, and has to be, since it is conjugate to the Troels' point limit set (fig 15). In actual fact the set looks closer to being conjugate to one of the right hand sets in fig 15, (which it is, since it is conjugate to 987/1597).

function [p0,q0,p1,q1]=fareyadd(s)

% e.g.[p0,q0,p1,q1]=fareyadd([0 1 0 1])

p0=0; q0=1; p1=1; q1=0;

for i=1:length(s)

if s(i)==1

p0=p0+p1;

q0=q0+q1;

else

p1=p1+p0;

q1=q1+q0;

end

end

To actually get a doubly degenerate group we aim to find a seqence of generator fractions *p _{n }*/

This limit set has intriguing properties linked to hyperbolic spaces and knot theory. In addition to the symmetries of *a*, *b* , there is an additional induced parabolic symmetry *c* , which arises because the LR move which transforms (*a*, *B*) into (*a*^{2}*B, aB*) must itself arise from a Mobius map which conjugates the limit set to itself, has the same fixed point as the parabolic commutator, and commutes with it, although performing a distinct 'translation' flipping the sets bounding the jagged line connecting the two principal centres (IPs p 339, 388). Robert Riley has shown that this system of transformations results in a gluing of 3D hyperbolic space via the discrete tiling to become the topological 3-sphere with the figure 8 knot (right) removed (doi:10.1016/j.exmath.2013.01.003). You can discover more of this knotty area by searching for "Not Knot" on YouTube.

Just as we did alternating left and right passes of successive fractions to get the Troels' point trace value, we can likewise look for sequences of moves which induce a spiral degeneracy. Starting at the pair of fractions 0/1, 1/0, we take ten left steps by means of Farey addition to arrive at 0/1, 1/10. Then a right step leads us to the pair 1/11, 1/10. Applying the same 10+1 steps again leads to 12/131, 11/120 and then to 143/1561, 131/1430. A Matlab and C script version of **fareyadd** are each provided to turn LR sequences into such pairs of associated fractions. Continuing this way several more times, we end up at –*iμ* = 1.9264340533 – 0.0273817919*i*. You can locate the sequence of fractional cusps on the boundary of the Maskit slice as at left and verify that they ascend near the highest cusp and then successive cusps rotate in exploded blowups of the boundary verifying their spiralling limit sequence. We can see why the engulfing takes place in terms of the sequence of Farey fractions, each stage of which leads to a smaller circle, which is broken by one further stage along the sequence seeking a smaller circle by breaking symmetry.This symmetry-breaking results directly from each higher fraction representing a deviation from the previous one so that the previous level horns do not unite but form a new engufling set of spirals.Two stages in construction of the spirally degenerate limit set are illustrated centre in fig 17. At each stage the limit set engulfs a step furtherinto the two large horizontal circles, with a spiral split again just beginning in the right hand image, confirming the process would continue in the same manner as for the golden mean Troels' point. Just as the LR pattern leading to the golden ratio is a limit of continued fractions approaching a quadratic irrational , so the value of the L^{10} R spiral degeneracy approaches a quadratic irrational giving *μ*(0.0916079) as the limit value of *Ta*:

**The Holy Grail at the Edge of Chaos**

We can also use the above techniques to define doubly spirally-degenerate limit sets based on the same development of double degeneracy applied to the spirally-degenerate example above. To see how this works, let's consider instead of alternating LRLR..., as in the doublydegenerate case above, the sequence LLRLLR..., which is a step to the L^{10} R in the above spiral case. The sequence LLR takes us as follows (*a,B*) > L > (*a*,*aB*) > L > (*a*,*aaB*) > R > (*aaaB*,*aaB*), so we need *TaaaB* = *Ta.TaaB* −* TaB* = *Ta, **TaaB* = *Ta.TaB* −* TB* = *TB *and *T*^{ 2}*a* + *T ^{ 2}B* +

Going another step further with L^{3}R... we have *Ta* = *TaTB* – *TaB*, *TB* = *TaTa *^{2}*B* – *TaB*, *Ta*^{ 2}*B* = *TaTaB* – *TB* plus the Markov identity, leading to *TB* = (*T *^{3}*a* – *Ta*)/(*T*^{ 3}*a* – 2*Ta* – 1), and *T*^{ 8}*a* - 4*T*^{ 6}*a* - 3*T*^{ 5}*a* + 4*T*^{ 4}*a* + 7*T*^{ 3}*a* + 3*T*^{ 2}*a* = 0, giving *Ta *= 1.727154868638 – 0.228442123457*i*, *TB * = 1.572495224561 + 1.575566700023*i*. The limit set fig 7a right shows the emergence of spiral degeneracy with both an external and an internal set of spiral fronds.

For L^{4}R ^{∞}, we have *TB* = (*T *^{4}*a* – 2*T ^{2}a*) / (

One can see from the repeating grandfather identity how this process extends readily to the doubly spirally-degenerate solution of L^{10}R^{ ∞}. The corresponding system of equations iteratively solves *Ta ^{n}B* back to

We can solve out the higher terms unsing the grandfather identity in terms of *Ta*, *TB* and* TaB* as follows:

syms a B aB

a2B=a*aB-B;

a3B=a*a2B-aB;

…

a10B=a*a9B-a8B;

a11B=a*a10B-a9B;

a11B=expand(symfun(a11B,[a B aB]));

a10B=expand(symfun(a10B,[a B aB]));

a11B(a, B, aB) =

a10B(a, B, aB) =

We can then assert *Ta* = *Ta*^{11}*B* and *TB* = *Ta*^{10}*B* , by setting the left-hand sides to *a* and *B* respectively. Since these equations are linear in *B *and *aB*, we can make *aB *the subject of both and cross multiply the two fractions to give the form below in which *TB *can again be made the subject of a fractional equation and we can then assert *Ta* = *Ta*.*TB* – *Tab, *or *TaB *= *Ta*(*TB –* 1) to get TaB as well.

We can again use Matlab to advantage to expand the cross product of the fraction above, to get* TB* as a fraction of powers of *Ta*:

Hence:

We can then substitute these back into the Markov identity and solve. In fact Matlab can solve the entire system outright for any L^{n }R, as shown right in the solvable Maskit 'slices' for
*Ta*, *TB*,
*n *≤ 53, coloured from red through to blue, with again the roots running along the top curve of *Ta *being the actual trace solutions. Notice how differently they behave!

syms a B aB

s=[0 0 0 0 0 0 0 0 0 0 1];

t0=a;

t1=aB;

t2=B;

for i=1:length(s)

if s(i)==0

h=expand(t1*t0-t2);

t2=t1;

t1=h;

else

h=expand(t1*t2-t0);

t0=t1;

t1=h;

end

end

aB1=solve(t0-a,aB);

aB2=solve(t2-B,aB);

B=solve(aB1-aB2,B);

aB=a*(B-1);

f=a^2+B^2+aB^2-a*B*aB;

g=expand(f);

h=simplify(g);

e=eval(vpasolve(g,a));

Setting this to zero and evaluating gives the somewhat surprising solution with large real values: *Ta* = 1.936872603726 + 0.021074052017*i*, Tb = 1.878060670271 – 1.955723310188*i*, but this is clearly the correct set of trace values.

We note however that the Matlab solution breaks down for 'mixed' sequences like LRL^{2} R because we end up with mixed powers of *a*, *B* and *aB*.

**Discreteness Enters the Chaotic Milieu**

When the DFS algorithm faces a limit set which is non-discrete and fails to tile the plane in an ordered fractal the search tree can fail to find closely connected points and so will appear to hang furnning to the end of levelmax each time before going onto the next branch. Many of these limit sets also have interesting properties, and can be explored using two different approaches. The DFS Viewer app has a timer function and an escape routine both of which can interrupt the process if too much time has elapsed or the iteration has become unbounded. This can throw up a surprising diversity of images, as shown in fig 17b where three slight variations on the example of fig 13 left are shown . However, because the strategy of the DFS algorithm is to draw lines beween closely identifiec points in the limit set, its depictions give unrealistic portraits of these examples because a large number of lines are drawn higgldy-piggldy only the end points of which are pitential membes of the limit set.

An alternative method to portray these processes is the stochastic algorithm below, which is simply remapping a set of limit points already found on the limit set. This has another drawback, that some regions of the limit set are visited exponentially less often, but the examples in figs 24 – 26 show a variety of non-discrete limit sets with chaotic iterative behavior, where the orbits are intermingled.

**The Non-Free Case and Automatic Group Tables**

The final examples in this section are non-free groups, in which the generators *a*, *b* are chosen so that *aba*^{-1}*b*^{-1}*aba*^{-1}*b*^{-1 }= I and hence these relations between the generators become trivial, causing the search tree to have new dead ends, which have to be eliminated by generating what is called an automatic group. The generators also no longer have parabolic commutators so we have to find the appropriate value of *Tab* from *Ta*, *Tb* by using the four alarm version of "grandma's algorithm" by solving for the trace *TC *= *T*(*aba*^{-1}*b*^{-1}) = 0 arising from the above relation, for the given values of *Ta*, *Tb*. If you examine line of the code with the comment %Grandma's three parameter double alarm
*tC* = 0 gives a quadratic in *Tab* and you can try out each solution. There is a short function below to do this. You can likewise find the p/q cusps for the non-free case. The commutator trace has moved from -2 to 0 so the only change needed to recurscoeffs() to calculate these is to the initial solution *TaB* = *Ta* + sqrt(2)*i, *giving the centre figure above for *μ*(1/10)*.* The central circles now no longer touch the outer ones, emphasized by their colouring in yellow, or grey in the coloured images.

We then have to set up an automaton(right) in the form of a table generating the dead ends in every part of the maze by recursively looking up successive words and terminating at neutralizing words such as *b*^{-1}*a*^{-1}*ba* and *aba*^{-1}*b*^{-1}*a*. You can see this right and encoded as FSA , where the 19 rows are given by I, *a*, *b*, *A=a*^{-1}, *B=b*^{-1}, *ab*, *AB, ba, bA, Ba, BA, abA, ABa, baB, bAB, Bab, BAb, abAB*, and *ABab*, with x counting as a death state. The table is used sequentially to take us from any state on the left by right multiplication by one of the four generators. The states are stored sequentially as we move down a level of the maze so they can be looked up and the next recalculated. The end result, centre, has two disconnected sets of circles, coloured white and yellow in the figure. It is also a double cusp in which *b* and *a*^{10}*b* are parabolic. Two further examples from Indra's pearls are shown left and right, each coloured by the first iterate, and three more below to illustrate the diversity, with the disconnected central region in yellow.

Here is a short function nonfree() to extract an appropriate *Tab* for the non-free case. Option selects which root (usually 0):

function tab=nonfree(A,B,opt)

%ta=1.924781-0.047529i;

%tb=2;

ta=A;

tb=B;

p=-ta*tb;

q=ta^2+tb^2-2;

if opt

tab=(-p+sqrt(p^2-4*q))/2;

else

tab=(-p-sqrt(p^2-4*q))/2;

end

You can use tracepq() in the non-free case to generate any fractional twist with *Tb* = 2, followed by non-free on the result to get the required *Tab *and then depthfirst to get the limit set. Below is a session for 3/7. The option is now 2 for non-free() and prt is 1 to print the identities. Non-free has option zero to pick the first root. The parameters in depthfirst() choose s(2) and now have option 2 for non-free and choose t for *Tab*.

>> s=recursecoeffs(3,7,2,1);

T(2/5)=T(1/2)*T(1/3)-T(0/1) 2.50 2.00 3.00

T(3/7)=T(1/2)*T(2/5)-T(1/3) 2.33 2.00 2.50

1
-1.61195573 -0.46319322

2
1.82867804 -0.60703463

3
-1.65805636 -0.71965515

4
-0.42695172 -0.16466462

5
-0.25623249 -1.15299930

6
1.03890744 -0.78635572

7
1.08561081 -0.34873804

>> t=nonfree(s(2),2,0);

>> depthfirst(s(2),2,0.001,1000,1/3,3,2,1100,t);

The limit sets show the differences in the non-free case where the commuator squared is the identity. The image on the right is the analogue of the right hand image of fig 4. The centre image displays periodic circles in the central region while the exterior is fronded, as in the left image of fig 19. You can also use tracepq() and then tacepqs() with option 2, 3 and nonfree() to produce limit sets in the same manner as fig 18 right.

These figures show the same set of Fibonacci fractions as in fig 15 trending from periodicity to non-free degeneracy. The circles are still evident in the centre image at the crossings.

The stochastic algorithm **qFstochastic()** is inferior at depicting the genuine cases above where we know that the limit set is a fractal quasi-circle of a genuint quasi-Fuchsian group or pinched examples like the Apollonian gasket, but it cones into its own when trying to explore the boundaries at the edge of degeneracy or dive into the chaotic regime where the depth first algorithm breaks down because it is assuming it can link lines scanning addresses in a circular search tree.

The function can depict limit sets in various parametrizations including the Jorgensen and Maskit versions.

The algorithm works by firstly finding attracting fixed points of the transofrmations, then choosing a partial orbit of the commutator and then repeatedly iterating these points in a random sequence of transfrmations *a, b, a*^{-1} and *b*^{-1}, collecting all these points and then plotting them in the same way as the above program. Somewhat paradoxically for those used to chaotic iterations of Julia sets, the random application of the four transforms preserves the limit set, but it fails to distribute evenly, visiting some points exponentially rarely. Nevertheless it does give an accurate depiction of a couldy version of the limit set and works irrespective of whether the example is quasi-Fuchsian, degenerate or chaotic, so one can explore a wider class of sets and use the algorithm t find new examples to apply to an image via the depth search algorithm.

The utility of the stochastic algorithm is emphasized by its ability to depict atypical limit sets such as the non-discrete examples below arising from the trace values *Ta *= 1 and *Ta* = 2^{1/2}, with *a *^{6} = I and *a* ^{8} = I respectively even under variations of the complementary trace, including two doubly non-free examples where *(aba*^{-1}*b*^{-1})^{2}^{ }= I as well.

Below are examples where the parameters have been chosen to take the set outside the regime of discrete groups, so that the limit set becomes unstable or fractal cantor dust, which can include akrbitrary commutator traces, or non-free situations requiring automatic groups for the depth first search algorithm.

The series of images below show how the limit sets vary with the trace of the commutator, when both *Ta* and *Tb *are* 2.* While only the values for -2 and 0 can be plotted using depthfirst() the stochastic algorithm can freely explore parameter space.

**Escape-time algorithm** **Download m-file here**

The orbital algorithm from Jos Leys web site works only in Maskit parametrization and only with tractable examples with Tb=2. The method scans all points in the field and iterates. The point is iterated by moving it via translations by 2 in the parallellograms bounded by the sloping line defined by the real and imaginary parts of Ta, until it ends in the centre section, then iterating it by a if it lies above the elliptical curve chosen below to separate the regions, or a if it is below, again translating if necessary and iterating again. If the point goes out the top, we colour it brown, if below, we colour it green and of it ends in a 2-cycle we colour it blue, separating the regions as shown in the diagram. It thus doesn't draw the limit set but the points on either side, highlighting the limit set by default as in level set methods for Julia sets. This works only for a few easily separated examples and require intial curve fitting of the separating sigmoid curve's slope and exponent.

**quasi3D** is a very simple example of a 3D transformation forming a limit set taken from Matlab central as the begining seed of the stochastic algorithm.

**Depthfirst Matlab Listing in Colour**

function outfl=depthfirst(A,B,ep,levmax,scal,sav,spec,siz,tab)

%depthfirst(2,2,0.01,1000,1,0);
%Standard option=0

%depthfirst(2,2,0.01,1000,1,0,1);
%Option=1 closing gaps

%depthfirst(1.87+0.1i,1.87-0.1i,0.001,1000,1/2.8,0,0);

%depthfirst(1.87+0.1i,1.87-0.1i,0.01,1000,1/2.8,0,1);
special words

%Example 11.1 non-free group
with relations (aba^-1b^-1)^2=I

%depthfirst(1.924781-0.047529i,2,0.005,1000,2,1,2,1100,1.9248-1.4617i);

%depthfirst(1.9247306-0.0449408i,
1.91+0.2i,0.005,3000,1/6,1,2,1100, 1.8569-1.2426i);

%depthfirst(2,2,0.005,1000,2,2,2,1100,2.0000-1.4142i);

%or
depthfirst(2,2,0.005,1000,2,2,2,1100,nonfree(2,2,0,0));

%Special words example 8.16
with aaB parabolic.

%the function speacialwords can
be reprogrammed to handle examples like

%the double cusps where longer
expressions like 2/5 where aaaBaaB is parabolic

%depthfirst2b(2+i,3,0.005,3000,1,2,4)

% A,B are traces of the Mobius
transforms generating a,b,a^-1,b^-1

% ep=epsilon in the tree search
using commutator attracting fixed pts

% levmax is the maximum depth
of the tree search (10000 if set to 0)

% scal is the scaling factor
for depicting the image

% sav=1,3 saves the file as a
png. if sav>=2 we also have a program timer

% sav=-1 test mode - exits on
completion outfl=1 or on reaching levmax outfl=0

% the timer reports tabs(1-4)
as the numbers cycle through branches

% the plot is black and white
but the png is coloured by the first iterate

% colfac and iterate can be set
in lines 53-55 to color by the ith iterate

% spec=0 standard method

% spec=1 special words used
also on fixed pts of a,b,a^-1,b^-1

% this refines whorls with
necks and when A=2 and B=2 closes the gasket gaps

% spec=2,3 as 0,1 but for non-free groups using automatic group
table.

% FSA is used to remove null
words when (aba^-1b^-1)^2=I

% spec=4 special words included
from the subfunction specialwords()

% siz is the image size default
1000x1100

% tab is the trace of A used in
the advanced grandma algorithm

% to test full exploration of
the tree uncomment 102-107, with levmax=4

outfl=true;

if nargin<8

siz=1100;

end

if nargin<7

spec=0;

end

if nargin<6

sav=0;

end

if nargin<5

scal=1;

end

ep=ep/scal;

tic;

clf;

hold on

M=255*ones(siz);

if nargin<9

[a,b]=grandma(A,B,0,0);

else

[a,b]=grandma(A,B,0,2,tab);

end

%[a,b]=grandma(A,B,0);

colfac=1; %Set to 0 for a black and white image 1 for colour-coded

iterate=1; %Which iterate to colour code by

gens=cell(1,4);

gens{1}=a;

gens{2}=b;

gens{3}=a^-1;

gens{4}=b^-1;

if levmax==0;

levmax=10000;

end

word=cell(1,levmax);

tags=zeros(1,levmax);

state=zeros(1,levmax);

FSA=[2 3 4 5;2 6 0 5;8 3 9 0;0 3 4
7;10 0 11 5;8 3 12 0;13 0 11 5;2 6 0 14;0 3 4 15;2 16 0 5;0 17 4 5;0 3 4 18;2
19 0 5;10 0 0 5;0 0 11 5;8 3 0 0;0 3 9 0;0 0 11 5;8 3 0 0];

%FSA=[2 3 4 5;2 3 0 5;2 3 4 0;0
3 4 5;2 0 4 5]; %table replicating std exs.

stfl=true;

oldpoint=0;

if spec<4

fix=zeros(3,4);

for i=1:4

fix(1,i)=fixp(gens{mymod(i-1)}*gens{mymod(i-2)}*gens{mymod(i-3)}*gens{mymod(i-4)});

fix(2,i)=fixp(gens{i});

fix(3,i)=fixp(gens{mymod(i+1)}*gens{mymod(i+2)}*gens{mymod(i+3)}*gens{mymod(i+4)});

end

else

[numfp,fix]=specialwords1(a,b);

points=zeros(max(numfp));

end

%stpt=0;

lev=1;

tags(1)=1;

state(1)=1;

word{1}=gens{1};

dotfl=false;

while ~(lev==0&&tags(1)==2)

fl=false;

if spec==2||spec==3

if state(lev)==0

fl=true;

end

end

while ~fl

plt=false;

switch spec

case {1,3}

newpoint=mob_on_point(word{lev},fix(1,tags(lev))); %branch termination?

oldpoint=mob_on_point(word{lev},fix(3,tags(lev))); %branch termination?

midpoint=mob_on_point(word{lev},fix(2,tags(lev)));

if ((lev==levmax) || (abs(newpoint-midpoint)<ep &&
abs(oldpoint-midpoint)<ep))

plt=true;

end

case {0,2}

newpoint=mob_on_point(word{lev},fix(1,tags(lev))); %branch termination?

oldpoint=mob_on_point(word{lev},fix(3,tags(lev))); %branch termination?

if ((lev==levmax) ||
(abs(newpoint-oldpoint)<ep))

plt=true;

end

case 4

pts=numfp(tags(lev));

fl4=true;

pts=numfp(tags(lev));

points(1)=mob_on_point(word{lev},fix(1,tags(lev)));

for i=2:pts

points(i)=mob_on_point(word{lev},fix(i,tags(lev)));

if abs(points(i-1)-points(i))>=ep

fl4=false;

break;

end

end

if lev==levmax

fl4=true;

end

if fl4

plt=true;

newpoint=points(1);

end

end

if lev==levmax

if sav>=2

fprintf('.');

dotfl=true;

end

if sav==-1

outfl=false;

return;

end

end

x=scal*abs(newpoint);

if x>50 %breakout
to portray chaotic limit sets

fprintf('Exiting...\n');

colormap(colorcube);

image(255*M');

return;

end

if spec<4 && x>1.6 && plt==true %too far offscreen, so break

fl=true;

break;

end

if plt

switch spec

case {1,3}

if A==2 && B==2

M=myplot(scal,M,oldpoint,midpoint,siz,colfac*tags(iterate),spec);

M=myplot(scal,M,midpoint,newpoint,siz,colfac*tags(iterate),spec);

else

M=myplot(scal,M,oldpoint,newpoint,siz,colfac*tags(iterate),spec);

end

case {0,2}

M=myplot(scal,M,oldpoint,newpoint,siz,colfac*tags(iterate),spec);

case 4

for i=2:pts

M=myplot(scal,M,points(i-1),points(i),siz,colfac*tags(iterate),spec);

end

end

if ~stfl

if spec==2 || spec==3

if abs(oldpoint-stpoint)<10*ep

M=myplot(scal,M,oldpoint,stpoint,siz,colfac*tags(iterate),spec);

end

end

end

stpoint=newpoint;

fl=true;

else

fl=false;

end

if fl

if stfl

%stpoint=oldpoint;

stfl=false;

end

break;

else

lev=lev+1; %go forward

tags(lev)=mymod(tags(lev-1)+1);

if lev==1

word{lev}=gens{tags(lev)};

if spec==2||spec==3

state(lev)=tags(lev);

end

else

word{lev}=word{lev-1}*gens{tags(lev)};

if spec==2||spec==3

state(lev)=FSA(state(lev-1),tags(lev));

if state(lev)==0

fl=true;

end

end

end

end

end

fl=true;

while fl

lev=lev-1; %go backward

if lev==0 ||
tags(lev+1)-1~=mod(tags(lev)+2,4) %available
turn

fl=false;

end

end

if (lev==0 && tags(1)==2)

%
M=myplot(scal,M,stpt,oldpoint,siz,colfac*tags(iterate),spec);

if dotfl

fprintf('\n');

dotfl=false;

end

fprintf('%d %d %d %d %4.4f\n', tags(1),tags(2),tags(3),tags(4),toc);

toc;

break; %Finished - exit routine

end

lev=lev+1;

tags(lev)=mymod(tags(lev)-1); %turn
and go forward

if sav>=2

if lev<4 %progress
clock for long iterations

if dotfl

fprintf('\n');

dotfl=false;

end

fprintf('%d %d %d %d %4.4f\n', tags(1),tags(2),tags(3),tags(4),toc);

end

end

if lev==1

word{1}=gens{tags(1)};

%state(1)=tags(1);

else

word{lev}=word{lev-1}*gens{tags(lev)};

if spec==2||spec==3

state(lev)=FSA(state(lev-1),tags(lev));

% if state(lev)==0

% fl=true;

% end

end

end

end

s=size(M); % Plot the matrix and save a png file

tit=[num2str(A),', ',num2str(B),', ',num2str(ep),', ',num2str(levmax)];

if spec==2 || spec==3

tit=[tit,', ',num2str(tab)];

end

z=zeros(1,s(1)*s(2));

k=1;

sc11=11/10;

for i=1:s(1)

for j=1:s(2)

if M(i,j)<255

z(k)=((i-s(1)/2)*2/s(1)+1i*(j-s(2)/2)*2/s(2))*sc11;

k=k+1;

end

end

end

plot(z(1:k-1),'.','MarkerSize',1,'Color','k');

axis([-sc11 sc11 -sc11 sc11]);

text(-1,1.05,tit);

if mod(sav,2)

% if colfac==0

% imwrite(255*M',[tit,'.png']);

% else

imwrite(M',colorcube(256),[tit,'.png']);

% end

%save([tit2,'.mat'],'M');

end

function n=mymod(m)

n=mod(m,4);

if n==0

n=4;

end

function M=myplot(sc,M,o,n,siz,c,spec) %Plots a scaled image in M

siz=siz/2;

sizs=10/11*siz;

fl=true;

big=10/sc;

if (abs(o))>big || abs(n)>big %added to avoid plotting line offscreen

fl=false;

end

switch c

case 0

co=241; %black

case 4

co=186; %magenta

case 1

co=233; %blue

case 2

co=208; %red

case 3

co=165; %cyan

end

% if spec==2 &&
abs(n-o)>0.1/sc %added to remove four extraneous lines

% fl=false;

% end

if fl

if abs(imag(n-o))>abs(real(n-o))

stps=max(ceil(abs(sc*imag(n-o)*siz)),2);

zs=linspace(sc*o,sc*n,stps);

else

stps=max(ceil(abs(sc*real(n-o)*siz)),2);

zs=linspace(sc*o,sc*n,stps);

end

for i=1:length(zs)

if ~isnan(zs(i))

if max(abs(real(zs(i))),abs(imag(zs(i))))<1

M(siz-floor(sizs*real(zs(i))),siz-floor(sizs*imag(zs(i))))=co;

end

end

end

end

function zo=mob_on_point(T,z) %Mobius
transforms a point

%D=struct();

a=T(1,1);

b=T(1,2);

c=T(2,1);

d=T(2,2);

zo=(a*z+b)/(c*z+d);

function z=fixp(T) %Finds an
attracting fixed point of a transformation

a=T(1,1);

b=T(1,2);

c=T(2,1);

d=T(2,2);

T2=T^2;

tr=T(1,1)+T(2,2);

tr2=T2(1,1)+T2(2,2);

k=((tr+sqrt(tr2-4))/2)^2;

if abs(k)>1

z=(-(d-a)+((d-a)^2+4*b*c))^(1/2)/(2*c);

else

z=(-(d-a)-((d-a)^2+4*b*c))^(1/2)/(2*c);

end

function [a,b]=grandma(ta,tb,posroot,opt,tab)

if nargin<5

opt=mod(opt,2);

end

if nargin<4

opt=0;

end

if nargin<3

posroot=false;

end

if opt<2 %Calculate
the trace Tab from the grandfather identity

p=-ta*tb;

q=ta^2+tb^2;

if posroot

tab=(-p+(p^2-4*q)^0.5)/2;

else

tab=(-p-(p^2-4*q)^0.5)/2;

end

z0=((tab-2)*tb)/(tb*tab-2*ta+2*1i*tab);

if opt==0 %Grandma's original recipe

ab=[tab/2
(tab-2)/(2*z0);(tab+2)*z0/2 tab/2];

b=[(tb-2*1i)/2 tb/2;tb/2
(tb+2*1i)/2];

a=ab*b^-1;

end

if opt==1 %Jorgensen's recipe

a=[ta-tb/tab ta/tab/tab;ta
tb/tab];

b=[tb-ta/tab -tb/tab/tab;-tb
ta/tab];

end

else

tC=ta^2+tb^2+tab^2-ta*tb*tab-2; %Grandma's three parameter double alarm

Q=sqrt(2-tC);

if abs(tC+1i*Q*sqrt(tC+2))>=2+1e-15 %e-15 added to avoid Ta = Tb =2 error

R=sqrt(tC+2);

else

R=-sqrt(tC+2);

end

%R=-R;

zo=(tab-2)*(tb-R)/(tb*tab-2*ta+1i*Q*tab);

a=[ta/2
(ta*tab-2*tb+2i*Q)/((2*tab+4)*zo);(ta*tab-2*tb-2i*Q)*zo/(2*tab-4) ta/2];

b=[(tb-1i*Q)/2
(tb*tab-2*ta-1i*Q*tab)/((2*tab+4)*zo);(tb*tab-2*ta+1i*Q*tab)*zo/(2*tab-4)
(tb+1i*Q)/2];

end

function [numfp, fix]=specialwords1(a,b)

A=a^-1;

B=b^-1;

numfp=zeros(1,4);

for i=1:4

switch i

case {1,3}

numfp(i)=4;

case{2,4}

numfp(i)=3;

end

end

fix=zeros(4,4);

fix(1,1)=fixp(b*A*B*a);

fix(2,1)=fixp(a*B*a);

fix(3,1)=fixp(B*a*a);

fix(4,1)=fixp(B*A*b*a);

fix(1,2)=fixp(A*B*a*b);

fix(2,2)=fixp(A*A*b);

fix(3,2)=fixp(a*B*A*b);

fix(1,3)=fixp(B*a*b*A);

fix(2,3)=fixp(A*b*A);

fix(3,3)=fixp(b*A*A);

fix(4,3)=fixp(b*a*B*A);

fix(1,4)=fixp(a*b*A*B);

fix(2,4)=fixp(a*a*B);

fix(3,4)=fixp(A*b*a*B);

**Tracepq Program Listing in Colour**

function [e,tt]=tracepq(p,q,opt,prt)

%opt=0 f==2 root 1 f==-2 2 non
free f==2 3 non free f==-2

syms t

g=gcd(p,q);

p=p/g;

q=q/g;

s=fareyseq(p,q);

t0=t;

if opt<2

t1=t+2i;

else

t1=t+sqrt(2)*1i;

end

t2=2;

for i=1:length(s)

if s(i)==0

h=expand(t1*t0-t2);

t2=t1;

t1=h;

else

h=expand(t1*t2-t0);

t0=t1;

t1=h;

end

%fprintf('.');

end

%fprintf('\n');

c=expand(t1); %Expand it into polynomial form

f=symfun(c,t); %Make it into a function we can solve

if mod(opt,2)==0

e=eval(vpasolve(f==2,t)); %Solve it for f=2

else

e=eval(vpasolve(f==-2,t));

end

q=q/gcd(p,q);

tt=zeros(length(e),3);

for i=1:q

tt(i,1)=i;

tt(i,2)=real(e(i));

tt(i,3)=imag(e(i));

end

tt=sortrows(tt,2);

if prt

for i=1:q

if tt(i,3)>=0

fprintf('%d\t%4.8f+%4.8fi\n',tt(i,1),tt(i,2),tt(i,3));

else

fprintf('%d\t%4.8f%4.8fi\n',tt(i,1),tt(i,2),tt(i,3));

end

end

end

function s=fareyseq(p,q)

p=p/gcd(p,q);

q=q/gcd(p,q);

s=[];

p0=0;

q0=1;

p1=1;

q1=0;

ph=p0+p1;

qh=q0+q1;

while p~=ph || q~=qh

if ph/qh>p/q

p1=ph;

q1=qh;

ph=ph+p0;

qh=qh+q0;

s=[s
0];

else

p0=ph;

q0=qh;

ph=ph+p1;

qh=qh+q1;

s=[s
1];

end

end

**Depth First C-script Listing in Colour**

//To run this file
from terminal put it in a folder, cd to the folder and input

//gcc -o depthf
depthfirst.c -lm && ./depthf

//To run successive
compiled versions type ./depthf

#include <stdio.h>

#include <math.h>

#include <complex.h>

#include <time.h>

typedef double _Complex arry[2][2];

void mult(arry a,arry b, arry c);

void inv(arry a, arry c);

double _Complex fixp(arry T);

time_t start,end;

int main(void)

{

FILE *mf;

const char base[] = "DF";

char filename [100];

//double _Complex ta=1.9264340533-0.0273817919*I,tb=2;

//double _Complex ta=1.80785524+0.136998688*I, tb=2;

//double _Complex ta = 1.924781-0.047529*I, tb=2, tab; //= 1.9248-1.4617*I //spec=2

//double _Complex ta = 2, tb=2; //=2.0000-1.41421356*I; //spec=2

//double _Complex ta = 1.87+0.1*I, tb = 1.87-0.1*I;

//double _Complex ta = 1.61891069-0.80846358*I,
tb=1.66692448-0.40921981*I;

double _Complex ta=1.95824516-0.01663690*I, tb=2;//spec=2

//double _Complex ta = 1.5-0.8660254038*I,
tb=1.5+0.8660254038*I;

int spec=2,
altroot=0; //spec=2 non-free
group altroot sets which sqare root for tab

int neck=1,iterate=0; //to aid drawing large whorls with narrow necks
usinf midpoint

const int levmax=1000,siz=1100;

double ep=0.001,scal=2;

double _Complex tab;

int M[siz][siz],siz2,sizs;

double _Complex gens[2][2][4];

double _Complex word[3][3][levmax+1];

int tags[levmax+1];

int state[levmax+1];

double dif;

int FSA[19][4]={{1,2,3,4},{1,5,-1,4},{7,2,8,-1},{-1,2,3,6},{9,-1,10,4},{7,2,11,-1},{12,-1,10,4},{1,5,-1,13},{-1,2,3,14},{1,15,-1,4},{-1,16,3,4},{-1,2,3,17},{1,18,-1,4},{9,-1,-1,4},{-1,-1,10,4},{7,2,-1,-1},{-1,2,8,-1},{-1,-1,10,4},{7,2,-1,-1}};

double _Complex fix[3][4];

int i,j,co;

double _Complex p,q,z0,z,Q,tC,R,zo;

arry a,b,ab,A,B,h;

arry *iv, *mu,h1,h2,h3;

int stfl,lev,fl,plt,stps;

double _Complex oldpoint,stpt,newpoint,midpoint;

double x;

sprintf(filename, "%s %e %e %e %e %d %e %e.txt",base,creal(ta),cimag(ta),creal(tb),cimag(tb),levmax,ep,scal);

mf=fopen(filename,"w");

time (&start);

ep=ep/scal;

for (i=0;i<siz;i++)

for (j=0;j<siz;j++)

M[i][j]=255;

siz2=siz/2;

sizs=(int)siz2*10.0/11;

//grandma's recipe

if (spec==0)

{

p=-ta*tb;

q=ta*ta+tb*tb;

if (altroot==0)

tab=(-p-csqrt(p*p-4*q))/2;

else

tab=(-p+csqrt(p*p-4*q))/2;

//tab=(-p-csqrt(p*p-4*q))/2;

z0=((tab-2)*tb)/(tb*tab-2*ta+2*I*tab);

ab[0][0]=tab/2;

ab[0][1]=(tab-2)/(2*z0);

ab[1][0]=(tab+2)*z0/2;

ab[1][1]=tab/2;

b[0][0]=(tb-2*I)/2;

b[0][1]=tb/2;

b[1][0]=tb/2;

b[1][1]=(tb+2*I)/2;

}

else

{

p=-ta*tb;

q=ta*ta+tb*tb-2;

if (altroot==0)

tab=(-p-csqrt(p*p-4*q))/2;

else

tab=(-p+csqrt(p*p-4*q))/2;

tC=ta*ta+tb*tb+tab*tab-ta*tb*tab-2; //Grandma's three parameter double alarm

Q=csqrt(2-tC);

if ((cabs(tC+I*Q*csqrt(tC+2))>=2+1e-15))

R=csqrt(tC+2);

else

R=-csqrt(tC+2);

zo=(tab-2)*(tb-R)/(tb*tab-2*ta+1i*Q*tab);

a[0][0]=ta/2;

a[0][1]=(ta*tab-2*tb+2*I*Q)/((2*tab+4)*zo);

a[1][0]=(ta*tab-2*tb-2*I*Q)*zo/(2*tab-4);

a[1][1]=ta/2;

b[0][0]=(tb-I*Q)/2;

b[0][1]=(tb*tab-2*ta-I*Q*tab)/((2*tab+4)*zo);

b[1][0]=(tb*tab-2*ta+I*Q*tab)*zo/(2*tab-4);

b[1][1]=(tb+I*Q)/2;

}

inv(b,B);

if (spec==0)

mult(ab,B,a);

inv(a,A);

for (i=0;i<2;i++) for (j=0;j<2;j++) gens[i][j][0]=a[i][j];

for (i=0;i<2;i++) for (j=0;j<2;j++) gens[i][j][1]=b[i][j];

for (i=0;i<2;i++) for (j=0;j<2;j++) gens[i][j][2]=A[i][j];

for (i=0;i<2;i++) for (j=0;j<2;j++) gens[i][j][3]=B[i][j];

mult(b,a,h1); mult(A,h1,h2);
mult(B,h2,h3); fix[0][0]=fixp(h3);

mult(A,b,h1); mult(B,h1,h2);
mult(a,h2,h3); fix[0][1]=fixp(h3);

mult(B,A,h1); mult(a,h1,h2);
mult(b,h2,h3); fix[0][2]=fixp(h3);

mult(a,B,h1); mult(b,h1,h2);
mult(A,h2,h3); fix[0][3]=fixp(h3);

fix[1][0]=fixp(a);
fix[1][1]=fixp(b); fix[1][2]=fixp(A);
fix[1][3]=fixp(B);

mult(B,a,h1); mult(A,h1,h2);
mult(b,h2,h3); fix[2][0]=fixp(h3);

mult(a,b,h1); mult(B,h1,h2);
mult(A,h2,h3); fix[2][1]=fixp(h3);

mult(b,A,h1); mult(a,h1,h2);
mult(B,h2,h3); fix[2][2]=fixp(h3);

mult(A,B,h1); mult(b,h1,h2);
mult(a,h2,h3); fix[2][3]=fixp(h3);

lev=0;

tags[0]=0;

state[0]=0;

for (i=0;i<2;i++) for (j=0;j<2;j++) word[i][j][0]=gens[i][j][0];

while (~((lev==-1)&&(tags[0]==1)))

{

fl=0;

if (spec==2)

if (state[lev]==-1)

fl=1;

while (fl==0)

{

z=fix[0][tags[lev]];

newpoint=(word[0][0][lev]*z+word[0][1][lev])/(word[1][0][lev]*z+word[1][1][lev]);

z=fix[2][tags[lev]];

oldpoint=(word[0][0][lev]*z+word[0][1][lev])/(word[1][0][lev]*z+word[1][1][lev]);

if (neck>0)

{

z=fix[1][tags[lev]];

midpoint=(word[0][0][lev]*z+word[0][1][lev])/(word[1][0][lev]*z+word[1][1][lev]);

}

plt=0;

x=scal*cabs(newpoint);

if (x>10) //Chaotic
instability

{

printf("exiting...\n");

return(0);

}

if (neck>0)

{

if ((lev==levmax) ||
((cabs(newpoint-midpoint)<ep)&&(cabs(midpoint-oldpoint)<ep)))

{

if (x>1.5) //Limit
set off screen

{

fl=1;

break;

}

plt=1;

}

}

else

{

if ((lev==levmax) || (cabs(newpoint-oldpoint)<ep))

{

if (x>1.5) //Limit
set off screen

{

fl=1;

break;

}

plt=1;

}

}

if (plt==1)

{

fl=1;

z=newpoint-oldpoint;

if (fabs(cimag(z))>fabs(creal(z)))

{

stps=1;

x=ceil(fabs(scal*cimag(z)*siz2));

if (x>1) stps=x;

}

else

{

stps=1;

x=ceil(fabs(scal*creal(z)*siz2));

if (x>1) stps=x;

}

for (j=0;j<=stps;j++)

{

z=scal*(newpoint*j/stps+oldpoint*(stps-j)/stps);

if ((fabs(creal(z))<1.0) && (fabs(cimag(z))<1.0))

{

if (iterate==0)

M[(int)(siz2-floor(sizs*creal(z)))][(int)(siz2-floor(sizs*cimag(z)))]=241;

if (iterate>0)

{

if (tags[iterate-1]==0) co=233; //blue

if (tags[iterate-1]==1) co=208; //red

if (tags[iterate-1]==2) co=165; //cyan

if (tags[iterate-1]==3) co=186; //magenta

M[(int)(siz2-floor(sizs*creal(z)))][(int)(siz2-floor(sizs*cimag(z)))]=co;

}

if (iterate==-1)

{

if (tags[lev]==0) co=233; //blue

if (tags[lev]==1) co=208; //red

if (tags[lev]==2) co=165; //cyan

if (tags[lev]==3) co=186; //magenta

M[(int)(siz2-floor(sizs*creal(z)))][(int)(siz2-floor(sizs*cimag(z)))]=co;

}

}

}

}

else

fl=0;

if (fl==1)

{

}

else

{

lev=lev+1; //go
forward

tags[lev]=(tags[lev-1]+1)%4;

if (lev==0)

{

for (i=0;i<2;i++) for (j=0;j<2;j++)

word[i][j][lev]=gens[i][j][tags[lev]];

if (spec==2) state[lev]=tags[lev];

}

else

{

for (i=0;i<2;i++) for (j=0;j<2;j++)

{

h1[i][j]=word[i][j][lev-1];

h2[i][j]=gens[i][j][tags[lev]];

}

mult(h1,h2,h3);

for (i=0;i<2;i++) for (j=0;j<2;j++)

word[i][j][lev]=h3[i][j];

if (spec==2)

{

state[lev]=FSA[state[lev-1]][tags[lev]];

if (state[lev]==-1) fl=1;

}

}

}

}

fl=1;

while (fl==1)

{

lev=lev-1; //go
backward

if ((lev==-1) || ((tags[lev+1]+3)%4 !=
((tags[lev]+2)%4))) //available
turn

{

fl=0;

}

}

if ((lev==-1)&&(tags[0]==1))

{

printf("%d %d %d %d\n", tags[0]+1,tags[1]+1,tags[2]+1,tags[3]+1);

time (&end);

dif=difftime (end,start);

printf("Elasped time is %.0lf
seconds.\n", dif );

break; //finished!!

}

lev=lev+1;

if (tags[lev]>0)

tags[lev]=(tags[lev]-1)%4; //turnandgoforward

else

tags[lev]=(tags[lev]+3)%4;

if (lev<3) //progress clock for long iterations

{

time (&end);

dif=difftime (end,start);

printf("%d %d %d %d ", tags[0]+1,tags[1]+1,tags[2]+1,tags[3]+1);

printf("%.1lf sec\n", dif );

}

if (lev==0)

{

for (i=0;i<2;i++) for (j=0;j<2;j++)

{

word[i][j][lev]=gens[i][j][tags[lev]];

//state[0]=tags[0];

}

}

else

{

for (i=0;i<2;i++) for (j=0;j<2;j++)

{

h1[i][j]=word[i][j][lev-1];

h2[i][j]=gens[i][j][tags[lev]];

}

mult(h1,h2,h3);

for (i=0;i<2;i++) for (j=0;j<2;j++)

word[i][j][lev]=h3[i][j];

if (spec==2)

{

state[lev]=FSA[state[lev-1]][tags[lev]];

lev=lev+1-1;

// if (state[lev]==-1 fl=1;

}

}

}

for(i=0;i<siz;i++)

for(j=0;j<siz;j++)

fprintf(mf,"%d ",M[i][j]);

fclose(mf);

return(0);

}

double _Complex fixp(arry T)

{

double _Complex z,tr,tr2,a,b,c,d;

double _Complex k;

a=T[0][0];

b=T[0][1];

c=T[1][0];

d=T[1][1];

tr=a+d;

tr2=a*a+2*b*c+d*d;

k=((tr+csqrt(tr2-4))/2)*((tr+csqrt(tr2-4))/2);

if (cabs(k)>1)

z=(-(d-a)+csqrt((d-a)*(d-a)+4*b*c))/(2*c);

else

z=(-(d-a)-csqrt((d-a)*(d-a)+4*b*c))/(2*c);

return(z);

}

void mult(arry a,arry b,arry c)

{

int i,j;

for (i=0;i<2;i++)

for (j=0;j<2;j++)

c[i][j]=a[i][0]*b[0][j]+a[i][1]*b[1][j];

}

void inv(arry a, arry c)

{

double _Complex d;

d=a[0][0]*a[1][1]-a[0][1]*a[1][0];

c[0][0]=a[1][1]/d;

c[0][1]=-a[0][1]/d;

c[1][0]=-a[1][0]/d;

c[1][1]=a[0][0]/d;

}

//To run this file
from terminal put it in a folder, cd to the folder and input

//gcc -o tracepqr
tracepqr.c -lm && ./tracepqr

//To run successive
compiled versions type ./depthf

#include <stdio.h>

#include <math.h>

#include <complex.h>

#include <time.h>

double _Complex fracfunct(int *s,int sl,int opt,double _Complex val);

int main(void)

{

//double _Complex val=2;

//int p=3,q=5,opt=0;

double _Complex val=1.61691788-1.29437374i;

int p=1597,q=2584,opt=0;

double ep=1e-8;

int g,hld,h1,h2,p0,q0,p1,q1,ph,qh,its,i;

int s[200],sl=0;

double _Complex vh,h,f,df;

h1=p; //Find gcd(p,q)

h2=q;

while (h2%h1>0)

{

hld=h2%h1;

h2=h1;

h1=hld;

}

g=h1;

p=p/g;

q=q/g;

p0=0; // Generate Farey sequence for p/q

q0=1;

p1=1;

q1=0;

ph=p0+p1;

qh=q0+q1;

while ((p!=ph) || (q!=qh))

if (ph*q>p*qh)

{

p1=ph;

q1=qh;

ph=ph+p0;

qh=qh+q0;

s[sl]=0;

sl=sl+1;

}

else

{

p0=ph;

q0=qh;

ph=ph+p1;

qh=qh+q1;

s[sl]=1;

sl=sl+1;

}

h=(1+I)*1e-7;

its=500;

for (i=0;i<its;i++) //Newton's iteration loop

{

vh=val;

f=fracfunct(s,sl,opt,val)-2;

df=(fracfunct(s,sl,opt,val+h)-fracfunct(s,sl,opt,val-h))/2/h;

val=val-f/df;

if (cabs(vh-val)<ep)

{

printf("OK
%d %4.12f %4.12f %4.12f\n",i,cabs(vh-val),creal(val),cimag(val));

break;

}

}

if (i==its)

printf("Bad:%d %4.12f %4.12f\n",i,creal(val),cimag(val));

}

double _Complex fracfunct(int *s, int sl,int opt,double _Complex val)

{

double _Complex t0,t1,t2,h;

int i;

t0=val;

if (opt<2)

t1=val+2i;

else

t1=val+sqrt(2)*1i;

t2=2;

for (i=0;i<sl;i++)

{

if (s[i]==0)

{

h=t1*t0-t2;

t2=t1;

t1=h;

}

else

{

h=t1*t2-t0;

t0=t1;

t1=h;

}

}

return(t1);

}