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

Chris King

Download the complete Matlab toolbox, C scripts and Mac app with web resource Toolbox contents page.

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

  
Fig 0: Left: Moving set of images of all the Kleinian limit sets of Ta = –iμ(n/61), Tb = 2 generated with the C-script depthfirst61.c This script is nearly 100 times faster than the extended feature Matlab version and will run on any platform. See explanation below.
Right: The Mac DFS Viewer app in operation, displaying in colour the 26-th iterate of a quasi-Fuchsian limit set. The fractal fragmentation of the iterates has only reached part way out from the centre and the tips are still single coloured.

Contents

  1. Overview of the Computational Package
  2. Diverse Mathematical Foundations
  3. The Depth First Search Algorithm
  4. Rational Maps and Period Counting
  5. Finding the Elusive Traces
  6. The Maskit Slice as a Fractal Solution Space
  7. Endless Varieties among Limits of Limits
  8. Lost in Four Dimensional Space
  9. Finding Conjugate Limit Sets
  10. Rational Periodicity becomes Irrational Degeneracy
  11. The Total Eclipse of Double Degeneracy
  12. The Spiral Nemesis
  13. The Holy Grail at the Edge of Chaos
  14. Discreteness Enters the Chaotic Milieu
  15. The Non-Free Case and Automatic Group Tables
  16. The Stochastic Algorithm
  17. Escape-time and Extra Dimensions
  18. Program Listings

Overview of the Computational Package

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

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

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

Ta and Tb set the real and imaginary parts of the two matrix traces, levmax sets an upper limit in the depth search tree, epsilon gives the fineness of the limit set resolution relative to the overall scale used to expand or shrink the image. The timer interrupts the process after a fixed time in seconds the next time it hits levmax, to exit the process if the trace values are super-critical. If only a part of the limit set draws, you need to set the timer for a longer period. Set timer to 0 if you want a high resolution image you know has good traces. Non-free chooses the non-free case discussed below. Necks performs the additional midpoint checks to deal with large narrow-necked fronds as in fig 4, Escape turns on an escape routine if the process begins to go unstable or chaotic. If super-critical values are input the app will then attempt to detect runaway and exit displaying the chaotic iteration up to this point. but can also disrupt depiction of large limit sets. Iterate, if set to a non-zero value n, colours each pixel by the n-th iterate of the four generators – a, b, a-1, b-1. Set to 0 it gives a black and white image, and set to -1 it gives the final level iterate for each pixel. Alt chooses the alternative square root for Tab, which sometimes differs between the Matlab and C versions. If you find you are getting different effectsfor a given set of trace parameters, try choosing alt root, and/or reversing the signs of the real or imaginary parts of the traces.

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

On a Mac, open depthfirst.c using Text Edit and edit it to change the parameters to suit your example. Then open Terminal and navigate to the folder containing depthfirst.c. You can do this by typing cd< folder name> until you get to the folder, or simply by using cd<space> and dragging the folder icon rfrom the Finder into the Terminal window. Then copy and paste the gcc command at the top of the C script into Terminal and the script will compile and run. After this make a copy of the .ppm file, call it "test.ppm" and copy and paste the sips command into Terminal and voila! – you have a full resolution highly compact colour image – test.png. Alternatively a Matlab user can execute replayC() on the original filename without .ppm. to get te same png image. If you don't have Apple's XCode installed you will need to download it, free from the Apple Store and then install the Terminal command line tools by executing xcode-select --install in Terminal to activate gcc. On older systems, a suitable version of XCode can be downloaded by signing on for a free Apple Developer account. The script can also be edited, compiled and run in XCode as a command line tool.

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

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

An Intersection of Diverse Mathematical Areas

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

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


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 (right) to 2 (left), 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 mapped over the white circles, forming a modular subgroup of the gasket group. The limit set (centre ) with traces Ta = 1.92 + 0.03i, Tb = 2 is colour coded by the first 29 successive iterates. The right-hand figure shows locations of some key elementary repeated elements.

Mobius transformations, as noted, correspond to complex versions of the rotations, scalings and translations of affine Euclidian maps.They can be decomposed into similar elementary transforms and are classified into elliptic, parabolic, hyperbolic and loxodromic types. Parabolic transforms are like transalations and have only one fixed point on the Riemann sphere, like the point at infinity of a translation. Loxodromic maps are complex dilations. They have a repelling and an attracting fixed point and the map expands in a spiral form from the repeller to the attractor. Elliptic maps are conjgugate to pure rotations. Pairs of Mobius transformation induce a group structure of transformations under composition. The limt sets consist of the asymptotic points the rest of the Riemann sphere is eventually mapped into under the group transformations. Conjugacy powerfully simplifies the classification, because any map M = C-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 of the limit set algorithm, were the exact action of each generator on the limit set becomes apparent, as in fig 1(right) and fig 4. The traces of parabolic maps, defined by p + s above (the sum of the diagonal elements) lie on the real line between –2 and 2. Those of elliptic maps lie on the unit circle excepting 1 and –1. Traces of hyperbolic maps lie on the real line outside [–2,2] and loxodromic maps lie anywhere in the complex plane outside these ranges, although hyperbolic are sometimes included in loxodromic.


Fig 2: Diverse mathematics underlies the forms of the limit sets. (b) shows Schottky circles of non-kissing circles as in (2), defined by loxodromic maps on the Riemann sphere (a), taking the inside of one circle over the outside of the paired one, resulting in fractal dust Cantor sets as the limit sets. (c) shows how this becomes a circle when pairs of circes touch. (e) shows the Riemann sphere with the two circle interiors removed ready for gluing in ways which can generate new surfaces. (h) is the construction for the modular group represented by the Escher tiling consisting of the group of matrices with determinant 1. (f) shows the arrangement for a quasi-Fuchsian group, where the touching circles are further drawn together, so that the surface becomes two punctured tori joined at their cusps and the limit 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 causing the pinching. As we move the circles, more can kiss, so that when all do, as in the gasket (d), two more transformations have become 'pinched' so that both a and b are parabolic and we end up with (g). The four tangent circles generating the Kleinian group of the Apollonian gasket partially portrayed in (d) are highlighted in black. One of the four circles is the line y = 0, with centre the point at infinity. The gasket group is called doubly cusped because we have pinched two extra handles because, a and b are also parabolic. It is also sometimes called maximally cusped, because, after all this squeezing, there are no more curves left to pinch. In many of the cases below with Tb = 2, only one generator b is parabolic, although in double cusp limit sets both Tb and a generator word defining the fractional motion are, and in some such as fig 3, neither are, although each quasi-Fuchsian limit set is still a fractal topological circle arising from the parabolic commutator..

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


Fig 3: Left and centre: Quasi-Fuchsian limit sets Ta = 1.91 ± 0.05i, Tb = 1.91 ± 0.05i  If we independently exchange the signs of the real and imaginary parts to find complementary points on the Maskit slice, which is symmetric under changes of sign of both x and y, we will find four symmetry variants of each of these sets.

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

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

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

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

The DFS algorithm has an ingenious strategy designed to work with transformations with parabolic commutators that produce quasi-Fuchsian groups which can be scanned in a circular manner by exploring the tree of words defined by compositions of the transformations. Firstly "grandma's algorithm" is used to generate Mobius transforms with parabolic commutators (aba-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 corkscrew-wise, making a right-hand turn after each level descent, and then exiting the level when we reach the inverse of the transform we entered the level by after three left turns to explore the other tunnels. It makes a branch terminal check by applying the composed transformation to the fixed point of the appropriate commutator defined by the last transform e.g. t(1)*t(2)*...*t(n-1)*a operates on the attracting fixed point of 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 4: 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. Secondlythis process is also invoked to avoid pinching off the whorls of limit sets with large appendages with narrow necks as in fig 4. To do this, in addition to the attracting fixed points of the four sets of commutators we also include those of the original four generators, which, as can be seen from the first image in fig 1 and the image on the right from Indra's pearls, are at the very tips of the whorls. Fig 4 illustrates the problem for the second example in fig 3, where two "new" and "old" candidates at the neck are ε-close, but the "tip" of the whorl is further away and joining them would cut the whorl at the neck. This means we can then decide to fill in a line only if the distances from our "new" and "old" points to the local whorl tip we call "mid" are both less than ε. In the case of the gasket we fill "new" to "mid" to "old" as the 'tip' is right in the centre of the gap, but for other sets just "new" to "old" provided both ε's are small enough.

Further examples aided by special words to get vey high reolution images efficiently arise when dealing with the cusp examples below, where other elements are chosen to be parabolic, such as 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 and the algorithm is already giving high reolution images of these in real time.

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


Fig 5: Lacy quasi-Fuchsian limit sets where Tb = 3 Left: Ta = 1.91 + 0.05i , Tb = 3 IPs 8.9 238   Right: Ta = 1.9656+0.99536i, Tb = 3 crop of spirals IPs 8.23 264     

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

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

Rational Maps and Period Counting

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

In fig 7 you can see examples of rational cusps 1/15, 2/19 and 7/43, with the central chain of circles numbered to highlight the relationship with the denominator. In the left figure the fixed point of the parabolic element 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 meaning there is one fixed point (parabolic) on the tangent point. b is parabolic with fixed point -i and pushes points out from its fixed point along clockwise circular trajectories. a is loxodromic and pushes the disks in the numbered spiral chain from right to left across the shrinking and spiralling into the sink on the left.The chain of 15 circles connecting to the opposite large circle 16 map to one another and the mauve circles are in turn emerging from the action of loxodromic a in a spiral from the repeller and proceeding in the orange numbered spiral to the attractor. In the 7/43 case, you can also see the way the chain is running in cycles of seven, corresponding to the numerator, and cycles of 2 in the 2/19 case.

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

More accurate plots of the small circles can be gained by including the parabolic elements such as 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 ε.


Fig 7: Double cusp 1/15 Ta = 1.958591030 + 0.011278560i, Tb = 2   both b and a15B are parabolic IPs 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  IPs 9.15 293  a10Ba9B and b are parabolic.  Right: 7/43 cusp  Ta=1.80785524 + 0.136998688i, Tb = 2  IPs 9.16 294 − 295 plotted to higher resolution using the associated C-program.   You can see the denominators clearly in the periodic mappings of the circles, and the numeratos both in the stepping periods along the chain and in the chains between the horizontal circle pair and the vertical pair.  

Finding the Elusive Traces

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

The traces combine in two key relations: (1) The "grandfather" identity Tmn + Tm-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.

Key to this is the nature of any composition of the generators which is parabolic. As an exercise in finding the trace coefficients of some of these examples lets look at the cases 1/n where the parabolic element is 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 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 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. Theory confirms that among all possible solutions to the trace equation (for a fixed p/q), exactly one will be discrete, so only one of these will have a credible root one can try. You can directly input the values such as s(7) into depthfirst as parameter, having run this function, giving full precision. Example values are 1/7 1.846276 – 0.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 11 centre.

Going beyond the 1/n fractions to p/q, we can also calculate the trace Ta for general fractions, because parabolic expressions like a3Ba3Ba2Ba3Ba2B can also be generated recursively e.g. 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).

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

A list of some fractional traces.

Coefficient recursion  tracepq() program listing in colour


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

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

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

2    -1.53746746-0.21870978i

5    -1.40222578-0.58650847i

4    -0.39303611-0.25387457i

3    -0.38479895-1.27170824i

1    0.79526188-1.15961696i

6    1.25146576-0.06152298i

7    1.67080064-0.44805900i

We can then follow with a portrait using depthfirst()

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

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

The Maskit Slice as a Fractal Solution Space

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

   
Fig 9: In the centre is shown a section of the "Maskit slice" - the plane of μ(p/q) where the generating p/q word, such as a15B is parabolic from Indra's Pearls (the Matlab version in in fig 10 below). A Matlab generated version is shown in fig 10. In the Maskit parametrization we have one parameter μ defining a Mobius transformation which is paired with translation by 2. In the standard parametrization we get the same limit set from Ta = -i μ(p/q), Tb = 2. 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 coloured by the final iterate. Notice the overlapping fractal dark circles in the right hand image. For an explanation of μ in relation to the trace values, see fig 23.

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


Fig 10: Left: Ta= –iμ(i /30) Tb = 2, i = 0-29 showing all the rotation states mod 30. This can be contrasted both with the version for prime 29 (800x800 moving gif for n/29) and fig 0 where prime 61 is shown. Center: "Maskit slice" as in fig 9 centre on a rectangle between -2 and 2 by generating the symmetrized plot of all roots up to q = 51 coloured from 1 red to 51 blue using tracepq() in the script slicec.m and plotting with pltmus(). Right: Conjugate roots for (i /29) and their superposition. The sequences for n/29 and n/61, where all the fractions are irreducible because they are prime can be contrasted with the 800x800 moving gif for n/30, where there are factors of 2, 3 and 5 in the iterating steps. Watch n/29 and n/30 move together.

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

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

Endless Varieties of Limits of Limits


Fig 11: Limit sets of μ(1/7), μ(1/15), μ(1/31) and μ(1/60) show they are not tending to the limit of μ(1/n) = μ(0/1), the Apollonian gasket, as 1/n tends to 1/∞ = 0/1. The limit set of μ(1/31) is pictured in blue with the gasket superimposed in black to highlight the relationship. The Apollonian gasket has trace Ta = 2 which is the limit of the traces of these limit sets and its matrix a is also the limit of the μ(1/n) 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 black on the limit set of μ(1/31) in blue (click to enlarge).

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

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

Hence we can see that (anB)μ(1/n) introduces in limit a new symmtery involving rotation by π / 2, not existant in the gasket group.

Fig 12: μ(2/(2n+1)) for n = 60, or μ(2/121), demonstrating a third limit of limit sets tending to μ(0/1) different from both the above sequence in fig 11 exemplified in limit by μ(1/60) and the Apollonian gasket μ(0/1). Note there are now pairs of circles separating the principal ones rather than one as in μ(1/60). This is the second in an infinite sequence of distinct limits of limits.

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, as shown left for μ(2/121). And we could do this for 3/(3n+1) and 3/(3n+2) and so on so we have an infinite number of different limits of limit sets, all tending to μ(0/1).

Lost in Four Dimensional Space

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


Fig 13: Kleinian limit sets where both traces are derived from parabolic examples where TB = 2. Left: Ta = μ(3/29), Tb= μ(10/29), Right: Ta = μ(7/43 Tb = spirlally degenerate.

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

   

Fig 13a: 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.

Finding Conjugate Limit Sets

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

aaB(aaBaB)-1aaB = aaB(bAbAA)aaB = aaBbAbAAaaB = a
(aaB)-1aaBaB(aaB)-1(aaB)-1aaBaB = (bAA)(aaBaB)(bAA)(bAA)(aaBaB) = B.

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


Fig 13b: Left: μ(8/13),2. Centre: (w1/2,w2/3), which is clearly conjugate to the left figure, as confirmed by the alternate root of the same traces illustrated right, giving an identical limit set to the left one except horizontally flipped, and rotated through 45 degrees.

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

Rational Periodicity becomes Irrational Degeneracy


Fig 14: Periodicity limiting in degeneracy: Fibonacci fractions μ(3/5), μ(8/13), μ(21/34), μ(55/89) and μ(144/233) tending to the Golden Mean γ = (–1+sqrt(5))/2. The limit sets tend in limit to the degenerate golden mean Troels' point illustrated in Fig 16.

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


Fig 15: Singly degenerate Troels' point Jorgensen's group Ta = 1.6168866453 − 1.29432650i, Tb = 2 coloured by first iterate to highlight the boundary lines between the complementary discs left and right being consumed by the degeneracy IPs 10.4 318-319. This value corresponds to the minimum imaginary value on the boundary of the Maskit slice fig 11 centre and corresponds closely to the μ value of the golden mean. It is the mirror image of Ta = 1.6168866453 − 0.7056734968i. See fig 14 for sequence convergence images to the golden mean. Right singly degenerate  Ta = 1.7847+0.73852i, Tb=3 Ips 10.1 311 

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

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

(Ta, TaB, Ta.TaBTB)   <L   (Ta, TB, TaB)    – R–>     (TaB, TB, TBTaBTa)

and these compose, so that for RL we get (Ta, TB, TaB)RL –> (Ta.TaBTB, TaB, T 2aB.TaTaB.TB – Ta). Hence these relations can also be used to decode other systems like L10 R. In these progressions, including the one above, the geometric scalings become determined by the eigenvalues of the transformation matrix. This leads to the discovery that the boundary of the Maskit slice is asymptotically self similar. You can see the nature of this phenomenon in fig 17 lower left inset, where successive sections of the corresponding spiral degeneracy are fractally scaled and rotated versions of the original. In our case above, one of the eigenvalues turns out to be (5 + 211/2)/2 ~ 4.791287848, so the traces of each successive pair of Fibonacci fractions are about five times closer. For example

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

The Total Eclipse of Double Degeneracy

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


Fig 16: Left: A singly degenerate conjugate of the Troels' point map above, forming an approximation to the doubly degenerate set right. As in fig 13b, the alternate grandma's root gives a circular limit set, again conjugate to the original even though Ta and Tb are not exact conjugates, and despite the mid-resolution iterates of the DFS algorithm (inset) proceding on an entirely different basis. Right: Doubly degenerate limit set LR coloured by the third iterate drawn on the same scale as left of the space-filling in the plane Ta,Tb = (3±31/2i)/2 = 1.5 ± 0.8660254038i IPs10.11 335. The fixed point o the commutator lies at the cusp between one of the principal centre pairs and the additional transformation, c, which is also parabolic, has the same fixed point and commutes with the commutator, although distinct from it.

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

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

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

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

for i=1:length(s)

    if s(i)==1

        p0=p0+p1;

        q0=q0+q1;

    else

        p1=p1+p0;

        q1=q1+q0;

    end

end

To actually get a doubly degenerate group we aim to find a seqence of generator fractions pn / qn and –qn / pn so that in the irrational limit r and –1/r , they are degenerate both on the inside and outside. For this relationship to happen, the words get extremely long and there is little variation near the center, which suggests the traces of the original (a, B) should be invariant in the sense that they are the same as those of the next in the list (aaB, aB), since the sequence LR takes us as follows (a,B) > L > (a,aB) > R > (aaB,aB). This system has three defining equations: TaaB = Ta.TaB TB = Ta from the grandfather identity, TB = TaB, plus the Markov equation T 2a + T 2B + T 2aB = Ta.TB.TaB, which provides three equations to solve for Ta, TB and TaB. Hence T 2a + T 2B + T 2B = Ta.TB.TB, or T 2a = (Ta.– 2)T 2B and Ta.TaB TB = Ta, or (Ta.– 1)TB = Ta. Solving for Ta gives T 2a – 3Ta + 3 = 0, or Ta=(3± 31/2i)/2, TB =Ta / (Ta-1) = (3-/+ 31/2i)/2, with conjugate trace values of 1.5 ± 0.866025i closely neighbouring the above values in the previous paragraph. We have now induced a doubly degenerate limit set LRengulfing the entire complex plane. The limit set displayed in fig 16 right clearly has convergence to degeneracy on both sides of the remaining white islands, confirming its double nature.

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

The Spiral Nemesis


Fig 17: Left: Free discrete group Ta = 1.8866–0.05i, Tb = 2, IPs 8.15 246 Centre and right: Spirally degenerate Ta = 1.9264340533 – 0.0273817919i, Tb = 2, IPs 10.6 322-323.The centre image shows partial engulfment and the right anoter stage further. It took 9.7 hours to run in the C script using a levmax of 10000 and an epsilon of 0.0000025 and would have taken 36.9 days in Matlab because of the very slow engulfing convergence for this limit set. One can see that this limit set is pervaded by two different types of double spirals, the main set consisting of the main central, four peripheral and so on as in the centre image and all other similar examples, such as figure 11, and a second set unique to this system, arising from the spiral degeneracy. Inset below left: Self-similarity of the Maskit slice on approaching the values leading to sprial degeneracy. In the asymptotic limit the scaling factor approaches 86.214 + 164.198i, accounting for both the scaling and rotation.

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

The Holy Grail at the Edge of Chaos

We can also use the above techniques to define doubly spirally-degenerate limit sets based on the same development of double degeneracy applied to the spirally-degenerate example above. To see how this works, let's consider instead of alternating LRLR..., as in the doublydegenerate case above, the sequence LLRLLR..., which is a step to the L10 R in the above spiral case. The sequence LLR takes us as follows (a,B) > L > (a,aB) > L > (a,aaB) > R > (aaaB,aaB), so we need TaaaB = Ta.TaaB TaB = Ta, TaaB = Ta.TaB TB = TB and T 2a + T 2B + T 2aB = Ta.TB.TaB. This system solves out via (1) Tab = Ta(TB – 1) and (2) TB =T 2a/(T 2a – 2), to give T 4a – 5T 2a + 8 = 0 or Ta = ±1.63224188 ± 0.40523272i. (2) then gives TB = ±1.5 -/+ 1.32287565Ii, giving non-conjugate traces, which grandma confirms satisfy TaaaB = Ta and TaaB = TB, producing the limit set, fig 17a left, which has similar characteristics of double degeneracy to the previous example in fig 16 which can be compared with the limit set of (LLR)5 = μ(153/418).


Fig 17a: Left μ(153/418),2 - LLR5 compared with the doubly degenrate LLR Right: μ(91/345),2 - LLLR4 compared with the doubly degenrate LLLR. Notice in the expanded voews, both the inner and outer whorls of the doubly degenerate cases reflect the limiting behavior approaching the horizontal 'circles' in the singly degenerate examples.

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


Fig 17b: Left: The corresponding motifs in single (L4R)3 35/169 and double L4R clearly show their relationship.
Right: Double degenerate L6R and single (L6R)3 63/433 generated using tracepqr() from the trace of 8/55 as seed.

For L4R , we have TB = (T 4a – 2T2a) / (T 4a – 3T2a) and TaB = Ta(TB-1) = T 3a / (T 4a – 3T2a), giving T 10a - 6T 8a + 8T6a + 4T 4a = 0 andTa = 1.79210587 - 0.141970345i, Tb = 1.64779887126 + 1.721433237i, which again compares with (L4R)3 35/169 Ta = 1.75300569-0.22184750i as shown above. This leads to a systematic way to resolve any Ln R .


Fig 17c: The Holy Grail: Left: Doublespiral degeneracy – genuine L10R – deduced as outlined below, compared centre with the previous image of spiral degeneracy.On finer resolutionwith a higher levmax, both the external and internal fronds to either side erupt further in space-filling spirals as shown in the moving gif, but the convergence is very slow. The last image took 28.9 hours running on C and would have taken 112 days on Matlab. Right: The limit set with Ta and TB swapped to show it has equivalent fractal structure, when the 'outside' is exchanged with the 'inside'.

One can see from the repeating grandfather identity how this process extends readily to the doubly spirally-degenerate solution of L10R. The corresponding system of equations iteratively solves TanB back to TaB and TB, giving three equations in Ta, TB, and TaB. Because the first two equations are linear in TB and TaB, one can eliminate these, to give TB and TaB as above, to get a higher degree polynomial in Ta.

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

syms a B aB

a2B=a*aB-B;

a3B=a*a2B-aB;

a10B=a*a9B-a8B;

a11B=a*a10B-a9B;

 

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

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

 

a11B(a, B, aB) = aB*a^10 - B*a^9 - 9*aB*a^8 + 8*B*a^7 + 28*aB*a^6 - 21*B*a^5 - 35*aB*a^4 + 20*B*a^3 + 15*aB*a^2 - 5*B*a - aB

a10B(a, B, aB) = aB*a^9 - B*a^8 - 8*aB*a^7 + 7*B*a^6 + 21*aB*a^5 - 15*B*a^4 - 20*aB*a^3 + 10*B*a^2 + 5*aB*a – B

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

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

o = expand( (B*a^9 - 8*B*a^7 + 21*B*a^5 - 20*B*a^3 + 5*B*a + a)*(a^9 - 8*a^7 + 21*a^5 - 20*a^3 + 5*a )
               
- (a^10 - 9*a^8 + 28*a^6 - 35*a^4 + 15*a^2 - 1 )*B*( a^8 - 7*a^6 + 15*a^4 - 10*a^2 + 2))
o = 2*B - 15*B*a^2 + 35*B*a^4 - 28*B*a^6 + 9*B*a^8 - B*a^10 + 5*a^2 - 20*a^4 + 21*a^6 - 8*a^8 + a^10

Hence:

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

syms a B aB

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

t0=a;

t1=aB;

t2=B;

for i=1:length(s)

    if s(i)==0

       h=expand(t1*t0-t2);

       t2=t1;

       t1=h;      

    else

       h=expand(t1*t2-t0);

       t0=t1;

       t1=h;      

    end

end

aB1=solve(t0-a,aB);

aB2=solve(t2-B,aB);

B=solve(aB1-aB2,B);

aB=a*(B-1);

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

g=expand(f);

h=simplify(g);

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

 

h = (a^2*(a^12 - 12*a^10 + 54*a^8 - 112*a^6 + 105*a^4 - 37*a^2 + 8))/(a^6 - 6*a^4 + 9*a^2 - 2)^2

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

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

Discreteness Enters the Chaotic Milieu

  
Fig 17d: The output of the DFS algorithm in the Mac app when the real part of Ta in the example of fig 0 right is changed to 1.41891069, 1.51891069 and 1.71891069 respectively.

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

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

The Non-Free Case and Automatic Group Tables


Fig 18: Left: Sequence of non-free group limit sets for fractions n/29. See both the n/29 rotations together. Right Non-free Maskit slice.

The final examples in this section are non-free groups, in which the generators a, b are chosen so that C2 = (aba-1b-1)2= I and hence these relations between the generators become trivial identities, 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. C2 = I implies TC = 0 since the off-diagonal elements of C2 i.e. (p+s)q = (p+s)r = 0, the off-diagonal elements of I, so TC = p+s = 0 unless r = q = 0. 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 tracepq() 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 19: 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

We then have to set up an automaton(right) in the form of a table generating the dead ends in every part of the maze by recursively looking up successive words and terminating at neutralizing words such as b-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 x counting as a death state. The table is used sequentially to take us from any state on the left by right multiplication by one of the four generators. The states are stored sequentially as we move down a level of the maze so they can be looked up and the next recalculated. The end result, centre, has two disconnected sets of circles, coloured white and yellow in the figure. It is also a double cusp in which b and 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 version of the 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

 


Fig 20: 1.906-0.06739i, 2 Near degenerate high period, 1.865-0.70567i, 2 Close to 1/3 of a revolution     μ(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.   Spirally degenerate L10R 143/1561 1.934141199416-0.036681260008i     

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

>> s=tracepq(3,7,2,1);
3 -1.65805636-0.71965515i
1 -1.61195573-0.46319322i
4 -0.42695172-0.16466462i
5 -0.25623249-1.15299930i
6 1.03890744-0.78635572i
7 1.08561081-0.34873804i
2 1.82867804-0.60703463i

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

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


Fig 21: 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 4. Centre: A limit set generated from μ(3/5),2 with tracepqs() and then coupling this as Tb with μ(1/7) using tacepqs() and nonfree() to get Tab coloured yellow to emphasize the periodic circles in the central region, Centre right: Limit set with Ta, Tb having values determined by μ(3/5), μ(8/13) showing two fractional parabolic traces can also be combined. Right: 1.91 – 0.1i, 3

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


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

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

The Stochastic Algorithm

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

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

Fig 23: 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, so is an example of a non-free group.

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

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


Fig 24: A variety of non-free limit sets induced by choosing traces within the Maskit slices. The coefficients are as shown with 1/n corresponding to the non-free traces shown below. The non-free example has Ta=21/2, Tb=2
with command qFstochastic(sqrt(2),2,10,10,1100,0.4,2,2,nonfree(sqrt(2),2,0,0));. Intriguingly the real solutions to 1/n for the non-free case T(aba-1b-1) = 0 and the standard case T(aba-1b-1) = -2 are the same, although all the other complex traces calculated by tracepq() differ.

There is an inriguing pattern in the non-free trace values giving an = I. For each trace solution μ(1/n) there turn out to be ceil(n/2)-1 real solutions, Ta, each of which is a trace solution for an = I, the most negative of which gives a discrete non-free limit set. The even n values have roots coming in positive and negative pairs reflecting n/2 as well, as one of the pair has an/2 = -I so an = I and the other an/2 = I. The solutions for n rapidly become intractable to solve explicitly. For n = 7, for example, the solutions factorize into a quartic and a cubic, x3 + x2 - 2x - 1 having entirely real roots, with one being as follows:

-1.8019= - (2^(2/3)*7^(1/3)*(1 + 3^(1/2)*3i)^(1/3))/12 - (2^(2/3)*7^(1/3)*(1 + 3^(1/2)*1i)*(1 - 3^(1/2)*3i)^(1/3))/12 + (2^(2/3)*3^(1/2)*7^(1/3)*(1 + 3^(1/2)*3i)^(1/3)*1i)/12 - 1/3.

It is thus intriguing why the μ(1/n) generated by the simpler recursive application of the Grandfather identity to the 1/n series always has such a solution for each n, while all the roots of μ(m/n) m >1 are complex and hence loxodromic Mobius transformations. As shown below, the factoring of the trace polynomial of Ta clearly splits recursively into a product of 2 (or more) polynomials, in which one or more of the lower dimensional ones have only real roots Ta which satisfy an = I.

factor(a{1}-2) = 0
n = 1   factor(a{2}-2)   = t - 2 + 2i
n = 2   factor(a{3}-2)   = t^2 + t*2i - 4
n = 3   factor(a{4}-2)   = [ t + 1, t^2 - t*(1 - 2i) - (2 + 2i)]
n = 4   factor(a{5}-2)   = [ t, t^3 + t^2*2i - 4*t - 4i]
n = 5   factor(a{6}-2)   = [ t^2 + t - 1, t^3 - t^2*(1 - 2i) - t*(3 + 2i) + (2 - 2i)]
n = 6   factor(a{7}-2)   = [ t + 1, t - 1, t^4 + t^3*2i - 5*t^2 - t*6i + 4]
n = 7   factor(a{8}-2)   = [ t^3 + t^2 - 2*t - 1, t^4 - t^3*(1 - 2i) - t^2*(4 + 2i) + t*(3 - 4i) + (2 + 2i)]
n = 8   factor(a{9}-2)   = [ t, t^2 - 2, t^5 + t^4*2i - 6*t^3 - t^2*8i + 8*t + 4i]
n = 9   factor(a{10}-2) = [ t + 1, t^3 - 3*t + 1, t^5 - t^4*(1 - 2i) - t^3*(5 + 2i) + t^2*(4 - 6i) + t*(5 + 4i) - (2 - 2i)]
n = 10 factor(a{11}-2) = [ t^2 - t - 1, t^2 + t - 1, t^6 + t^5*2i - 7*t^4 - t^3*10i + 13*t^2 + t*10i - 4]
n = 11 factor(a{12}-2) = [ t^5 + t^4 - 4*t^3 - 3*t^2 + 3*t + 1, t^6 - t^5*(1 - 2i) - t^4*(6 + 2i) + t^3*(5 - 8i) + t^2*(9 + 6i) - t*(5 - 6i) - (2 + 2i)]
n = 12 factor(a{13}-2) = [ t, t + 1, t - 1, t^2 - 3, t^7 + t^6*2i - 8*t^5 - t^4*12i + 19*t^3 + t^2*18i - 12*t - 4i]
n = 13 factor(a{14}-2) = [ t^6 + t^5 - 5*t^4 - 4*t^3 + 6*t^2 + 3*t - 1, t^7 - t^6*(1 - 2i) - t^5*(7 + 2i) + t^4*(6 - 10i) + t^3*(14 + 8i) - t^2*(9 - 12i) - t*(7 + 6i) + (2 - 2i)]

In the case of n = 5 with x2 + x - 1, the discriminant is 5 > 0 and the cubic of n = 7 with x3 + x2 - 2x - 1 likewise has disriminant 49 > 0, guaranteeing these both have entirely real roots, all of which satisfy a5 = I and a7 = I respectively. The key to this pattern is that the Trace equations above are generated in each case as solutions for Ta satisfying TanB = 2. But for this 1/n word, unlike the general p/n case, there is clearly the possible solution an = I, since then TanB = TB = 2. Such solutions, have to be real because complex traces lead to loxodromic Mobius transformations, for which no power of a can become the identity. Notice the real roots of non-prime n include the roots of their factors, expressed as sub-factors.

Attempting to solve the trace equations for an = I directly leads to a partial overlap with the above polynomials which generates more potential solutions than allowed under the constraints. For example a3 = I leads via the Grandfather to Ta2a + Ta2a-1= Ta2Ta, and Ta2+ Taa-1= TaTa, so Ta2= ( Ta)2 -1 since Taa-1= TI = 1. Hence we have Ta3= ((Ta)2 -1)Ta - Ta = 1, giving (Ta+1)((Ta)2 - Ta - 1)= 0. This throws out the correct solution Ta = -1 plus 2 real roots which are violated by other conditions.

Similarly a2 = I since Ta2 = (Ta)2 - 1 = 1, leads to the false solutions ±21/2, which nevertheless appear later as a correct solution to a8 = I for n = 8. On the other hand, we know from the previous section that a2 = I can imply Ta = 0 since the off-diagonal elements of a2 i.e. (p+s)q = (p+s)r = 0 since the off-diagonal elements of I are 0, so Ta = p+s = 0 unless r =q = 0. Again the trace 0 input to Grandma's recipe gives a2 = -I or a4 = I and this solution is also illegitimate for n = 2, although it also reappears for n = 4. In fact n = 2 has no real solutions as listed above.

These examples also show us that looking for limit sets having higher power non-free commuator relations such as (aba-1b-1)3= I would give a commutator trace e.g. of -1, as noted for a above, and for all n would range between 0 and -2, generating intermediate limit set series, as illustrated for Ta = Tb = 2 in fig 26.

The valid real roots form a sequence of traces tending to 2, with even numbers giving factored values as follows (n, Tn) = (3,-1), (8, -21/2 = -1.414), (5, -(1+51/2)/2) = -1.618, (12, -1.732), (7, -1.802), (16, -1.847), (9, -1.879) (20, -1.902), (11,-1.919 ),(24, -1.932), (13, -1.942), ..., (61, -1.99734), with the number of central circles between the top and bottom increasing by 1 at each stage. However the values of n remain faithful to their non-free relation, so while Ta = T(1/16) lies between Tb = T(1/7) and Tc = T(1/9), and b7 = I and c9 = I, a8 = -I so we still get a16 = I. Where the positive root (opt=4) is chosen for the negative trace values or we can alternatively use the ususal negative root (opt=0) with the minus of the real value. As noted in fig 24, these all appear to give elliptical variants of the limit sets, where a parabolic trace of 2 is replaced by one of these, and the same values intriguingly apply to both the standard case where the commutator has trace -2 and the non-free case where it has trace 0.

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

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

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 26: 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.99 showing the evolution of the limit sets towards non-discrete chaotic instability. There is a transition into chaos as we cross 1. The last value 1.99 approaches 2 where the extended Grandma recipe becomes singular and returns 0/0 Nans or infinities for a, b. The scale also fluctuates extensively as we move along the sequence. Two neighbouring complex trace values are also shown to illustrate the dynamics when the commutator is loxodromic. One should note however that the traces do not form a spectrum as such, but rather traces of -1, 0 and 1 depend on powers of the commutator which equal the identity ( 3, 2 and 6 respectively). The right-hand example in fig 26a shows that discrete limit sets can exist in the interval (1,2).

The limit set below with a commutator such that abAB14 = I, and b3 = I shows both close relationship with the hyperbolic plane and inversion of the exterior space of the limit set in the centre of every whorl. The corresponding system with positive trace and abAB7 = I is chaotically unstable in the stochastic algorithm but is nicely depicted by a variant fo the depth first search in help tutorial 4 in the MAF package introduction (see below). Alun Williams, who developed MAF notes in this case that the algorithm struggles, because the usual generators of the group don't have a nice automatic structure. As lesson 4 discusses, changing to using three elliptic generators [a, b and x] gives a much more elegant structure, which is easy using MAF, which he originally started developing to enable him to generate the required automatic structures for such groups 'on the fly'. In MAF, this process is spelled out clearly as noted below, although subsequent steps are quite tricky, and further complexities ensue, involving extremely long word lists and indirect inferences, to arrive at a tractable automatic system.

"In order to find a correct presentation it is necessary to know that we can find a matrix x with trace 0 and determinant 1 such that xx=-I and x*a*x=-a^-1 and x*b*x=-b^-1. (One can find a suitable matrix x by computing the Lie bracket (a*b-b*a) of a and b and normalising the resulting matrix to have determinant 1). This means that the matrix a*b*a^-1*b^-1 = a*b*x*a*x*x*b*x = -(a*b*x)^2. The reason our original input file does not work is that it does not capture the fact that a*b*x should have order 7, rather than 14. It is perhaps not obvious that this should be the case, but the particular value chosen for the trace of a*b*a^-1*b^-1 guarantees that it is". Although x has order 4 as a matrix, as a Moebius transformation it is an involution. The group we want is the subgroup of this group that is generated by a and b, so we can find a presentation for the group by processing this input file as a coset system..


Fig 26a: Left the limit set of Ta = 1.91+0.5i, Tb = 1 with TabAB = -2cos(2π/7) using the stochastic algorihm. Centre: A blow up showing the exterior domain of the set is inverted inside the centre of each whorl. Right: the limit set of Ta = 1.91+0.5i Tb = 1 with TabAB = +2cos(2π/7) using a variant of the DFS with an amended automaton devised using MAF by Alun Williams, who developed the package. Significantly 2cos(2π/7) = 1.247 is one of the real solutions to tracepq(1,7,0,1).

Life at the Edge of Chaos: As can be seen from fig 26, the edge of non-discrete chaos appears to be at TabAB = 1. Fig 26b thus shows some of the standard examples in this case, portrayed using the DFS algorithm. Having initially generated these by the stochastic algorithm, I realized that TabAB = 1 can have the solution TabAB 6 = I, just as a did previously. This presents another type of generator solution then TabAB = 0 because, here the commutator has r = q = 0, rather than p + s = 0. The traces for increasing powers of abAB thus follow the sequence TabAB = 1 -> -1 -> -2 -> -1 -> 1 -> 2 = TI. Since the extended Grandma recipe turns out to generate a, b satisfying TabAB 6 = I, using the extended Markov, the TabAB = 0 automaton works well enough for this case to depict the limit sets quite well, although convergence is limited at high resolution. As noted aove, the packages KBMAG or MAF can be used generate exact automata for these and other Kleinian groups, all of which have been proved to be automatic, although this isn't the case for groups in general.


Fig 26b: A series of examples for TabAB=1 including periodic 1/2, 1/7, 1/19, close to degenerate 144/233 and a quasi-Fuchsian dragon.

Escape-time algorithm   m-file

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 27: 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 28: 3D Plots of quasi3D

Depthfirst Matlab Listing in Colour

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

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

%depthfirst(2,2,0.001,1000,1,0,1); %spec=1 closing gaps

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

%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,-2,1.9248-1.4617i);

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

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

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

 

%Special words example 8.16 with aaB parabolic.

%the function specialwords can be reprogrammed to handle examples like

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

%depthfirst(2+i,3,0.005,3000,1,0,8);

 

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

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

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

% scal is the scaling factor for depicting the image

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

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

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

% you can either produce a black line plot or a coloured image

% iterate = -2 plot n>0 colour image  by nth iterate, -1 by final level, 0 black

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

% If the positive grandma root is to be chosen add 4 to spec.

% To include special words from the subfunction specialwords() add 8 to spec

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

% siz is the image size default 1100x1100

% if you add a tenth parameter tab you can generate non-free limit sets

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

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

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

if nargin<9

    iterate=-2;

    colfac=1;

end

if iterate==-2

    colfac=0;

else

    colfac=1;

end

 

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

% iterate=1; %Which iterate to colour code by -1 colours by the final level

 

ext=0; %set to 1 to retun if image out of range (large or chaotic)

dotc=0;

outfl=true;

if nargin<8

    siz=1100;

end

if nargin<7

    spec=0;

end

if nargin<6

    sav=0;

end

if nargin<5

    scal=1;

end

ep=ep/scal;

tic;

clf;

hold on

M=256*ones(siz);

if mod(spec,8)>=4

    posrt=1;

else

    posrt=0;

end

if nargin<10

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

else

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

end

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<8

    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

%spec=mod(spec,4);

%stpt=0;

lev=1;

tags(1)=1;

state(1)=1;

word{1}=gens{1};

dotfl=false;

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

    fl=false;

    if spec==2||spec==3

        if state(lev)==0

            fl=true;

        end

    end

    while ~fl

         plt=false;

         switch spec

             case {1,3,5,7}

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

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

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

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

                    plt=true;

                end

             case {0,2,4,6}

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

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

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

                    plt=true;

               end

             case 8

               pts=numfp(tags(lev));

               fl4=true;

               pts=numfp(tags(lev));

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

               for i=2:pts

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

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

                       fl4=false;

                       break;

                   end

               end

               if lev==levmax

                   fl4=true;

               end

               if fl4

                   plt=true;

                   newpoint=points(1);

               end

          end

          if lev==levmax

             if sav>=2

                fprintf('.');

                dotfl=true;

                dotc=dotc+1;

                if dotc>=200

                    dotc=0;

                    fprintf('\n');

                end

             end

             if sav==-1

                outfl=false;

                return;

             end

          end

         if spec<8

             x=scal*abs(newpoint);

             if ext && x>50 %breakout to portray chaotic limit sets

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

                   colormap(colorcube);

                   image(255*M');

                   return

             end

         end

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

             fl=true;

             break;           

         end

         if iterate>0

             tagpos=tags(iterate);

         end

         if iterate==-1

             tagpos=tags(lev);

         end

         if iterate==0 || iterate==-2

             tagpos=0;

         end

         if plt

             switch spec

                  case {1,3,5,7}

                      if A==2 && B==2

                        M=myplot(scal,M,oldpoint,midpoint,siz,tagpos,spec);

                        M=myplot(scal,M,midpoint,newpoint,siz,tagpos,spec);

                      else

                        M=myplot(scal,M,oldpoint,newpoint,siz,tagpos,spec);

                      end

                  case {0,2,4,6}

                        M=myplot(scal,M,oldpoint,newpoint,siz,tagpos,spec);

                  case 8

                        for i=2:pts

                            M=myplot(scal,M,points(i-1),points(i),siz,tagpos,spec);

                        end

              end

              if ~stfl

                 if spec==2 || spec==3 || spec==6 || spec==7

                    if abs(oldpoint-stpoint)<10*ep

                    M=myplot(scal,M,oldpoint,stpoint,siz,tagpos,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 || spec==6 || spec==7

                    state(lev)=tags(lev);

                end

            else

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

                if spec==2 || spec==3 || spec==6 || spec==7

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

                    if state(lev)==0

                        fl=true;

                    end

                end

            end

         end

    end

    fl=true;

    while fl

        lev=lev-1; %go backward

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

            fl=false;

        end

    end

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

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

        if dotfl

            fprintf('\n');

            dotfl=false;

        end

        if sav>=2

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

        end

        toc;

        break; %Finished - exit routine      

    end

    lev=lev+1;

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

    if sav>=2

         if lev<4 %progress clock for long iterations

            if dotfl

                fprintf('\n');

                dotfl=false;

            end

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

         end

    end

    if lev==1

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

        %state(1)=tags(1);

    else

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

        if spec==2 || spec==3 || spec==6 || spec==7

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

%             if state(lev)==0

%                 fl=true;

%             end

        end

    end

end

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

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

if spec==2 || spec==3 || spec==6 || spec==7

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

end

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

k=1;

sc11=11/10;

if ~colfac  

    for i=1:s(1)

        for j=1:s(2)

            if M(i,j)<255

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

                k=k+1;

            end

        end

    end

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

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

    text(-1,1.05,tit);

else

    image(M');

    colormap(colorcube(256));

    sz=size(M');

    axis([0 s(2) 0 s(1)]);

end

if mod(sav,2)

%     if colfac==0

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

%     else

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

%     end

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

end

 

function n=mymod(m)

n=mod(m,4);

if n==0

    n=4;

end

 

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

siz=siz/2;

sizs=10/11*siz;

fl=true;

big=10/sc;

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

    fl=false;

end

switch c

    case 0

        co=241; %black

    case 4

        co=186; %magenta

    case 1

        co=233; %blue

    case 2

        co=208; %red

    case 3

        co=165; %cyan

        end

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

%     fl=false;

% end

if fl

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

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

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

    else

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

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

    end

    for i=1:length(zs)

        if ~isnan(zs(i))

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

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

            end

        end

    end

end

       

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

%D=struct();

a=T(1,1);

b=T(1,2);

c=T(2,1);

d=T(2,2);

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

 

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

a=T(1,1);

b=T(1,2);

c=T(2,1);

d=T(2,2);

T2=T^2;

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

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

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

if abs(k)>1

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

 else

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

 end

 

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

if nargin<5

    opt=mod(opt,2);

end

if nargin<4

    opt=0;

end

if nargin<3

    posroot=false;

end

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

    p=-ta*tb;

    q=ta^2+tb^2;

    if posroot

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

    else

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

    end

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

    if opt==0 %Grandma's original recipe

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

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

    a=ab*b^-1;

    end

    if opt==1 %Jorgensen's recipe

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

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

    end

 else

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

  

    Q=sqrt(2-tC);

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

        R=sqrt(tC+2);

    else

        R=-sqrt(tC+2);

    end

    %R=-R;

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

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

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

end

 

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

A=a^-1;

B=b^-1;

numfp=zeros(1,4);

for i=1:4

    switch i

    case {1,3}

        numfp(i)=4;

    case{2,4}

        numfp(i)=3;

    end

end

fix=zeros(4,4);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 

 

 

 

Tracepq Program Listing in Colour

 

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

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

syms t

g=gcd(p,q);

p=p/g;

q=q/g;

s=fareyseq(p,q);

t0=t;

if opt<2

    t1=t+2i;

else

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

end

t2=2;

for i=1:length(s)

    if s(i)==0

       h=expand(t1*t0-t2);

       t2=t1;

       t1=h;      

    else

       h=expand(t1*t2-t0);

       t0=t1;

       t1=h;      

    end

    %fprintf('.');

end

%fprintf('\n');

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

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

if mod(opt,2)==0

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

else

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

end

q=q/gcd(p,q);

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

    for i=1:q

        tt(i,1)=i;

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

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

    end

    tt=sortrows(tt,2);

if prt

    for i=1:q

        if tt(i,3)>=0

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

        else

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

        end

    end

   

end

 

function s=fareyseq(p,q)

p=p/gcd(p,q);

q=q/gcd(p,q);

s=[];

p0=0;

q0=1;

p1=1;

q1=0;

ph=p0+p1;

qh=q0+q1;

while p~=ph || q~=qh

    if ph/qh>p/q

        p1=ph;

        q1=qh;

        ph=ph+p0;

        qh=qh+q0;

        s=[s 0];

    else

        p0=ph;

        q0=qh;

        ph=ph+p1;

        qh=qh+q1;

        s=[s 1];

    end

end

 

Depth First C-script Listing in Colour

 

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

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

//To run successive compiled versions type ./depthf

 

#include <stdio.h>

#include <math.h>

#include <complex.h>

#include <time.h>

 

typedef double _Complex arry[2][2];

 

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

void inv(arry a, arry c);

double _Complex fixp(arry T);

 

time_t start,end;

 

int main(void)

{

    FILE *mf;

    const char base[] = "DF";

    char filename [100];

   

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

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

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

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

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

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

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

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

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

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

    const int levmax=1000,siz=1100;

    double ep=0.001,scal=2;

    double _Complex tab;

   

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

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

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

    int tags[levmax+1];

    int state[levmax+1];

 

    double dif;

   

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

   

    double _Complex fix[3][4];

    int i,j,co;

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

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

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

    int stfl,lev,fl,plt,stps;

    double _Complex oldpoint,stpt,newpoint,midpoint;

    double x;

   

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

    mf=fopen(filename,"w");

    time (&start);

    ep=ep/scal;

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

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

            M[i][j]=255;

    siz2=siz/2;

    sizs=(int)siz2*10.0/11;

    //grandma's recipe

    if (spec==0)

    {

        p=-ta*tb;

        q=ta*ta+tb*tb;

        if (altroot==0)

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

        else

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

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

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

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

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

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

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

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

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

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

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

    }

    else

    {

        p=-ta*tb;

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

        if (altroot==0)

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

        else

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

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

        Q=csqrt(2-tC);

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

            R=csqrt(tC+2);

        else

            R=-csqrt(tC+2);

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

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

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

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

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

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

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

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

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

       

    }

    inv(b,B);

    if (spec==0)

        mult(ab,B,a);

    inv(a,A);

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

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

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

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

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

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

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

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

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

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

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

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

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

    lev=0;

    tags[0]=0;

    state[0]=0;

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

   

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

    {

        fl=0;

        if (spec==2)

            if (state[lev]==-1)

                fl=1;

        while (fl==0)

        {

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

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

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

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

            if (neck>0)

            {

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

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

            }

            plt=0;

            x=scal*cabs(newpoint);

            if (x>10) //Chaotic instability

            {

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

                return(0);

            }

            if (neck>0)

            {

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

                {

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

                    {

                        fl=1;

                        break;

                    }

                    plt=1;

                }

            }

            else

            {

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

                {

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

                    {

                        fl=1;

                        break;

                    }

                    plt=1;

                }

            }

            if (plt==1)

            {

                fl=1;

                z=newpoint-oldpoint;

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

                {

                    stps=1;

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

                    if (x>1) stps=x;

                }

                else

                {

                    stps=1;

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

                    if (x>1)  stps=x;

                }

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

                {

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

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

                    {

                        if (iterate==0)

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

                        if (iterate>0)

                        {

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

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

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

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

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

                        }

                        if (iterate==-1)

                        {

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

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

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

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

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

                        }

                       

                    }

                }

            }

           

            else

                fl=0;

            if (fl==1)

            {

            }

            else

            {

                lev=lev+1; //go forward

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

                if (lev==0)

                {

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

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

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

                }

                else

                {

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

                    {

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

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

                    }

                    mult(h1,h2,h3);

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

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

                    if (spec==2)

                    {

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

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

                    }

                }

            }

        }

        fl=1;

        while (fl==1)

        {

            lev=lev-1; //go backward

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

            {

                fl=0;

            }

        }

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

        {

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

           time (&end);

           dif=difftime (end,start);

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

           break; //finished!!

        }

       

        lev=lev+1;

        if (tags[lev]>0)

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

        else

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

       

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

        {

            time (&end);

            dif=difftime (end,start);

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

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

        }

        if (lev==0)

        {

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

            {

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

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

            }

        }

        else

        {

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

            {

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

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

            }

            mult(h1,h2,h3);

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

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

            if (spec==2)

            {

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

                lev=lev+1-1;

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

            }

        }

       

    }

   

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

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

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

    fclose(mf);

    return(0);

}

 

 

double _Complex fixp(arry T)

{

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

    double _Complex k;

    a=T[0][0];

    b=T[0][1];

    c=T[1][0];

    d=T[1][1];

    tr=a+d;

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

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

    if (cabs(k)>1)

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

    else

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

    return(z);

}

 

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

{

    int i,j;

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

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

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

}

 

void inv(arry a, arry c)

{

    double _Complex d;

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

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

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

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

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

}

 

C-script Tracepqr

 

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

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

//To run successive compiled versions type ./depthf

 

#include <stdio.h>

#include <math.h>

#include <complex.h>

#include <time.h>

 

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

 

 

int main(void)

{

    //double _Complex val=2;

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

    double _Complex val=1.61691788-1.29437374i;

    int p=1597,q=2584,opt=0; //opt=0 free opt=2 non-free

    double ep=1e-8;

   

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

    int s[200],sl=0;

    double _Complex vh,h,f,df;

   

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

    h2=q;

    while (h2%h1>0)

    {

        hld=h2%h1;

        h2=h1;

       h1=hld;

    }

    g=h1;

    p=p/g;

    q=q/g;

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

    q0=1;

    p1=1;

    q1=0;

    ph=p0+p1;

    qh=q0+q1;

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

        if (ph*q>p*qh)

        {

            p1=ph;

            q1=qh;

            ph=ph+p0;

            qh=qh+q0;

            s[sl]=0;

            sl=sl+1;

        }

        else

        {

            p0=ph;

            q0=qh;

            ph=ph+p1;

            qh=qh+q1;

            s[sl]=1;

            sl=sl+1;

        }

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

    its=500;

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

    {

        vh=val;

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

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

        val=val-f/df;

        if (cabs(vh-val)<ep)

        {

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

            break;

        }

    }

    if (i==its)

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

}

 

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

    {

        double _Complex t0,t1,t2,h;

        int i;

        t0=val;

        if (opt<2)

            t1=val+2i;

        else

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

        t2=2;

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

        {

            if (s[i]==0)

            {

                h=t1*t0-t2;

                t2=t1;

                t1=h;

            }

            else

            {

                h=t1*t2-t0;

                t0=t1;

                t1=h;

            }

        }

        return(t1);

    }