**Kleinian and Quasi-Fuchsian Limit Sets in Matlab
Genotype:
1.1.4**

Chris King

**Download** the **complete toolbox** with web resource

This page is a first edition of a suite of programs that depict limit sets of Kleinian and quasi-Fuchsian groups generated by a pair of Mobius transformations in the complex plane, which we will call *a *and *b*, along with their inverses *A = a*^{-1}* *and *B* = *b*^{-1}.

The aim is to provide a free 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 code could also readily br transferred to an open-source Octave format if desired. This makes the algorithm freely available to maths students and researchers, as well as the general public. These images seek to explore the systems themselves rather than produce extensive ray-traced artwork, which is widely covered elsewhere, but anyone seeking to produce 3D ray-traced versons could easily extend thsee algorithms within Matlab to do this.

There are four algorithms: (1) Depth first search with special words and automatic groups, (2) Stochastic and (3) Escape-time (4) A simple 3D iteration

**Depth First Search **Program listing in colour

The depth first search algorithm **depthfirst()** is the quintessance of depiction methods of quasi-Fuchsian limit sets and is the key instrument for studying diverse key examples. The quasi-code for the depth search algorithm below is derived from Indra's Pearls.

The Matlab function does a full depth search of all the key examples in the book, including the last example below, which requires 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 below 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.

Fig 1: The limit sets for real joint traces greater than 3 are a circle (Fuchsian). As the traces *Ta, Tb* reduce 3 through quasi-Fuchsian 2.2 (left) colour coded by the first iterate to 2 (centre), the fractal circle becomes a Kleinian Apollonian Gasket. The subset of circles shaded light grey are mapped from the largest left circle by the transformations *a *and *BAb *and never intersect the white circles, forming a modular subgroup of the gasket group. The quasi-Fuchsian limit set (right ) with traces *Ta* = 1.92 + 0.03*i*, *Tb* = 2 is colour coded by the second iterate. The left figure shows locations of some key elementary repeated elements.

The page also includes stochastic algorithms which can depict chaotic cases, a colour filling escape-time algorithm, and a 3-D example, but the DFS approach is key to the diverse mathematics of the limit sets, which becomes a marriage of group theory, 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.

Fig 1b: Diverse mathematics underlies the forms of the limit sets. The little topolgical drawings are from Indra's beads. (5) shows Schottky circles of non-kissing circles as in (2), defined by loxodromic maps on the Riemann sphere (1), taking the inside of one circle over the outside of the paired one, resulting in fractal dust Cantor sets as the limit sets. (6) shows the Riemann sphere with the two circle interiors removed ready for gluing in ways which can generate new surfaces. (7) is the construction for the modular group represented by the Escher tiling consisting of the group of matrices with determinant 1. (8) shows the arrangement for a quasi-Fuchsian group, where the four touching circles as in (3)are further drawn together, so that the surface becomes two punctured tori joined at their cusps and th elimit set becomes a fractal tpological circle, or Jordan curve. This is manifest algebraically in choosing mappings in which the trace* *of the commutator is parabolic, with *TabAB* = –2. As we move the circles, more can kiss, so that when all do, as in the gasket (4), two more transformations have become 'pinched' so that both *a *and b are parabolic and we end up with (9). The four tangent circles generating the Kleinian group of the Apollonian gasket partially portrayed in (4) are highlighted in black. One is the line *y* = 0, with centre the point at infinity. In many of the cases below with *Tb* = 2, only one generator *b* is parabolic, and in some such as fig 2, neither are, although the limit set of each quasi-Fuchsian limit set is still a fractal topological circle.

Mobius transformations are the analog of affine Euclidean transformations (translations, scalings, rotations) on the Riemann sphere (1) , 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 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. 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 fo the limit set algorithm, were the exact action of each group on the limit set becomes apparent, as in fig 1(left) and fig 1c. The traces of parabolic maps, defined by *p* + *s* (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 1b.Mobius mappings preserve circles and lines in the complex plane Lines are simply circles going through the north pole at infinity. When two Mobius transforrmations are combined, the actions of the matrices and their inverses generate a symmetry group. like the 3 rotations and three reflections of an equilateral tringle 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 (5), 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 (2) forming a fractal dust or Cantor set – the limit set of the group. As the circles come together and kiss or touch as in (3) the limit set can become a circle or a fractal quasi-circle consisting of 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 (3). The topology of the surface can undergo other pinching and puncturing transfomations as in (6), (7), and (8), resulting in Kleinian groups such as the Apollonian gasket (partial Schottky computation in 4) and limit set fig 1 centre, which is the limiting case of the 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 (9) 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 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 corckscrew-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.

Fig 1c: Special words and the words defining a frond of a limit set.

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 and secondly to avoid pinching off the whorls of limit sets with large appendages with narrow necks. 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. The image right illustrates the problem for the first example in fig 2, 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 requiring special words to get full 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.

The three examples in fig 9 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 degenerate, 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 degenerate word combinations.

Fig 2: Left and centre:* Ta *= 1.91 ± 0.05*i*, *Tb *= 1.91 ± 0.05*i* Notice that *Ta* and *Tb* have differing effects with *Ta* applying to the 'interior' around zero and *Tb* to the 'exterior' around infinity and that the conjugates (here + / +) can interact to produce surprising changes.
Right: *Ta *= 1.92 + 0.03*i*, *Tb* = 1.88 + 0.12*i*

In fig 2 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 other figures and is also about ten times larger.

Fig 3: Left:* Ta* = 1.91 + 0.05*i *, *Tb* = 3 Fig 8.9 238 Centre:* Ta *= 1.6421 − 0.76659*i*, *Tb *= 1.91+0.2*i* Right: *Ta *= 1.9656+0.99536*i*, *Tb *= 3 crop of spirals Fig 8.23 264

Fig 4: Near degenerate *Ta *= 1.7847+0.73852*i*, Tb=3 Fig 10.1 311 Space-filling in the plane *Ta *= 1.5 − 0.86603*i,* *Tb *= 1.5+0.86603*i, *i.e. (1±sqrt(3)*i*)/2 Fig 10.11 335 Degenerate *Ta *= 1.6169 − 0.70567*i*, *Tb* = 2 Fig 10.9 330

Group theorists take a special interest in the spectrum of properties of limit sets where one trace factor is 2, so below are several key examples.

Fig 5: Doubly degenerate Troell's point Jorgensen's group *Ta* = 1.6168866453 − 1.29432650*i**,* *Tb* = 2 Right: coloured by first iterate to highlight the boundary lines between the complementary discs left and right being consumed by the degeneracy Fig 10.4 318-319. This value corresponds to the minimum imaginary value on the boundary of the Maskit slice fig 8 centre and corresponds closely to the *μ* value of the golden mean. Actually it is the mirror image of the original from Indra's beads *Ta* = 1.6168866453 − 0.7056734968*i* printed for consistency with the book. See fig 8d for sequence convergence images to the golden mean.

Below 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 meanaing there is one (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. The resolution of the 7/43 case is a little too low to see all the smallest circles because of the limits on computation time but the counting is clear and you can also see the way the chain is running in cycles of seven corresponding to the numerator, and cycles of 2 i the centre image.

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 *ε*. We can also 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* =

Fig 5b: Left: Special word arrays for the parabolic element *aaB*. Right" Special words algorithm working on Ta = 2 + *i*, Tb = 2 8.16 247.

I have 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, such as those below 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 5b, I have included the one from Indra's Pearls 8.16 designed to show the method in the case of *aaB*. This ends up with the look up structure shown above 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, with the exception of the Schottky circles method. The images in fig 9 all correspond closely to those in Indra's Pearls, but the early stages of the centre example looks distinctly different in version 2b, at the same time as being an order of magnitude faster, because the linkages are defined in a slightly different way.

Fig 6: Double cusp 1/15 *Ta* = 1.958591030 + 0.011278560*i*, *Tb* = 2 both *b* and *a*^{15}*B *are parabolic Fig 9.1 269. In general solving for a given denominator *q* involves iteratively solving a *q*-th degree polynomial. Centre: 2/19 cusp *Ta*=1.9038 + 0.03958*i*, *Tb*=2 Fig 9.15 293 *a*^{10}*B*^{}*a*^{9}*B*^{} and *b* are parabolic. Right: 7/43 cusp *Ta*=1.80785524 + 0.136998688*i*, *Tb* = 2 Fig 9.16 294 − 295. You can see the denominators clearly in the periodic mappings of the circles, and the numeratos both in the steeping periods along the chain and in the chains between the horizontal circle pair and the vertical pair.

Fig 7: Left: Free discrete group *Ta *= 1.887+0.05*i,* *Tb *= 2 Fig 8.15 246 Centre: *a*^{3}*Ba*^{2}*B* is parabolic for *Ta *= 1.64213876 − 0.76658841*i*, as is *b* since *Tb *= 2 Fig 9.3 272. Right: Spirally degenerate *Ta *= 1.9264 − 0.027382*i*,* Tb* = 2 Fig 10.6 322-323.

Fig 8: In the centre is shown the "Maskit slice" - the plane of *μ*(*p*/*q*) where the generating *p*/*q* word, such as *a*^{15}*B* is parabolic from Indra's Pearls. Left and right are shown the values of -*i*(*μ*(*p*/*q*), 2 for two known values 1/2 and 3/10 with values 1 + sqrt(3)*i* and (1 + sqrt(11)*i*)/2. Notice the overlapping fractal dark circles in the right hand image. For an explanation of *μ* in realtion to the trace values, see fig 10.

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, only one or two of which will have a credible root close to 2 with a small imaginary value we can try. You can directly input the values such as s(7) into depthfirst as parameter A, giving full precision rather than the default 4-digit precision of the printed value. The function will happily solve equations up to 100, but values of very high denominators may nt portray correctly due to numerical errors in the root search. I had to spend a lot of time to find a large denominator which would compute well! 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 8 centre.

Fig 8b: Limit sets of* μ*(1/7), *μ*(1/15), *μ*(1/23), *μ*(1/31) pictured in blue with the gasket superimposed in black and *μ*(1/60) , show they are not tending to the limit of *μ*(1/*n*) = *μ*(0/1), as 1/*n* tends to 1/∞ = 0/1. This limit set is the Apollonian gasket whose trace *Ta* = 2 and matrix *a* is also the limit of the *μ*(1/*n*) traces and matrices. But the limit sets themselves tend in limit to another limit set altogether whose group contains an additional factor as well as that of the Apollonian gasket.You can see this with the gasket superimposed in green on the limit set of *μ*(1/31) (click to enlarge).

*Note that in the right superimposed figure for 1/31 in fig 8b, 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 6 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. 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, just as 2/19 does in fig 6 center.

We can also calculate the trace *Ta* for other cusp fractions because parabolic expressions like *a*^{3}*Ba*^{3}*Ba*^{2}*Ba*^{3}*Ba*^{2}*B* can also be generated recursively 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*)

This is all very fine for the Fibonacci numbers, but other fractions may require a little more fiddling. For example to get 2/3, we have 1/1 *TaB* = *Ta* +2*i,* 1/2 TaaB = *Ta** ^{2}* + 2

**Coefficient recursion** Program listing in colour

We can thus set up **recursecoeffs()** – an extension of the Matlab function above to find all the fractional *Ta*'s required to explore the dynamics. The function does this recursively, which is natural, because the grandfather identity defines a given trace polynomial in terms of three polynomials of lower fractional complexity, so we need to create a function which calls itself three times over aain and again until it has linked all the identities down to the root level. This is far less trivial than it might appear, because many ways of splitting the fraction lead to unfeasable solutions. In particular reducable fractions with common factors have to be excluded as they will make the polynomial degenerate, as well as improper fractions, which are out of [0,1]. More subtly, we need to make sure the way we are splitting ensures the words have compatable patterns of *aaaBaaB* etc, so that when we multiply the inverse *Tm*^{-1}*n, *the cancellations give us another fraction word, rather than a compound word that is not in the series of fractional words. The algorithm does this by splitting so that the two fraction reciprocals lie in the same integer step range. It works beginning with two smaller fractions that evenly divide the previous one and shuffles to try to avoid illegitimate decompositions, and has successfuly generated the traces for the fractions shown right that run up to the limits of computability on a personal computer. However, to get fractions of successive Fibonacci numbers accurately and efficiently, a fourth parameter fib provides the option to decompose in terms of previous Fibonacci sequence members.

Fig 8c: Three higher fractions with rotation rates close to 1/2, 1/3 and 1/4 10/19 1.70408769 – 1.01161324*i* 13/40 1.68667404 – 0.57810830*i* 10/39 1.71314918 – 0.35283835

Below is a sample session for depicting *μ*(4/9). The parameters in recursecoeffs(4,9,0,1) determine the fraction, option 0 is for standard rather than non-free groups, and 1 is to print out the fractional splits. The parameters in depthfirst select s(9) the 9th 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 standard group without special words for the generators (see listing), and 1100 is the image size. This will print out a series of stages of exploring the tree so you know progress is taking place for a detailed plot.

>> s=recursecoeffs(4,9,0,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

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

T(4/9)=T(1/2)*T(3/7)-T(2/5) 2.25 2.00 2.33

1
1.16618403 -1.05050570

2
1.18265280 -0.52328487

3
-0.78549275 -0.40614104

4
-0.44290892 -1.30596770

5
-1.45944586 -0.71925791

6
-1.59174796 -1.00910194

7
0.07875458 -1.82409535

8
0.19925605 -0.24362771

9
1.65274803 -0.91801778

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

Fig 8d: Periodicity limiting in degeneracy: Fibonacci fractions
* μ*(3/5)

The only noticeable difference between the right hand image in fig 8d and the limit set of Troell's point in fig 5 is right on the horozontal circles' boundary, where *μ*(89/144) still has infinitessimal hammer heads defining the horizontal circles..

Fig 8e: Left: Two sets with *Ta* = 1.8385 – 0.0964*i*, *Tb* = 1.7667 – 0.2201*i* and * Ta* = 1.8385 – 0.0964*i*, *Tb* = 1.7667 + 0.2201*i*. Right: Two limit set with *Ta*, *Tb* having values determined by * μ*(3/5)

The first two sets in fig 8c were generated by finding *Ta for Tb* = 2, using recursecoeffs() for 1/7 and then letting *Tb* = *Ta *and using a generalization recurscoeffs2() of the introductionary version above for *Tb* ≠ 2 to get a new *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. The right set is generated by combining the traces of * μ*(3/5)

Fig 9: Left: *Ta*=1.9247306-0.0449408*i*, *Tb*=1.91+0.2*i*, *Tab*=1.85692801+-1.24257380*i *11.3 363 Centre: *Ta* = 1.92478160 – 0.04752913*i*, *Tb* = 2 and *Tab* = 1.92480000-1.46174256*i* 11.1 354 *μ*(1/10) Right: *Ta*=*Tb*=2 *Tab*=2.0000-1.41421356*i *11.3 363

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.

Fig 9b: 1.865-0.70567i, 2 Spirally degenerate? *μ*(1/60) *Ta* = 1.9973 – 0.00025756*i*, showing the same limiting behaviour as in the standard version does not converge to the non-free gasket above. 1.906-0.06739i, 2 Near degenerate high period

We then have to set up an automaton 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 0 counting as a death state. 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 recursecoeffs() 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);

Fig 9c: Left: 1.92+0.05*i, *1.89+0.12*i*, Centre Left: 1.94+0.05*i*, 1.94+0.05*i* with similar non-free coefficients to the left hand images of fig 2. Centre right:Aa limit set generated from * μ*(3/5)

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 2. The centre image displays periodic circles in the central region while the exterior is fronded, as in the left image of fig 9. You can also use recursecoeffs() and then recurscoeffs2() with option 2, 3 and nonfree() to produce limit sets in the same manner as fig 8e right.

Fig 9d: Transition from periodicity to degeneracy: Non-free versions of *μ*(3/5), *μ*(8/13), *μ*(21/34) and *μ*(89/144) tending to the Golden Mean. For values see table above.

These figures show the same set of Fibonacci fractions as in fig 8d 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 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 function can depict limit sets in various parametrizations including the Jorgensen and Maskit versions.

Fig 10: Left: Standard, Jorgensen and Maskit representations of *Ta *= 1.91 + 0.1*i*, Tb = 2.
Research in the Maskit parametrization has focused on the two transforms *
a*(

Below are three examples where the parameters have been chosen to take the set outside the quasi-Fuchsian regime where the limit set becomes unstable or fractal cantor dust.

Fig 11: *Ta*=1.78+0.05*i*, *Tb*=3, *Ta*=1.6, *Tb*=1.6, and *Ta*=1.85+0.05*i*, *Tb*=1.85-0.2*i*

You can also plot examples requiring three trace parameters which are no longer quasi-circles, as shown left, and depict approximations to limit sets requiring automatic group representations such as the last example above, as shown right.

Fig 12: *Ta *= 1.991 + 0.05*i*,* Tb* = 1.93 - 0.05*i*, *Tab *= 2.2 - 1.8*i * *Ta *= 1.924781 - 0.047529*i*,* Tb *= 2, *Tab* = 1.9248-1.4617*i*

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.

Fig 12b: The function non-free() can also be used with a fourth parameter to find *Tab* for any value of the commutator trace *TabAB*. Above are shown a series of values spanning the real line from -2.5 to 1.8 showing the evolution of the limit sets towards chaotic instability.

**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.

Fig 13:* Ta* with real and imaginary parts 1.912, 1/30 and 1.74, 0.24

**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.

Fig 14: 3D Plots of quasi3D

**Depthfirst Program Listing in Colour**

function 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,1,2,1100,2.0000-1.4142i);

%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

% 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

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=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=0; %Set to 0 for a black and white image 0.041 for colour-coded

%slight changes to colfac alter
the colour scheme see: colormap(colorcube)

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};

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)<10/scal) ||
(abs(newpoint-midpoint)<ep && abs(oldpoint-midpoint)<ep))

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)<10/scal) ||
(abs(newpoint-oldpoint)<ep))

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 plt

% ep=0; %Test to check full tree us explored with ep=0

% for
ii=1:lev %set levmax ~ 4 for a quick view

% fprintf('%d ',tags(ii));

% end

% fprintf('\n');

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);

toc;

break; %Finished - exit routine

end

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

if sav>=2

if lev<4 %progress
clock for long iterations

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

end

end

if lev==0

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

level(1)=tags(1);

else

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

if spec==2||spec==3

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

end

end

lev=lev+1;

end

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

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

%tit2=[num2str(A),'-',num2str(B),'-',num2str(ep),',
',num2str(levmax)];

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)<1

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(255*M',colorcube,[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

% 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))))=c;

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 k>1

%if k>1 &&
abs(k)>1 fprintf('!'); end

if abs(k)<1

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

else

z=(-(d-a)-sqrt((d-a)^2+4*b*c))/(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
thre 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

R=sqrt(tC+2);

else

R=-sqrt(tC+2);

end

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);

**Recursecoeffs Program Listing in Colour**

function s=recursecoeffs(p,q,opt,prt,fib)

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

%e.g. s=recursecoeffs(7/43,0,1)

%Which you can follow by

%depthfirst(t(13),2,0.001,1000,1,3,0,1100);

%fib is designed to work
directly with two successive Fibonacci numbers

if nargin<5

fib=false;

end

if nargin<4

prt=false;

end

syms ta

opt1=(opt>=2);

t=tapoly(p,q,opt1,prt,fib); % Find the polynomial

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

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

if mod(opt,2)==0

s=eval(vpasolve(f==2,ta)); %Solve it for f=2

else

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

end

q=q/gcd(p,q);

for i=1:q

if abs(imag(s(i)))>0

fprintf('%d %4.8f %4.8f\n',i,real(s(i)),imag(s(i)));

end

end

function t=tapoly(p,q,opt,prt,fib) %Recursive polynomial definition

if nargin<4

prt=false;

end

if nargin<3

opt=false;

end

syms t

g=gcd(p,q);

if g>1

p=p/g;

q=q/g;

end

if p==1 %Calculate 1/q any q by
iteration

t=orderone(q,opt);

else

if q-p==1 %if we are at the other end e.g. 6/7 split to 1/1 and 5/6

p0=1;

q0=1;

p1=p-1;

q1=q-1;

else

if ~fib

p0=split(p,0); %Try a
near-even fraction split

p1=split(p,1);

q0=split(q,0);

q1=split(q,1);

else

[p0,q0,p1,q1]=findfracs(p,q);

end

st=0;

%Avoid common factors improper fractions and
incompatibile word ratios

while gcd(p0,q0)>1
|| gcd(p1,q1)>1 || p1-p0>q1-q0 || ~cf(q0/p0,q1/p1)

st=st+1;

if st==1

if q0>1 && q1<q %Adjust split fractions to find a good one

q1=q1+1;

q0=q0-1;

end

end

if st==2

if p0>1 && p1<q

q1=q1-1;

q0=q0+1;

p1=p1+1;

p0=p0-1;

else

st=0;

end

end

if st==3

if q0>1 && q1<q

q1=q1+1;

q0=q0-1;

end

end

if st==3

st=0;

end

end

end

tp=tapoly(p0,q0,opt,prt,fib); %Recurse
into tapoly 3 times

bt=tapoly(p1,q1,opt,prt,fib); %to
find the next fractions

if p1-p0==0 %to put in the grandfather
identity

if ~opt

df=tapoly(1,1,opt,prt,fib)-2i;

else

df=tapoly(1,1,opt,prt,fib)-sqrt(2)*1i;

end

% df=ta;

else

df=tapoly(p1-p0,q1-q0,opt,prt,fib);

end

t=tp*bt-df; %Apply the Grandfather identity

if prt

fprintf('T(%d/%d)=T(%d/%d)*T(%d/%d)-T(%d/%d)
%4.2f %4.2f %4.2f\n',p,q,p0,q0,p1,q1,p1-p0,q1-q0,q/p,q0/p0,q1/p1);

end

t=expand(t);

end

function [p0,q0,p1,q1]=findfracs(p,q)

p1=floor(p*q/(p+q));

q1=floor(q*q/(p+q));

p0=p-p1;

q0=q-q1;

while p0>=p1

p0=p0-1;

p1=p1+1;

end

while q0>=q1

q0=q0-1;

q1=q1+1;

end

if p0==0

p0=1;

p1=p1-1;

end

function a=split(c,typ) %Make one
of the four split test fractions

if c==2

a=1;

else

if mod(c,2)

if ~typ

a=floor(c/2);

else

a=ceil(c/2);

end

else

if ~typ

a=c/2-1;

else

a=c/2+1;

end

end

end

function t=orderone(stps,opt) %Iteratively
expand 1/n poly

%opt=0 normal 1 non-free

syms ta

a=cell(1,stps+1);

a{1}=2;

tic;

if ~opt

a{2}=ta+2i;

else

a{2}=ta+sqrt(2)*1i;

end;

for i=3:stps+1

a{i}=expand(ta*a{i-1}-a{i-2});

end

t=a{stps+1};

function f=cf(a,b) %Test if
the fractions have compatible word ratios

f=false;

if isint(a) && isint(b) && abs(b-a)==1

f=true;

end

if isint(a) && ~isint(b) && abs(b-a)<1

f=true;

end

if ~isint(a) && isint(b) && abs(b-a)<1

f=true;

end

if ~isint(a) && ~isint(b) && floor(b)==floor(a)

f=true;

end

function f=isint(a) %Test for
integer value

if ceil(a)==floor(a)

f=true;

else

f=false;

end