Kleinian and Quasi-Fuchsian Limit Sets on Mac, in C and in Matlab

Kleinian and Quasi-Fuchsian Limit Sets: An Open Source Toolbox
Program Listings in Colour
Chris King

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

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

if lev==levmax

fl4=true;

end

if fl4

plt=true;

newpoint=points(1);

end

if lev==levmax

if sav>=2

fprintf('.');

dotfl=true;

dotc=dotc+1;

if dotc>=200

dotc=0;

fprintf('\n');

end

if sav==-1

outfl=false;

return;

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

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

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

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

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

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

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

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

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

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

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

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

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

}