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.03i, 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-1NC 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-1b-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 bABa, 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 a15B 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-1b-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.05i, Tb = 1.91 ± 0.05i     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.03i, Tb = 1.88 + 0.12i 

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.05i , Tb = 3 Fig 8.9 238  Centre: Ta = 1.6421 − 0.76659i, Tb = 1.91+0.2i  Right: Ta = 1.9656+0.99536i, Tb = 3 crop of spirals Fig 8.23 264     


Fig 4: Near degenerate  Ta = 1.7847+0.73852i, Tb=3 Fig 10.1 311  Space-filling in the plane Ta = 1.5 − 0.86603i, Tb = 1.5+0.86603i, i.e. (1±sqrt(3)i)/2 Fig 10.11 335 Degenerate Ta = 1.6169 − 0.70567i, 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.29432650i, 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.7056734968i 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 a15B is in green and circles 1 and 2 are thus mapped to themselves. B carries C1 to C16 and then a15 carries it back to C1. In a similar way, B carries C2 to C17 and then a15 carries it back to C2. This means that a15B maps both the disks C0 and C1 to themselves, the fixed points of a15B 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 a15B for 1/15 or a3Ba3Ba2Ba3Ba2B 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-1n = TmTn and (2) the extended Markov identity TabAB = (Ta)2 + (Tb)2 + (Tab)2 Ta.Tb.Tab – 2. These form the founding identities to enable the decoding of the traces of the examples. We can then find Ta by solving a q-th degree polynomial e.g. by Newton's method, after solving for TaB and iteratively asserting constraints such as TaaB = ± 2.


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.011278560i, Tb = 2   both b and a15B 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.03958i, Tb=2  Fig 9.15 293  a10Ba9B and b are parabolic.  Right: 7/43 cusp  Ta=1.80785524 + 0.136998688i, 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.05i, Tb = 2 Fig 8.15 246 Centre: a3Ba2B is parabolic for Ta = 1.64213876 − 0.76658841i, as is b since Tb = 2 Fig 9.3 272. Right: Spirally degenerate Ta = 1.9264 − 0.027382i, 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 a15B 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 anB. 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 + 2i . Then using the grandfather identity on a and B , with TB = 2, we get a series of equations TanB = Ta.Tan-1B Tan-2B, making it possible to write a recursive process to define Ta as an nth 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 % Ta is most likely to be in the f=2 set.

    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.09628i, 1/15 1.958591 – 0.01128i, 1/23 1.98179 – 0.00320i, 1/31 1.98987 – 0.00131i, and 1/60 1.99726820 – 0.00018236, tending towards Ta = 2, or μ = 2i 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/(2n+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 a3Ba3Ba2Ba3Ba2B can also be generated recursively from a3B and a2B by similar relations. To be more precise, examination shows that the word w(p+r)/(q+s) = wp/r.wq/s, in the manner of mediants and the Farey tree, so we can use the grandfather identity as follows:
Tw(p+r)/(q+s) = Twp/r.Twq/sT(wp/r-1.wq/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) = Tw8/21 = Tw3/8.Tw5/13T(w3/8-1.w5/13) = Tw3/8.Tw5/13Tw2/5, because w3/8-1.w5/13 = (a3Ba3Ba2B)-1.a3Ba3Ba2Ba3Ba2B = a3Ba2B = w2/5.=w(5–3)/(13–8).

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 +2i, 1/2 TaaB = Ta2 + 2iTa - 2, but the grandfather idenity will need (aB)-1.aaB = baB. We thus need to go back to the Markov identity to find Tbw = Tw ± 2i, since Tb = 2 as previously. In our case, we have TbaB = (Ta +2i) - 2i = Ta, which makes sense since 1/1 with Ta = 2 – 2i is also the Apollonian gasket. Hence we get Tw(1+1)/(1+2) = Tw2/3 = Tw1/1.Tw1/2T(baB) = Tw1/1.Tw1/2Ta, giving a value of 1.6939-1.4188i, by comparison with the 1/3 Ta value of 1.6939 - 0.5812i, with imaginaryparts neatly adding to –2i and forming mirror images. However, in the process of root searching, you will also find some unstable values which still have the appearance of legitimate circle packings but will compute only with the stochastic algorithm below, such as Ta= sqrt(2) - 2i, fig 10 which appears to be an unstable point off the Maskit slice.

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-1n, 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.01161324i 13/40 1.68667404 – 0.57810830i 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), μ(8/13), μ(21/34) and μ(89/144) tending to the Golden Mean γ = (–1+sqrt(5))/2. The limit sets tend in limit to the doubly degenerate Troell's point illustrated in Fig 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.0964i, Tb = 1.7667 – 0.2201i and Ta = 1.8385 – 0.0964i, Tb = 1.7667 + 0.2201i. Right: Two limit set with Ta, Tb having values determined by μ(3/5), μ(8/13) and their double conjugates, showing two fractional parabolic traces can also be combined, but here the inner circles have now been lost.

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), μ(8/13) both of which are parabolic since they each form gaskets with trace Tb = 2.


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

The final examples in this section are non-free groups, in which the generators a, b are chosen so that aba-1b-1aba-1b-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-1b-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  you will see setting 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.00025756i, 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-1a-1ba and aba-1b-1a. 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 a10b 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.05i, 1.89+0.12i, Centre Left: 1.94+0.05i, 1.94+0.05i with similar non-free coefficients to the left hand images of fig 2. Centre right:Aa limit set generated from μ(3/5),2 with recursecoeffs() and then coupling this as Tb with μ(1/7) using recurscoeffs2() and nonfree() to get Tab coloured yellow to emphasize the periodic circles in the central region, Right: Limit set with Ta, Tb having values determined by μ(3/5), μ(8/13) showing two fractional parabolic traces can also be combined.

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.

Stochastic Algorithm

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.1i, Tb = 2. Research in the Maskit parametrization has focused on the two transforms
a
(z) = μ + 1/z, b(z) = z + 2. In the standard view this corresponds to the traces Ta = - and Tb = 2. Right: Ta = sqrt(2)-2i, Tb = 2 appears to be a circle packing in an unstable island off the Maskit slice.

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.05i, Tb=3, Ta=1.6, Tb=1.6, and Ta=1.85+0.05i, Tb=1.85-0.2i

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.05i, Tb = 1.93 - 0.05i, Tab = 2.2 - 1.8i      Ta = 1.924781 - 0.047529i, Tb = 2, Tab = 1.9248-1.4617i

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

3D Script   

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